Solving and Learning Partial Differential Equations with Variational Q-Exponential Processes
Authors: Guangting Yu
Shiwei Lan
Venue: NeurIPS 2025 poster ()
Keywords: Probabilistic PDE Solvers
Bayesian Inverse Problems
Data Inhomogeneity
Modeling Derivatives
Uncertainty Quantification
Primary area: probabilistic_methods
TL;DR: We propose a novel probabilistic solver for PDE and inverse problems based on a recent generalization of Gaussian process named Q-exponential process. The method yields more accurate solutions and provides meaningful uncertainty quantification.

Abstract: Solving and learning partial differential equations (PDEs) lies at the core of physics-informed machine learning. Traditional numerical methods, such as finite difference and finite element approaches, are rooted in domain-specific techniques and often lack scalability. Recent advances have introduced neural networks and Gaussian processes (GPs) as flexible tools for automating PDE solving and incorporating physical knowledge into learning frameworks. While GPs offer tractable predictive distributions and a principled probabilistic foundation, they may be suboptimal in capturing complex behaviors such as sharp transitions or non-smooth dynamics. To address this limitation, we propose the use of the Q-exponential process (Q-EP), a recently developed generalization of GPs designed to better handle data with abrupt changes and to more accurately model derivative information. We advocate for Q-EP as a superior alternative to GPs in solving PDEs and associated inverse problems. Leveraging sparse variational inference, our method enables principled uncertainty quantification -- a capability not naturally afforded by neural network-based approaches. Through a series of experiments, including the Eikonal equation, Burgers’ equation, and an inverse Darcy flow problem, we demonstrate that the variational Q-EP method consistently yields more accurate solutions while providing meaningful uncertainty estimates.

Solving and Learning Partial Differential Equations with Variational Q-Exponential Processes
Guangting Yu
Shiwei Lan
School of Mathematical & Statistical Sciences
Arizona State University, Tempe, AZ. 85287
Abstract
Solving and learning partial differential equations (PDEs) lies at the core of physics- informed machine learning. Traditional numerical methods, such finite differ- ence and finite element approaches, are rooted in domain-specifet techniques and often lack scalability. Recent advances have introduced neural networks and Gaus- sian processes (GPs) as flexible tools for automating PDE solving and incorporating physical knowledge into learning frameworks. While GPs offer tractable predictive distributions and a principled probabilistic foundation, they may be suboptimal in capturing complex behaviors such sharp transitions non-smooth dynam- ics. To address this limitation, we propose the use of the g-exponential process (Q-EP), recently developed generalization of GPs designed to better handle data with abrupt changes and to more accurately model derivative information. We advocate for Q-EP as superior alternative GPs in solving PDEs and associated inverse problems. Leveraging sparse variational inference. our method enables principled uncertainty quantification a capability not naturally afforded by neural network-basee approaches. Through series of experiments, including the Eikonal equation, Burgers' equation, and an inverse Darcy flow problem, we demonstrate that the variational Q-EP method consistently yields more accurate solutions while providing meaningful uncertainty estimates.
Keywords: Probabilistic PDE Solvers, Bayesian Inverse Problems, Data Inhomogeneity, Modeling Derivatives, Uncertainty Quantification
Introduction
It is of fundamental importance in science and technology solve mathematical models represented as system differential equations and to learn such a complex system by identifying crucial physical quantities with estimated uncertainty (a.k.a. inverse problems). Over centuries, theoretic foundations and computational methods have been developed for solving and learning partial differential equations (PDEs), which is facilitated by the development of modern computers. Traditional numerical algorithms such as the finite element method remain demanding for both domain knowledge and computing resources. Recently, there has been increasing interest and effort efficiently automate this process using machine learning techniques.
The surge of physics-informed machine learning is driven by two main thrusts: neural network-based algorithms and Gaussian process (GP)-based probabilistic solvers. The former works are represented by innfrse- neural networks [PINN 49), the deep Ritz method [9), the deep Galerkin
method [44, and operator learning (20 methods including the Fourier neural operator [FNO 25],
deep operator networks [DeepONet 26 and the neural inverse operator [NIO 28
See [9 for
slan @ @asu.edu
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

review of recent advances. The basic idea is to parametrize the solution with a neural network and to minimize certain loss with respect to network parameters to obtain the solution. Despite empirical successes, these neural network-type approaches typically require large samples but lack convergence guarantee or uncertainty quantification (UQ). On the other hand, GP has been introduced to solve and learn ordinary differential equations (ODEs) 45 39
4 43
[5] and PDEs [32 33 35| [8] with
theoretic guarantee [34] and UQ [6]. More recent development on GP for solving and learning
PDEs includes [6] 27 14 13
29 3 Different from neural network approaches, these probabilistic
solvers model the solution as a GP conditioned on PDE constraints and identify the solution as the maximum a posteriori (MAP).
Due to the tractability of conditional and predictive densities, GP has been widely adopted in machine learning and scientific computing [40]. However, as an L2 regularization, GP tends to be over-smooth for modeling certain objects with abrupt changes or sharp contrast. For example, it is known in imaging analysis that GP may not detect or preserve edges very well in an image [ 77. On the other hand, researchers [4 notice that total variation regularization penalizes the I,2 norm of derivatives and yields edge-preserving reconstructions. However, the total variation prior degenerates to GP prior with increasingly finer discretization mesh [23 and hence loses its edge-preserving feature. Therefore, [22] propose the Besov prior an Ly regularization and prove its discretization-
invariant property. (24] further develop the g-exponential process (Q-EP) as a probabilistic definition of the Besov process with tractable posterior prediction and demonstrate it a superior generalization of GP (with q = 2) in modeling inhomogeneous data with sharp transitions.
In this paper, we discover that Q-EP (with : =) GP, and hence presents as a preferable candidate for solving PDEs. Heuristically, this is attributed to Q-EP's enhanced ability to model inhomogeneous objects with sharp variations, resulting in better regularization of large derivatives. Theoretically, this can be explained by a faster posterior convergence rate in Bayesian modeling with Q-EP priors. Unlike optimization-based approaches 6 27), we adopt sparse variational inference [46 47 for Q-EP [30 5 to solve and learn PDEs, allowing natural UQ. An emerging challenge is that in addition to mapping the Q-EP mean function by the nonlinear PDE dynamics, one also needs to propagate the whole variational distribution through, which no longer renders Q-EP. We solve this difficulty by linearizing the complicated PDE mapping. We also extend the resulting variational Q-EP solver for inverse problems.
1) is better in modeling derivative information than
Connection the literature Our work is motivated by 66] which optimizes the log-posterior for MAP as the PDE solution. Our proposed method replaces GP with more general Q-EP and adopts variational Bayes for UQ We investigate Q-EPs in solving various forward and inverse PDEs for a spectrum of q with g = 2 corresponding to GP. As a probabilistic solver, Q-EP may not be best compared with neural network-based approaches. However, we still include PINN [37 and Bayesian PINN [B-PING PINN 49] as baselines in our comparison. We emphasize that our algorithms rely only on limited data, e.g. boundary values or interior observations, while providing meaningful UQ Our work is also related to the recently proposed physics-informa informed state-space GP [33, which however focuses on dependenn PDEs. It adopts a variational spatiotemporal state-space GP and can be regarded as related special case of ours for q = 2. Our work solving and learning PDEs makes multiple contributions to the field of physics-informed machine learning:
1. It is novel probabilistic PDE solver based on Q-EP with superior capability of modeling data inhomogeneity and derivative information.
2. It theoretically justifies the preference of Q-EP over GP in solving and learning PDEs.
3. It provides efficient UQ for solving forward and inverse PDE problems.
The remainder of the paper is organized as follows. Section 2 reviews Q-EP as flexible prior in Bayesian models for inhomogeneous data and introduces an extension to incorporate derivative information. Section [3 explains the details of applying Q-EP to solve the forward and inverse problems of PDEs. We follow [6] to model the solution as MAP of Q-EP but highlight the challenges of variational inference including distribution propagation and variational lower bound. In Section
we justify the preference of Q-EP for q = over GP (q 2) in solving PDEs. In Section 5 we demonstrate the numerical advantages, particularly faster convergence, of Q-EP compared with alternatives using forward problems involving Eikonal equation and Burgers' equation and inverse problems of identifying permeability in the Darcy flow. Section 6 concludes with discussion of limitations and future improvements.

---

2 Bayesian Modeling with Q-Exponential Process
2.1 Q-Exponential Process
The univariate q- g-exponential distribution 7 has density [u) x exp [/])), whose logarithm yields Ly regularization term. (24] generalize the univariate g-exponential random variable to multivariate random vector, based on which a stochastic process can be defined. Suppose a function u(x) is observed at N locations, "1.'.' XN Ω C TR. (24] define the multivariate q-exponential distribution for u == (u(xi),.... ... "((x)) )) as a member of the family of elliptic distributions [18]. Definition 1, A multivariate random vector E IR- follows the g-exponential distribution, denoted as u q- q-EDy (I, C), if it has the following density:
pl C, 2 2 ( (c) (1-1) exp {- r {-2 2 S
-u- (u
(I)
Remark 1. The negative log density of q-ED in (T) yields quantity dominated by some weighted --l From the optimization perspective, q-ED, when used as a prior, imposes Lgg regularization in obtaining the maximum a posteriori (MAP).
Lg norm of μe, i.e. 1r
Li et al. 24] prove that the above multivariate g-exponential distribution satisfies the conditions of Kolmogorov's 's extension theorem (3] and thus can be generalized to stochastic process. For this purpose, we scale u ~ q-ED v (0, C) by a factor N so that the scaled q-exponential random variable u == - tu q- yyE (0, C) has covariance asymptotically finite [Proposition 3. of 24). With a covariance (symmetric and positive-definite kernel C : S2 x + I, we define the following g-exponential process (Q-EP) based the scaled g-exponential distribution.
Definition 2. A (centered) q- q-exponential process u(x) with a kernel q-E C). is a collection of random variables such that any finite sel, g-exponential distribution q-ED' (0, C), where = [C(x;,x;)][<> x
(u(xi),- ..ux))))
=u(x))), follo)). scaled multivariate
Remark 2. When q = 2, g-ED C) reduces to NN C) and q-EP C) becomes GP(O, C). When q E (0, 2). q-EP(O, C) lends flexibility to modeling functional data with more regularization than GP. In practice, q is often adopted for faster posterior convergence 22 and the capability of preserving inhomogeneous features (rough functional data, edges in image, etc).
The covariance kernel C associated with a (22),, u(-) + Sq ((,,((((((())) (x)))(x'))(x())) which has eigen-pair (de, (bx(-))22 such that for VE N. Tx))(x) = dbr()) and ||øel12 = 1. Assume Te is trace-class, i.e. tr(Tc) == Ei= de < co. Theorem 3.4 of [24] presents a series representation of Q-EP similar to GP and the Besov process {.
(HS) integral operator Te: TI2 ->
Theorem 2.1 (Karhunen-Loéve), uf-) ~ q-EP(O, C) with a trace-class HS operator Te having eigen-pairs ()e, de(-)) then we have the following series representation for u(x):
u(x)
(x)) ue = S u(x))e((x)
ind q-ED (0,)
where E[ue] and Cov ug, ue Agdee with Dirac function dee if = l' and 0 otherwise. Moreover, we have Bi]lu(-) -==1 Elu = tr(Te) < co.
2.2 Bayesian Regression with Q-EP Priors
Given X [xx)) and yxx1 {yn J.Y we consider the Bayesian regression model:
y = u (X)+) + E (0,,E,
(3)
CP.
Li et al. 24 Theorem 3.5] show that the posterior (predictive) distribution is analytically tractable when both the prior and the likelihood are qeexponential.
Theorem 2.2. For the regression model
3 the posterior distribution of ul u(x.) at x, is
u( ((x,)) )ly, X, x, q- F " C.
y, C C.. (C+
where C = C), C, C(X, and C.. C(x.,X.).
3

---

35
32
14
65
3.0
25
3.5
3.2
15
1.5
16-
32
36
23
3.
20
21
20-
14-
Contrasting Q-EP (q- = 1.0, middle row) with GP (q- 2.0, bottom row) against the truth (top row) for modeling function values and derivatives of Rosenbrock (left) and Rastrigin (right).
Figure I:
2.3 Modeling with Derivative Information
Let u q-EP(O,CC. Denote the function and its derivatives by ä (u. ax to order k. Because linear operation preserves elliptic distributions [88 [0], ū ~ q-EP(O, is also Q-EP if C in Definition 2 is differentiable up to order k, where the augmented kernel, č, has structure illustrated in Table A.1 For example, the (1, block of Ć is interpreted as 0 Cov(u(x), dx u(x)) 7x xx- C(xxx) For solving the second-order PDEs in Section 5 we could adopt Matérn kernel (matern52 for 5/2), C(x, x') = o² +1+ v5r + exp( -vr).. r
2
8x u) up
V VELL (z. x/) pp for which the process is twice differentiable in the mean- mean-square sense.
To prepare for solving PDEs, we define || ||s.g for based on (2) with smoothness parameter
and an integrability parameter 7 2 22 7: [lu()) = (E (T(())y) ( ( (8) =
Consider the Banach space ( = {u Ω - IR [1|1(()))l < o... q = 2 and { J form
the Fourier basis, then (S2) reduces to the Sobolev space H (22). For u to be regular enough, we
make the assumption on č so that E (22) almost surely by the following proposition (21.
Assumption 1. Suppose X = { are eigenvalues of HS operator Tg for the kernel Ć. We assume X E E i.e. ||al Σi= Aa < xc.
Proposition 2.1. if)(.) ~ q-EP(O, C) with trace-class HS operator Té satisfying Assumption then ūi-) E ((()) = {ū R R R[ < oo) and El]((()) |ka|| <0.
Proof. See Appendix B
X)]. For brevity, we denote D = + kd. Then Y YNxa D and ENXD are the corresponding observations and errors respectively in the model (3). We apply this model to Rosenbrock (xx = Σ [100(2(41 - x²?² + (1 - =2)) and Rastrigin
Let Ú = ū(X) = [u(X), 8x X
g²
(f(x) = 10d + Diml" observed on a 20 x 20 grid, i.e. N = 400, = 2, 1. Figure contrasts Q-EP (q 1.0) and GP (q 2.0) in predicting function and derivative values the 50 x 50 grid. Q-EP outperforms GP in yielding a more accurate recovery, especially in the more challenging example of Rastrigin function.
10 cos(2rz;) test functions with function values and their derivatives
Heuristically, the superiority of Q-EP in modeling derivatives over GP comes from its improved ability to handle inhomogeneous data with sharp variation. From the perspective of the I. norm of the gradient, the L1 norm imposes stronger regularization than the Le norm on large values in V

---

resulting in better model for prediction. In Section 4 the following, we take advantage of Q-EP's capability of modeling derivative information and apply it to solve PDEs.
we will give more rigorous justification. In
3 Solving Partial Differential Equations with Q-EP
3.1 Bayesian Solver
Consider the following general PDE defined on a bounded domain Q C RR:
(x)) = f(x), (xu) g(x),
x E Ω
(4)
x E an.
where D (S)) - (S2) a differential operator and B: (32) (892) is boundary
operator with data (S2 and (aQ)). Here we assume sufficiently large smoothness index such that the PDE 4 well-defined pointwise and has unique strong solution [6] Let 0 12 02. For convenience of exposition, we denote the joint operator as P (D,
B* > (12), and the right-hand side function h (f.9) E L9()).
set of collocation points X {xn \. consists of Na interior points X {xI,... XN₂ Q} and N₂ boundary points X, {x (XN++I''' X. 02]. i.e. X X U Xs. and N = Na + No. Regarding the evaluation of Puu) on X, we make the following assumption so that we can properly define the likelihood model.
Assumption 2. There exists a differentiable function P: IR - R such that P(u)(x) = P(ū(x)).
And further there a constant C > O such that [| S C.
Then (X) becomes a nonlinear function of ū( XX) denoted P(ū(X)) = (X). Let h h()). The probabilistic solver seeks obtain U (X) based on observations (P(a (((((()))))))))))) h).
Even if we model (XX) (u, S) with ü, S to be specified (8) in Section 3.2 the nonlinear mapping P would not render ((ū(X) another q-ED random variable. Therefore, to properly define likelihood model, we propose the following distribution propagution by linearizing P:
(()( P(ūn) P((o (i)) q-FD(m,, I),
(5)
m = P(ūo) + P ( ūo), r = VP( + ÓlN.
where the Taylor expansion of P is about üg which can be chosen as ú, from the previous training epoch simply ú, and 8 >> O is a small nugget to ensure positive definiteness of r.
Let Y xx1 == (X)(.)) The Q-EP solver aims to minimize the discrepancy potential (negative log- lieliioe)))od) function, + ( x R - R, can then be defined:
Y h. The
(ū; Y) (r; (r: == (m- r- h).
(6)
where (r;; r, N) == - log r + x 1) logr- -r. Under Assumption continuous ü. which is used in the convergence theorem in Section 4
2 + is Lipschitz
Proposition 3.1. Suppose that the PDE mapping P satisfies Assumption 2 Let (0,2). Then for every r > 0, there exists L = L(r) > such that for every Y E EN and for all ūn ūg E (9)) with max [ 11: [(ا2lly]] <r, ( Y)- Φ(ū2; Y) < L u22ly..
Proof. See Appendix
B
Therefore, the Bayesian model for the solution u (4 can be summarized as
Y (Iū(X), h q- E> (h,
u CO.
Our goal is to infer the posterior p YÚ PY Note that because the extended function a(X) enters the likelihood model nonlinear way, Theorem 2.2 does not apply. the following, we solve the inference problem using variational Bayes
5

---

3.2 Variational Inference
We approximate the posterior p YU with some variational distribution (() using the varia- tional Bayes method, which aims to minimize the Kullback-Leibler divergence (U)) Because log P(Y) KL(q(U))P((()))) + (U)), it reduces to maximizing the lower bound L()(T)). The sparse variational approximation 46 47 is adopted by introducing inducing points X E pddx with their function values V ū(
RM
With the variational distribution for inducing values (VV) q-EDarD(u, E), the marginal varia- tional distribution (() can be obtained as 17 30
q( (Ū) - / Vi / M
~ q-ED(ü,,
(8)
= CNAN CMs vec(H),
S= CNN + CNM CNMC CM (E- ČMas) (MM))OM CM
The final evidence lower bound (ELBO) (g is (Refer to Section A.1 for more details.)
J (q(U) =p
log Σ (((vec(()) T CM + tr( (EECM CMM,
The variational solution (T) can be obtained by maximizing the ELBO (9 with respect to the variational parameters (n. Σ, X) and typer-parameters in the kernel C. By introducing the M inducing points, the overall computational complexity is reduced from O(N3) to ONNM2) 47
3.3 Bayesian Inverse Problems
The above Bayesian framework can be readily extended to solve inverse problems. The following adaptation enables us to obtain both forward and inverse PDE solutions simultaneously.
Suppose that the PDE (4) contains quantity of interest, a(x), which could appear in the differential equation D as part of the boundary condition B. The task of Bayesian inverse problems to find a true solution, at, with proper UQ based on observations. Suppose a is differentiable enough
gK aa) to the order k' < k. Now the joint operator P applies to
and we denote ä = (a, Xx
... u.... both u and a, which produces a nonlinear function of ū (). ((X) when evaluated on X, denoted as ((((X))(X)) P(u, (X).
In addition, there is an observation operator O such that observations, ( = u ((.),, are obtained on some set of No observation points, X, C with X., = No, This can be achieved by solving 4 with true at in simulations or simply modeling measurement data as noisy realization of (4) in real-world applications. Let Ñ = N N.. We supplement the joint equation operator P with the observation operator to form an augmented operator P = =P,O) O). Therefore, we have Y Xxx P(u, (X) = ((((X))(())
((() (( Similarly, we augment the right-hand side data h with observed u(X.) to make hixi [h u(X.)T!'.
If we model a(X) using q X)) q-ED(a, Sa)) with variational mean and covariance a, S. respectively, then (()(X(())X)))) propagates through the PDE dynamics similarly as in 5). Finally,
we summarize the Bayesian inverse model for a(x) as follows:
([(X), h q- ED &D (h,
f VV(u. an)
S.
VV(inn an) T + 5I (10)
S.
i~ i~(0 à ~ C.P(CC))
where åg can be similarly chosen as a. from the previous training epoch or simply a. The variational Bayes procedure in Section can be modified accordingly to obtain the variational solution of ā(X))Y. Y. Meanwhile, obtain the variational solution of ū ū(X) Y as byproduct.

---

4 Convergence Theorem
In this section, we study the posterior contraction of the Bayesian model (7) in the infinite data limit.
Similar theory can be developed for the model (1O). We focus on q E [1. 2 and leave the technically
more challenging case q E (0, 1) to future study. For brevity, we denote ü, C = C, and n N.
Consider the separable Banach space X = (L*(Ω), ll I1g)) D (B8))(1)) lle.g) for 8 > ( - 2).
and Y = RR. Define the concentration function of Q-EP measure IT at u = ut as
inf Put (e) be hEB*M(():] h- 11 2 hh 11.. log [T se).
(11)
Let P(P)) be the measure of the observations ( on (Y, B, po) having density Pa and corresponding potential function ((u;)) with respect to the Lebesgue measure No. i.e.
dr(n) (Y) = (Y) o duo
exp(-$(u: Y)). Define the Hellinger distance as dh,, (u, u') = [(VPu - -Pw) ddu. We have the following posterior contraction theorem.
Theorem 4.1 (Posterior Contraction). Let C-EP(O,C) with C satisfying Assumption | in
:= (9(()) and IP(") is the measure of ( parameterized by with PDE (4) satisfying Assumption
2 If the true value ul θ is in the support of u, and en satisfies the rate equation Pui (en) S ne
with en 2 - then there exists On C O such that TT E On : dr, (u, ut) ) 2 Mn Y(x) -
p(m) in at -probability for every Mn - co.
Proof. See Appendix B
O. By solving the inequality Put (En) & ne? for the minimal 6n, we obtain the posterior contraction rate as follows.
Denote a A b ma, b). a V b ma b). and x4
Theorem 4.2 (Contraction Rate). Let u ~ q-EP(O, C) with C satisfying Assumption 7in = (9()... The rest of the settings are the same as in Theorem 4.1 If the true value ut E B (S2) with sf >
s' +
d
d S.J and g¹ q E [1. 21, then we have the rate of the posterior contraction as
an sto where 08, St ,²))))))))))))))) = =(0 - --5
(
9 - )
4-41
En = n 2(oln.g.
Proof. See Appendix
B
Remark 3. we set the smoothness parameter =t d (4 and allow the integrability
parameter q S qt, then the contraction rate En is maximized as et, that this optimal rate is achieved regardless of the value of modeling regularization parameter q as long q q². This implies that when modeling inhomogeneous data and derivative information, under-smood smoothing (with smaller regularization parameter q) is preferred to over-smoothing. When the true integrability q² is at least Lı, setting q guarantees the fastest convergence of the posterior. In the following section, we will demonstrate such optimal choice with various numerical examples.
3. Note
5 Numerical Experiments
In this section, we test variational Q-EP for solving and Icarning PDEs. Using the Eikonal equation, Burgers' equation, a nonlinear elliptic equation (Section C.3 and an inverse problem involving the Darcy flow, we investigate the proposed method for a variety of different qs. We include PINN 37 and B-PINN 49 as baselines and the variational GP as a special case with 2 and mainly compare them using the relative error of the estimated solution â with reference to the true solution ut: RLE-p ||u | for p 1,2, sc. Q-EP with q outperforms PINN and variational GP in all ||=|||
cases and attains the best results in most comparisons. The numerical evidence also supports that Q-EP converges the fastest with the optimal choice of q 1.
7

---

SoMMinn
Uncemalrcy 1a
Uncemtrry
04
64
68
61 [pw]
620
Ci
618
03
618
33
Figure 2: Solving Eikonal equation 12) using high-resolution finite difference method (upper left), PINN
(lower left), Q-EPs (middle two upper row: q = 0.8; right two upper row: 4 = 1.0; middle two in lower row: q = 1.5), and GP (right two in lower row: = 2.0) respectively. Blue crosses are learned inducing points.
In most experiments, we choose the interior collocation points on 24 x mesh grid and the corresponding 100 boundary collocation points unless otherwise stated. For the sparse variational inference, M 256 inducing points are randomly initialized and learned by optimizing the ELBO ( The kernel C of Q-EP/GP is chosen to be matern52 with the hyperparameters, e.g. the correlation strength, automatically tuned in the python package GPyTorch [TT] implemented based on PyTorch. PINN is configured to have neural network parameters (weights and biases) of similar size to the variational Q-EP model. The computer codes are publicly available at https //gg::/gtttttt com/ lanzithinking/Diie
5.1
Eikonal Equation
First, we consider the following regularized Eikonal equation on f2 = [00, 1|2 also considered in 6
Vu(x)) <x)(x) = ff f(x)? x X E2 (xx) 0,
(12)
aΩ.
where / = and 0.1. Based on the set-up in Section 3.1 we define the nonlinear function
D(u, diu, dzu, dfu dzu) = (dlu)² + (dzu)² - E(dRu + d.)). Then the observations are Y =
((ū(X)) = Duu Xx), Des a (Xa), Jwz u (a), 0² u (Xa), 0ug u ((xa)) and h
[(x))² We ON.
(X))
obtain Na = 24 interior and N₂ = 100 boundary collocation points. Then we apply the variational Bayes in Section 3.2 to solve (2) with M 256 inducing points. Though not required for
convergence, we train each algorithm for 5000 iterations for fair comparison (PINN and B-PINN need much more training epochs than Q-EP solvers to converge).
Figure 2 compares solutions and uncertainty estimates generated by variety of Q-EP solvers for 0.8, 1.0, 1.5 and 2.0 (GP) respectively. We solve the equation using highly-resolver finite difference method with the Cole-Hopf transformation 6 and use as the true solution for comparison (upper left). We notice that the solution by PINN (lower left) is much worse and UQ is not available.
All Q-EP solvers (q < 2) produce better solutions than GP (q = 2) which does not precisely characterize the pyramid feature of the true solution. Only two solvers with q 0.8 and 1.0 yield solutions that correctly match the range of true solution. GP also manifests higher uncertainty in its generated solution. Table C.I further verifies the superior accuracy of Q-EP solvers compared to GP and PINN (B-PINN) in terms multiple error metries including MAE, MSE and RLE by repeating the experiments for 10 times with different random seeds. In this example, Q-EP solver with = 1.0 attains the result comparable to the most accurate solution. Note that with similar size of collocation points, our best results (4.43e-2 in L2 error and 1.03e-2 in Lxc error) comparable to those reported (1 (1.64e-2 in La error and 7.76e-2 in Lax error) in [6] from which UQ is absent.
5.2 Burgers' Equation
Next, we test our Q-EP solvers the Burgers' equation 6 with 0.1:
8

---

Peration 30 (e=1.0)
100 (q-1-00
beration 1000 lg-1.00
Reletive
+0.23
90c
-0.30
2.
10
1.0
Renddor 10 (umy
Heraton 100 (gm2200)
A.71
1.0
Bort 1000 fum 01
1.00
Pidulliv Error
616
610
A30
636
031
A11
640
-6.1
1026
1/0
-010
0.4
1.00
Meratioo (PINN)
190
Iorrtoon 1000 (PINNO
Rwative Error
6.-
0.43
10
0.8
64
04
-c.1
0.4
02
-0.38 12
4.5
0.
05
100
Figure 3: Comparing convergence of Q-EP (left three in top row: q= q= 1.0), GP (left three in middle
row: q = 2.0) and PINN (left three in bottom row) in solving Burgers' equation (13), with the right column illustrating the error reducing in L1 norm (top). L2 norm (middle), and Lo norm (bottom) respectively. Blue crosses are learned inducing points. Shaded regions are standard errors based on 10 repeated experiments.
0 ØL
82
u 0, (I, E (-1, x (0,1), O1
Ox
(13)
u(x, 0) sin(T.)), x (-,,1),
u(-1,1) - -(1,()=1))) / (0.
We use the same experiment setup as above. Figure compares the solutions by a highly-resolved finite difference method (upper left, treated as true solution for comparison purpose). PINN (lower left), Q-EP (right three in upper row: g = 1.0), and GP (right three in lower row: q = 2.0). Q-EP about one order of magnitude more accurate than PINN and GP, which can be verified from the pointwise error plots and Table [ The plots of posterior standard deviation on the right column also indicate meaningful uncertainty in the middle area around x = where the shock is difficult to resolve. Q-EP still achieves one order of magnitude lower uncertainty compared with GP.
Table I: Comparing accuracy of various solvers for Burgers' equation 13 in terms of mean absolute error (MAE), mean squared error (MSE), and relative errors in 71 norm (RLE-I), ), L2 norm (RLF-2), and I norm (RIF-oo) respectively. Result in each cell are averaged over 10 experiments with different random seeds; values after + standard deviations of these repeated experiments.
Model (g)
MAE
MSE
RLE-1
RLE-2
RL-0c
PINN
5.81e-2 + 1.03c-4
7.12c-3 + 4.963-3 3.96c44
0.1508 + 0.0017
0.1896 + 0.0052
0,2842 = 0.0071
B-PINN
2.94e-2 :: 1.57e-2
1.67e-3 1.58e-3
0.0785 + 0.0420
0.0833
0.0409
0.1327 = 0.0446
0.5
7.77e-2 4: 7.099e-
3.50e-2 + 7.15e-2
0.2018 0.3751
0.2056 0.3764
0.2177 0.3838
0.8
1.66c-244.966
0.0405
0.0780
0,0434
0.0804
0.0529 = 0.0924
9.33e-3 1: 1.76
1.77e- 2.16e-4
0.0242 0.0128
0.0266
1 0.0140
0.0324 - 0.0145
1.87-22+
5.92c-4 4.06e-4
0.0485 0.0151
0.0522 0.0168
0.0560 0.0181
1.5
2,64e-2 8.58e-4
1.31e-3 + 1.23e-3
0.0684 0.0262
0.0754 0.0316
0.0854 0.0430
2.0 2.uussa))
7.13e-2 8.68--3
1.08c-2 # 1.28c-2
0.1848 0.0908
0.2068 0.1118
0.2376
0.1456
2.5
2.26e-1 3.20e-1
2.49e-1 1 6.04e-1
0.5879 0.8296
0.6669 0.9576
0.7558 1.0521
Table compares comparess the accuracy of solutions in terms of MAE, MSE and RLE, for which Q-EP with 1.0 attains the most accurate one, verifying its best performance. The error for Q-EP models does not decrease monotonically with q > 0, but reaches the lowest at q 1.0. To further investigate the convergence, we plot solutions by Q-EP (top row: q = 1.0), GP (middle row: q 2.0), and PINN (bottom row) in Figure 3 at (first column), 100 (second column), and 1000 (third column) iterations respectively. We can tell that Q-EP with g = 1.0 converges the fastest to the best estimate, which confirmed by the error-reducing plots on the right column. Note that in this example, our algorithm does not need random perturbation of mesh grid as required by (6].

---

Irveme
Invene (MNN)
04
04
24
04
6.3
C4
64 C4
C2 Franee to
13
Koe
ICe Sulibon O9 [N)
6:0 verria 11-17
08
com
cme
M0
con
300
6000
6.00
Figure 4: Solving inverse Darcy flow (14) using PINN (second column), Q-EP (right three in upper row: q 1.0), and GP (right three in lower row: y = 2.0) respectively. Upper left: true inverse solution at; lower left: fine-resolution finite element solution u1 to (14) with a². Blue crosses are learned inducing points, and red dots indicate locations of observations.
5.3 Inverse Darcy Flow
Now we consider an inverse problem that involves the following Darcy flow:
(di((exp())Vu)G = f(x), x s == [0, 1]2,
(14)
u(x) 0, x E On.
where the true coefficient such that exp(al (x)) = exp(sin( 2 +in (2Tx2)) ) +exp(- sin (2xx1)
sin(2xrr2) [6] is plotted in the upper left panel of Figure 4 We generate data by solving (14 with a on a mesh 80 x 80 to obtain ul using the finite element method, illustrated in the lower left panel in Figure 4 On a (coarser) mesh N₂ = 20 x 20 used for inference, we randomly select N, 100 points X. Q and obtain observations as u! (X.) + E with noise N N(0, (0,,721g 22IN.) for 10-3.
We solve the inverse problem finding in 14 given these observations. With N₂ 400 interior and N₂ = 84 boundary collocation points and M 256 inducing points, we train Q-EP solvers for 2000 iterations. F Iillustrates the forward PDE solution (third column), the inverse solution (forth column), and the uncertainty of (rightmost column). Compared with the true at, Q-HP (q = 1.0) recovers a more faithfully than GP (g 2.0) and PINN (upper in the second column). Meanwhile, Q-EP also generates the solution much closer to that by the finite element method, as shown in the lower left panel. Higher uncertainty is observed by Q-EP (q = 1.0) around the corners with less data (Figure C.3p. reflecting the configuration of observations. If run for longer (e.g. 5000) iterations, PINN may improve its forward solution, but still gets a poor inverse solution (Figure C.4).
Table C.5 compares the relative errors (RLEs) of forward and inverse solutions. Q-EP with 1.0 also achieves the best or comparable solutions. Here we emphasize that all the results are based on only No = 100 observations from one solution, as opposed to the thousands of PDE solutions typically used in operator learning algorithms. Figure compares the inverse solutions of Q-EP with those of GP in increasingly finer meshes. If we view training in finer mesh as process to see more data, Q-EP (g = 1.0) has already converged fast to better estimates in coarser mesh, leaving less room to improve compared to GP (q = 2.0).
6 Conclusion
In this paper, we propose variational Q-EP to solve and learn PDEs. We advocate Q-EP with q = 1.0 over GP (q 2.0) for modeling derivative information and hence better probabilistic solver for PDEs. The fastest convergence at 1.0 is theoretically justified and empirically verified using two nonlinear forward PDE problems and an inverse Darcy flow problem.
One of the limitations might be the variational inference adopted this paper. The highly nonlinear nature of some PDEs imposes challenges on the quality of variational approximation to the resulting posterior, which in turn may undermine both the solution accuracy and the associated UQ. A possible remedy could be more flexible inference methods, such as normalizing flow [4] 36).
10

---

References
[1] Sergios Agapiou, Masoumeh Dashti, and Tapio Helin. Rates of contraction of posterior distributions based on p-exponential priors. Bernoulli, 27(3)1616 1642, 2021.
[2] Frank Aurzada. On the lower tail probabilities of some random sequences in lp. Journal of Theoretical Probability, 20(4):843-858 Dec 2007.
[3] Ricardo Baptista, Edoardo Calvello, Matthieu Darcy, Houman Owhadi, Andrew M. Stuart, and Xianjin Yang. Solving roughly forced nonlinear pdes via misspecified kernel methods and neural networks. arXiv preprint arXiv:2501.17110 01 2025.
[4] Ben Calderhead, Mark Girolami, and Neil Lawrence. Accelerating bayesian inference over nonlinear differential equations with gaussian processes. In D. Koller, D. Schuurmans, Y. Bengio, and L. Bottou, editors, Advances in Neural Information Processing Systems, volume 21. Curran Associates, Inc., 2008
[5] Zhi Chang, Chukwudi Obite, Shuang Zhou, and Shiwei Lan. Deep q-exponential processes. In Proceedings of the 7th Symposium Advances in Approximate Bayesian Inference. AABI, 10 2025.
[6] Yifan Chen, Bamdad Hosseini, Houman Owhadi, and Andrew M. Stuart. Solving and learning nonlinear pdes with gaussian processes. Journal Computational Physics, 447:1 10668, 2021.
[7] Masoumeh Dashti, Stephen Harris, and Andrew Stuart. Besov priors for bayesian inverse problems. Inverse Problems and Imaging, 6(2):183-200 may 2012.
[8] Masoumeh Dashti and Andrew M. Stuart. The Bayesian Approach to Inverse Problems. In Roger Ghanem, David Higdon, and Houman Owhadi, editors, Handbook of Uncertainty Quantificatin, pages 311 1-148 Springer International Publishing, Cham, 2017.
[9] Weinan E and Bing Yu. The deep ritz method: A deep learning-bassee numerical algorithm for solving variational problems. Communications in Mathematics and Statistics, 6(1):1-12, Mar 2018.
[10] K. Fang and Y.T. Zhang. Generalized Multivariate Analysis. Science Press, 1990.
[11] Jacob Gardner, Geoff Pleiss, Kilian Q Weinherger, David Bindel, and Andrew G Wilson. Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration. In Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
[12] S. Ghosal and A.W. van der Vaart. Convergence rates of posterior distributions for non-i.i.d. observations. Annals of Statistics, 35(1):192 2007. MR2332274.
[13] Oliver Hamelijnck, Arno Solin, and Theodoros Damoulas. Physics-informed variational state- space gaussian processes. In The Thirty-eighth Annual Conference on Neural Information Processing Systems, 2024.
[14] Marc Harkonen, Markus nongg-HHgrmmann. and Bogdan Raita. Gaussian process priors for systems of linear partial differential equations with constant coefficients. In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett, editors, Proceedings of the 40th International Conference Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 12587-12615. PMLR, 23-29 Jul 2023.
[15] Markus Heinonen, Cagatay Yildiz, Henrik Mannerström, Jukka Intosalmi, and Harri Lähdesmäki. Learning unknown ODE models with Gaussian processes. In Jennifer Dy and Andreas Krause, editors, Proceedings of the 35th International Conference on Machine Learning. volume 80 of Proceedings of Machine Learning Research, pages 1959-1968. PMLR. 10-15 Jul 2018.
[16) Philipp Hennig, Michael A. Osborne, and Mark Girolami. Probabilistic numerics and uncertainty computations. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471121179):015504422 July 2015.

---

[17] James Hensman, Alexander Matthews, and Zoubin Ghahramani. Scalable Variational Gaussian Process Classification. In Guy Lebanon and S. V. N. Vishwanathan, editors. Proceedings of the Eighteenth International Conference on Artificial Intelligence and Statistics, volume 38 of Proceedings of Machine Learning Research, pages 351-360, San Diego, California, USA, 09-12 May 2015. PMLR.
[18] Mark E. Johnson. Multivariate Statistical Simulation, chapter 6 Elliptically Contoured Distribu- tions, pages 106-124, Probability and Statistics. John Wiley & Sons, Ltd, 1987.
[19] George Em Karniadakis, Ioannis G. Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informan machine learning. Nature Reviews Physics, 3(6):422 -40,22244042. May 2021.
[20] Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications to pdes, Journal of Machine Learning Research, 24(89):1-97, 2023.
[21] Shiwei Lan, Mirjeta Pasha, Shuyi Li, and Weining Shen. Spatiotemporal besov priors for bayesian inverse problems. arXiv preprint arXiv:2306. 16378, 06 2023.
[22] Matti Lassas, Eero Saksman, and Samuli Siltanen. Discretizattion-inaaiin bayesian inversion and besov space priors. Inverse Problems and Imaging, 3(1):87-122, 2009.
[23] Matti Lassas and Samuli Siltanen. Can one use total variation prior for edge-preserving Bayesian inversion? Inverse Problems, 20(5):1537, 2004.
[24] Shuyi Li, Michael O' Connor, and Shiwei Lan. Bayesian learning via q-exponential proees. In
A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, editors, Proceedings of the 37th Conference on Neural Information Processing Systems, volume 36, pages 72867-72887. Curran Associates, Inc., 2023.
[25] Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadeneshell, Burigede liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations, 2021.
[26] Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nature Machine Inielligence, 3(3):218-2) 8-229, March 2021.
[27] Rui Meng and Xianjin Yang. Sparse gaussian processes for solving nonlinear pdes. Journal of Computational Physics, 490: 1000:123:0 12340, 2023.
[28) Roberto Molinaro, Yunan Yang, Björn Engquist, and Siddhartha Mishra. Neural inverse operators for solving PDE inverse problems. In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato. and Jonathan Scarlett, editors, Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 25 25105-25133 1399. PMLR, 23-29 Jul 2023.
[29] Carlos Mora, Amin Yousefpour, Shirin Hosseinmardi, and Ramin Bostanabad. A gaussian pro- cess framework for solving forward and inverse problems involving nonlinear partial differential equations. Computational Mechanics, (14))12213-1333 2025.
[30] Chukwudi Paul Obite, Zhi Chang, Keyan Wu, and Shiwei Lan. Bayesian regularization of latent representation. In The Thirteenth International Conference on Learning Representations, 2025.
[31] Bernt Øksendal. Stochastic Differential Equations. Springer Berlin Heidelberg, 2003.
[32] Houman Owhadi. Bayesian numerical homogenization, Multiscale Modeling and Simulation, 13(3)=812-828, January 2015.
[33] Houman Owhadi. Multigrid with rough coefficients and multiresolution operator decomposition from hicrarchical information games. SIAM Review, 59(1)999-149 January 2017.
12

---

[34] Houman Owhadi and Clint Scovel. Operator-Adapted Wavelets, Fast Solvers, and Numerical Homogenization: From a Game Theoretic Approach to Numerical Approximation and Algorithm Design. Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, 2019.
[35] Houman Owhadi and Lei Zhang. Gamblets for opening the complexity-bottlened of im- plicit schemes for hyperbolic and parabolic odes/pdes with rough coefficients. Journal of Computational Physics, 347:99-128, 2017.
[36] George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan. Normalizing flows for probabilistic modeling and inference. Journal of Machine Learning Research, 22(57):1-64, 2021.
[37] M. Raissi, P. Perdikaris, and G.E. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
[38] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Numerical gaussian processes for time-dependent and nonlinear partial differential equations. SIAM Journal on Scientific Computing, 40(1 40(1):A172 11722--1998 January 2018.
[39] O. Ramsay, G. Hooker, D. Campbell, and J. Cao. Parameter estimation for differential equations: a generalized smoothing approach. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69(5):7417766, 2007.
[40] Carl Edward Rasmussen and Christopher K. I. Williams, Gaussian Processes for Machine Learning. The MIT Press, 2005.
[41] Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In Francis Bach and David Blei, editors, Proceedings of the 32nd International Conference on Machine Learning. volume 37 of Proceedings of Machine Learning Research, pages 1530-1538, Lille, France, 07-09 Jul 2015. PMLR.
[42] Leonid I. Rudin, Stanley Osher, and Emad Fatemi. Nonlinear total variation based noise removal algorithms. Physica D: Nonlinear Phenomena, 60(1):259 2285, 1992.
[43] Michael Schober, David Duvenaud, and Philipp Hennig. Probabilistic ode solvers with runge- kutta means. In Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K.Q. Weinberger, editors, Advances in Neural Information Processing Systems, volume 27. Curran Associates, Inc., 2014.
[44] Justin Sirignano and Konstantinos Spiliopoulos. Dgm: A deep learning algorithm for solving partial differential equations. Journal of Computational Physics, 375: 339-1364, 2018
[45] John Skilling. Bayesian solution of ordinary differential equations. In C. Ray Smith, Gary J. Erickson, and Paul O. Neudorfer, editors, Maximum Entropy and Bayesian Methods: Seattle, 1991, pages 23-37. Springer Netherlands. Dordrecht, 1992.
[46] Michalis Titsias. Variational learning of inducing variables in sparse gaussian processes. In David van Dyk and Max Welling, editors, Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics, volume 5 of Proceedings of Machine Learning Research, pages 567-574, Hilton Clearwater Beach Resort, Clearwater Beach, Florida USA, 16-18 Apr 2009. PMLR.
[47] Michalis Titsias and Neil D. Lawrence. Bayesian gaussian process latent variable model. In Yee Whye Teh and Mike Titterington. editors, Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statisties, volume of Proceedings of Machine Learning Research, pages 844-851, Chia Laguna Resort, Sardinia, Italy, 13-15 May 2010. PMLR.
[48] C.R. Vogel and M.E. Oman. Fast, robust total variation- based reconstruction of noisy, blurred images. IEEE Transactions on Image Processing, 7(6):813 824, 1998.
[49] Liu Yang, Xuhui Meng, and George Em Karniadakis. B-pinns: Bayesian physics-informed neural networks for forward and inverse pde problems with noisy data. Journal of Computational Physics, 109911, 2021.
13

---

Technical Appendices for "Solving and Learning Partial Differential Equations with Variational Q-Exponential Process"
A Calculations
Table A.I: The structure of kernel Ć with derivatives.
C-v(.,-)
u(x')
@ u(x'
8(x') u(x')
8x
u(x)
C(x, x')
a C(x, x')
C(x, x')
0x
8x ur(x xx
C(x, x')
C(x, x')
pxprzey C(x, x')
8x
5x2 u(x)
0x2 C(x,x')
0x20x x(x,x')
0a
0x20(x'22 C(x, x')
A.1 Variational Lower Bound
Because PYY (L(((ū)IP(U)))) (()()), it reduces to maximizing the lower bound C(q(U)). The sparse variational approximation [46] 47 is adopted by introducing inducing points E R MxD Hence the joint distribution of and
X E RMxd xd with their function values V úXX)
U can be augmented by including V:
p Y, Ū) x p
x, X.). (V)))
where we have vec( [ MD (0, mMM) and the conditional distribution
vec(Ú)] x, ND Cy (V), ČNN CNMC CMN).
Now we approximate the joint posterior p(ū, Y with the following variational distribution
(Ū V) - ((V), q(V) g-ED MD (M, E),
where the covariance of size MD MD and can be chosen as a (block) -diagoal matrix for convenience. A standard variational Bayes procedure yields the variational bound:
logp(Y)> > / ((ū, V) log (( X)P(V VV abaV
(Ú.
/ ((V)) log P(V)
VVV
dV
((V)
= (Ú) log (YJŪ) - (L(G(VV
where the marginal variational distribution q( ÚÚ) can be obtained as 17 30
(()) ((O VaVV M(UMV ~ Q-ED(ū,S)
û ú= CNCM C vec(H),
S = ČNN + CNMO - CMM) MM CMN-
Denote by (r(; C) = - log CI + N [ 1) logr 1 which is convex for (0, 2]. Let spo(r) = ((;; P, N) and r(Y) = ( h) r-- - h) quadratic form of random variable Y. Then log logu(Y\() = :p((r(Y)). Therefore, by Jensen's incquality, we can bound from below as
(log p( YY (Y())) = (Po(r(Y))) 2 (ol((r())) (r(Y)) = - tr(
Now compute the K-L divergence KL (L(I(V)
KL - logg - - log = (Vg(V) - -llo (V))(())
14

---

Denote by r(V) vecl V - n))
T Iー veeve (.). Then log v(r(V); Σ, MD). From
[Proposition A.1. of 24 we know that rř (MDD). Therefore
Hz (V) - log Σ
MD
q
-1) (1x H(22(M + MD
2
1
MD
2
MD
MD
(1
MD
MD
MD
= log |1 +
1
2 log 2r
4
2
2
2
Let (r)))) = sr;;; CMr,, DD. Then by Jensen's 's inequality
(log (V)) vV; = (421 (V) T MMs voc (V))) (() 2 PI((voc(V) (vee(V) T O MM vec v) (vec(u) C tr (r(ECwm)
T " "MM vec V)) ((v))
Therefore, the final evidence lower bound (ELBO) (q())) is
logp(E) L (q(
= u tr(T="S):I,,) I)
+ log Σ + (((vec(H) TO vec(H)) + tr tt(EErt CMM, MD). 2
B Proofs
Notations: S means "less than or approximately equal to"; on S ba implies an < Cbb, for some
constant C > means "asymptotically equal to"; on x bn implies 1im
On on C for some
constant c.
Proof of Proposition 2.1
Note r(ūs) = X W (1) for all E E N by Proposition A. of std 24]] Denote x² == x @ ~
A X becomes an infinite mixture of chi-squared random variables whose density is analytically intractable. Yet we have
(1). Hence | Σ
E B
Σ
Σ>
x
<0,
Assumption Taaiooooooooo holds. Thus it completes the proof.
Proof of Proposition 3. Based on 5 and 6 the potential Φ(ū; Y) = -p(r; [,r;; N) where
<p(;: r, N) = log r + (2 log r r is convex in r if g (0, and = r(ū)
with
VP(.)), b -(ū.) + u th, r ASAT T +8I. Since convex function is Lipschitz continuous over compact domain, it suffices to prove that r(ū) is bounded (both bounds achievable) and Lipschitz.
r(ū) -) r- (Aū-b)
Note that rūū) S A bl2 S ([ ||| [[b]]) where represents PDE solution the
compact domain 52 2 and is therefore bounded by some M > By Assumption 2 ||| S C. Therefore, S rū) < 8- 8-1 (CM + ||b]))
On the other hand, we have the gradient of the quadratic form r(ū) bounded
|vr(ū]] = 2||A r- - b)ll < | bl s 28- + [[b]).
Hence r(ū) is Lipschitz.
Lastly, there exists L,, La > 0 such that
[0; Y) (ū Y) = (p(r((r): ((((()): T,N)
< [[r((ū)) - r(ū)) S 1 T- ū1 u21l4
1 | أل ll
2L
[
assuming the collocation points are uniformly sampled. Therefore, it completes the proof.
15

---

According to Theorem 3.1 and Lemma 5.14 of []. the following general contraction conditions hold for Q-EP and will be used in the proof of posterior contraction Theorem
Theorem B.1. Let μ be q-EP(0, C) measure satisfying Assumption | in the separable Banach space (I.)))), || llq). where q E [1. Let a ~ and the true parameter ut (.t).. Assume
En > O such that Pue (En) S ne, where ne Z 1. Then for any C > 1, there exists measurable
set Bn C (S2) and a constant R > 0 depending on C and q, such that
log (4en: Ba-11-1la) 1 Rne
(15)
1uuu & B.))
5 exp(-Cneal)
(16)
p 11lu-wlll 1 2
2 exp(-nen))
(17)
where (4e,, BR: 11 (1.)) is the minimal number of bly-balls of radius 4e, to cover B
We need the following lemma bound the Hellinger distance, Kullback-Leibler (K-L.) divergence, and K-L variation to complete the proof of Theorem 4.1 21
Lemma B.I. Suppose the potential function Then we have
(6) satisfies Lipschitz continuity in u as Proposition
3.
dı (Pas Pu') S ||0- 11a
"K(pu I ||u ut 11
VV(u.Pw') S I|u - u' v'l1a-
Proof. First, consider K-L divergence:
K Pw) = / log Pi; Pw
Pu Pudp / (d((u's)) ((rr))),duuv) S Ll| - u u'lle
by Proposition 3.1 Similarly, we have for K-L variation:
V( Par) / (1og Pu
Pu' ) Pudd (a(ull ((u:)) L2 u 11....
Lastly, we bound the Hellinger distance:
f (Vu)- / [1-eex -ep exp ( f (20(1:)) (ud' )}² pudu(y)
2df (Pus par)
/ C 4
CL
v) -0( ((u
1|u - u' 11
where the inequality holds for I|u w' 11 small enough.
Proof of Theorem 4.7 Based on [Theorem [f suffices to verify the following two conditions
(the entropy condition (2.4), and the prior mass condition (2.5)) for some universal constants 77, K > 0
and sufficiently large k E N.
sup log N (TJ/22 {u On : dr, (u, ut) < dn.t)) E>en
<
ne
(18)
(u E On ken < dr, (u, ul, I 2kEn)
< Kne2,2112
(19)
II (B, (ut,en))
where the left side (18) is logarithm the minimal number of dh, -balls of radius <e/2 needed to cover ball of radius around the true value ut; B, (ut, .n)) = {u E _ u) s
: V.V u) < {2) with ((ut, u) = K ((put. Pw) and V(ut, u) V(put, pw).
Since u(-) E (S2) satisfy conditions for Theorem B.I there exists B, C (*( such that 15)
(17) holds. Now we set θ B.. For Vu, w' E O such that [|u(-) - u' u'( S Ens we have (n, w') X ||u u'll, S &. by Lemma B. Therefore by [15 we have the following global entropy bound holds
log N O dn, Rne,
16

---

which is stronger than the local entropy condition 18]
Now by Lemma B.T and 177. have
II (L, ..,) (n ullg S E ||u|
ull3 S (Ilu =ull - nch
Then the prior mass condition 19 is satisfied because the numerator is bounded by The proof is hence completed.
The following lemma studies the small ball probability in the concentration function
Lemma B.2 (Small ball probability). Let 11 be a q-EP(O, C) prior on B (S2) with < 8 Then as E - O, we have
log [[[[u]]ass S 5c) xe -
Proof. We can compute
[[ = =P
Σ ( wel) s
where P is the probability measure on the infinite product space (L9())) From the proof of Proposition 2.1 we know ||u]| Ee > x² is an infinite mixture of x2) random variables, so the condition [Theorem 4.2 of is trivially met and we have
8 Σ (g7(6)))7)0)))
log P
11 S
xe 7000)
The second lemma gives upper bound of the first term of the concentration function
21
Lemma B.3 (Decentering function), Assume u² Bs (52) for some > s' and qt E [1, 2]. Then as + O, have the following bounds
(i) If 2 q. require sf s':
iʃs < st = if
loge) 1- J. st (angt)
inf (Q): w
|
(ii) If < q. require s s d q.
1
if + +
inf 2. (EB*
|/ I
9
if 8l + 4
Proof. We identify ut E Bs oo' with fu/ e (t))r Then we follow I to approximate
ul with hI. {u] where uf 0 for all l > L. Note his (1,7(()) for any finite L N. Identifying hTr with (R), we could get
||utt1g] st.
1, -'=
isq =
|/ - -uull 1111
Σ
)' -'
if >q.
E=
( ++))
| gt
if <1

---

Therefore, to have ||/ ut llor. S E we let
ifal 2 q
I こ
(20)
8
if <9
On the other hand, the infimum is less than ||A]],,, with above h, which can be bounded as follows. If
1¹ 9.
| = E en Juglo' 1
|1u oloptt || llg st.gg L 1
ifs S st
if 8 s
If q² q, by similar argument using Hölder inequality,
Cllutll gt
ifs <
| =
grolo
|| (log L) -
if ==
|| 11 st.gt L
i58> 8
If q" q, by similar argument,
| 11 at.
iss s st + d
31
<
1=1
|
at
ifs> st +
Substituting L in 20
to the above equations yields the conclusion.
Proof of Theorem
4.2 By Lemmas B.2 and
B.3
we have the following bounds for the concentration
function
as +
0, q² 2 q.
+
4
if st
Put (
loge) 1- St
+
m-
if st
(ana")
+
- if8> st
For S st, the bound is dominated by - For the last case, need to determine balancing point of 8 for the two terms by setting their powers equal. The calculation shows that if 8 S sf + d
9'
but otherwise is dominated by 3 Therefore,
the bound is still dominated by have
- ifs S
Pat S
T-2"
ifs> 8² d
We need to determine minimal en such that (ut
S ne Hence for qt 2 g.
w's-d/ied n 111
ifs 5 +
En
n 2ل (
ifs>st+
Now if < by similar argument we have the concentration function as S + 0
I+e E
iʃa 58
( S
اييع -4-
-واد t
+e
- if > st + d d
18

---

Thus the contraction rate for q¹ < becomes
71 Zul () 2)d
2d d if S
En 2
26a 141s
s1 2)
ifs> s
24 d
n
Rewriting the equations into one yields the conclusion.
C More Numerical Results
C.1 Eikonal Equation
Table C.I:
Comparing accuracy of various solvers for Eikonal equation 122 in terms of mean absolute (MAE), mean squared error (MSE), and relative errors L1 norm (RLE-1), L2 norm (RLE-2), and Loe norm (RLE-oo) respectively. Result in each cell are averaged over 10 experiments with different random seeds; values after + are standard deviations of these repeated experiments.
Model (g)
MAE
MSE
RLE-1
RLE-2
RL-0c
PINN B-PINN
1.77e-2 :: 9.03c-4
4.65e-4 4.08c-5
0.1120 0.0057
0.1167
0.0051
0.1354 + 0.0061
1.81e 1: 2.59e-2
1.24e:3 3.52e-3
0.1253 0.1792
0.1 0.1700
0.2143 1 0.2193
0.5
5.09c-
6,20c-6 + 2.87c-6
0.0124 0.0032
0.0132 0.0029
0.0387 + 0.0067
0.8
3.50ee- + 1.28e-6
0.0091 + 0.0022
0.0100 + 0.0018
0.0339 + 0.0074
1.0
4.42e-6 = 2.26c-6
0.0106 + 0.0034
0.0110 = 0.0030
0,0321 = 0.0071
1,2
3.835
4.11e-6 1. 1.52e-6
0.0107 T 0.0024
0.0108 1 0.0020
0.0284 1 0.0057
1.5
1.68--3 3.31ec4
4.65e-6 + 1.35e-6
0.0106 0.0021
0.0116 + 0.0018
0.0398 + 0.0080
2.0 (Gaussian)
9.64e-3 1.33e-3
1.30c-4 +: 3.61e-5
0.0610 + 0.0084
0.0612 = 0.0087
0.1009 + 0,0124
2.5
7.71e-2 4: 3.18e-2
8.36e-3 1: 5.29e-3
0.4881 1 0.2011
0.4668 0.1743
0.4442 + 0.1
To verify the systematic superiority of Q-EP with = 1.0 than GP (q = 2.0), we extend our comparison to more kernels including the radius basis function (rbf), C(x, x') 0² exp( 00.5r2)
and rational quadratic (rq). C(x, x') = ( + 2a 2 a 0, for r v (x (x: x1)2 in Table Within each type of kernel, Q-EP with = 1.0 consistently outperforms GP (q = 2.0). This indicates that the advantage of Q-EP over GP is independent of the choice of kernels (as long as they are compared with the same kernel).
Table C.2: Comparing accuracy of Q-EP (q against GP (q 2) solvers with various kernels for Eikonal equation [22 terms mean absolute error (MAE), mean squared error (MSE), and relative emrors L1 norm (RLE-I), La norm (RLE-2), and Los norm (RLE-oc) respectively. Result in each cell are averaged over 10 experiments with different random seeds; values after + are standard deviations of these repeated experiments.
Model (g)
kemel
MAE
MSE
RLE-1
RLE-2
RL-20
1.0
Matern
1.68c3 = 5.41e44
4.42c-6 2.26c-6
0.0106
0.0034
0.0110 = 0.0030
0.0321 + 0.0071
2.0 (Ganssiun)
Matern
9.64e- + 1,33e-3
1.30e 3.61e-5
0.0610 0.0084
0.0612
0.0087
0.1009
0.0124
1.0
rbf
4.39e-3 1.05e-3
3.91e-5
2.52e-5
0.0278
0.0066
0.0327
0.0107
0.0648 0.0324
2.0 (Gaussiun)
rbf
1.55e-2 = 4.97c-
3,07e- +: 2.87e-5
0.0982 = 0.0031
0.0949 = 0.0044
0.1236 + 0.0048
1.0
rq
186c-3 + 5.40e44
5.51e-6 3.44e-6
0.0118 : 0.0034
0.0122 + 0.0038
0.0168
1: 0.0064
2.0 (Gaussian)
rq
1.70e-3
1.89e-5 2.83e-5
0.0194
0.0200
0.0134
0.0325
C.2 Burgers' Equation
Define the nonlinear function D(u, di d2u, dzu,dRu) dfu) dru + ud, udfu. The observations Y can
then be expressed as P(ū(X)) = ((u(X.), ar u( Xx), au Bu (X.a):
uz ul (Xa)) and h 0
Ø²
w (Xx))
where g is a vector of size N₂ whose elements are sin(TT) or 0 depending on the order of the corresponding elements in Xs. We use the same experiment setup as above.
In Table C.3, we also observe a consistently better performance Q-EP (g = 1.0) compared to GP (q = 2.0) solving Burgers' equation (13) using various kernels.
One can possibly find a GP with fine-tuned kernel, e.g. rational quadratic (rq), to have better (Q-EP with rbf in Table C.2) matching (Q-EP with Matern52 in Table C.3 results. However,
19

---

FD Solution
Soluton (g-1.0)
Error 10
Uncertainty (q=1.0)
3.24
0416
625
825
0000000
-3.9
-1.00
-1.0
-0.40
-615
1.5
Solution (PINN)
Solution (q=2.0)
6.00
1.0
Error 1em2.0)
Muace
1.
Uncertainty (qm2.0)
00032
030
0.
850
68
0.050
018
0h
044
0.4
-0.2
01
-6,23 M41
0.2
022
02
00334
60033
55
03
-1.00
05
-0.30
0.0
00038
Figure C.I: Solving Burgers' equation L
using high-resolution finite difference method (upper left), PINN (lower left), Q-EP (right three in upper row: g = 1.0), and GP (right three in lower row:
q = 2.0) respectively. Blue crosses are learned inducing points.
Table C.3: Comparing accuracy Q-EP (q = 1) against GP (q = 2) solvers with various kernels Burgers' equation (K) terms mean absolute error (MAE), mean squared error (MSE), and relative errors L₁ norm 1). La norm (RLE-2), and Lo norm (RLE-oc) respectively. Result in each cell are averaged over 10 experiments with different random seeds; values after + are standard deviations of these repeated experiments.
Model (g)
kemnel
MAE
MSE
RLE-1
RLE-2
RL-20
1.0
Matern
9.33e-3 = 1.76c-4
1.77e-4 2.16c-4
0,0242
0.0128
0,0266 = 0.0140
0.0324 : 0.0145
2.0 (Gaussiun)
Matern
7.13e-2 + 8.68e-3
1.08e-2 1.28e-2
0.1848 + 0.0908
0.2068 + 0.1118
0.2376 +: 0. 0.1456
1.0
rbf
2.51e-3 1.30e-3
1.37e-5
1.59c-5
0.0065
0.0074 0.0040
0.0150
2.0 (Gaussian)
rbf
2.49e-2 = 3.72c-3
1.06e-3
2.63e-4
0.0646 = 0.0097
0.0726 = 0.0094
0.1144 : 0.0073
1.0
rq
2.57e-3 = 5.52e-4
1.39e-5
3.15e-6
0.0067 : 0.0014
0.0083 = 0.0010
0.0199 : 0.0049
2.0 (Gaussian)
rq
1.97e-3
7.02e-5
0.0283
0.0057
0.0551
principled regularization over function spaces, Q-EP facilitates effective modeling of derivatives and differential equations, thereby alleviating the struggle of GP on meticulous kernel engineering.
C,3 Nonlinear Elliptic Equation
Now we consider nonlinear elliptic equation with Dirichlet boundary condition on = 1]::
-Au(x) + 7 (u(x)) f(x), x e 52,
(21)
u(x) =((x,, 0Ω.
where we choose (u)) = w., g(x) 0, and f such that the following true solution at satisfies equation (2l
6
wt (x) sin(TTx1) sin()Tx2) + sin (allee))snnnninnnn
We follow the variational procedure in Section 3.1 and use the same experiment setup in Section 5 Figure C.2 compares the solutions by a highly-resolved finite difference method (upper left, treated as true solution), PINN (lower left), Q-EP (right three in upper row: = 1.0), and GP (right three in lower row: q 2.0). Compared the true solution, Q-EP with q 1.0 yields an estimate more accurate than PINN and GP. Based on Table C.4 the best solution is obtained by Q-EP with q = 0.5. Since the true solution is highly fluctuating over the domain 12, imposing more challenges on the boundary, where the largest pointwise errors occur. Adding more boundary points might help with Q-EP solvers.
C.4 Inverse Darcy Flow
Finally, we consider more realistic Darcy flow data used by Fourier Neural Operator (FNO) and Physics-Informet Neural Operator (PINO) available IVIDIA PhysicsNeMo This dataset contains thousands of permeability-solution pairs that reflect realistic porous media. Since our method is
20

---

True Solutio
Solution (3=1.))
Error 10-1.0)
Uncertainty (g=1.0)
00024
ax
600335
25
13
00015
00012
0.
6.022
Solution (PINN)
Solution (qm2.0)
Error 1e=2.0)
Uncertainty (qm2.0)
50
00012
10035
0.
52
02
02
06
04
06
04
0.8
Figure C.2: Solving nonlinear elliptic equation (2T) using high-resolution finite difference method
(upper left), PINN (lower left), Q-EP (right three in upper row: 1.0), and GP (right three in lower row: q = 2.0) respectively. Blue crosses are learned inducing points.
Table C.4: Comparing accuracy of various solvers for nonlinear elliptic equation 21 terms of relative error in norm (RLE- 1), Le norm (RLE-2), Loe norm (RLE-3o), and per iteration. Result each cell are averaged over 10 experiments with different random seeds; values after + are standard deviations these repeated experiments.
Model (g)
RLE-1
RLE-2
RLE-x
time/iteration
PINN
1.6428 + 0.7133
1.5979 0.6536
1.5595 1: 0.4462
0,0052 1: 0.0023
B-PINN
0.7026 = 0.2879
0.7868 = 0.2872
1.2395 = 0.3480
0.0108 = 0.0013
0.5
0.0341 + 0.0203
0.0514 0.0322
0.1617 +: 0.0880
0.2530 + 0.0008
0.8
0.3020 = 5.0915
0.3159 0.0749
0.5549 + 0.1092
0.2604 + 0.0016
1.0
0.2835 + 0.0532
0.2958 + 0.0434
0.5122 + 0.0870
0.2614 + 0.0005
1.2
0.3379 = 0.0845
0.3342 0.0689
0.5229 + 0.1094
0.2714 + 0.0004
1.5
0.3953 = 0.1058
0.4429 + 0.1591
0.8938 = 0.5714
0.2645 0.0003
2.0(Gaussian)
0.9876 = 0.2777
1.1028 = 0.3254
2.0255 + 1.1255
0.2562 = 0.0004
not to train surrogate model, take only one pair and impose the data 200 X 200 mesh to obtain 2000 randomly sampled observations. We then train Q-EP solvers on 60 x 60 mesh (with collocation points taken from the grid) and with 512 inducing points. This problem has much larger scale (~ 20 times larger) than all the previous examples. As illustrated in Table C.6 Q-EP with q still achieves remarkable advantage compared with GP (g and PINN. Note, Q-EP with q = 1 has comparable accuracy with & = 0.5 in forward solution. With only one training pair in the NIVIDA FNO-Darcy data, it understandably more challenging obtain the inverse solution, yet for which Q-EP with attains the best accuracy. In both solutions, GP (q 2) is much worse by metries.
Table C.5: Comparing accuracy of various solvers for inverse Darcy flow 14 in terms of relative error in Li (RLE-I), L2 norm (RLE-2), and Los (RLE-co) respectively. Result each cell averaged 10 experiments with different random seeds; values after are standard deviations of these repeated experiments.
Forward Solurion
Inverse Solution
Model (g)
RLE-1
RLB-2
RLE-00
RLE-1
RLE-2
RLE-20
PINN
0.3145
0.0163
0.3331
0.0137
0.4225
0.0250
6.6922
0.5800
6.2856
0.5432
4,1070
0.3189
0.5
0.2397
0.0975
0.2689
0.37
0.3445
0.1490
0.4827
0.2011
0.5217
0.1916
0.8711
0,144
D.8
0.1666
0.0242
0.1784
0.0248
0.2155
0.0321
0.3577
0.0758
0.4154
0.0675
0.810
+ 0.1239
1.0
0,1354
0,0286
0,1432
0.0253
0,1810
0.0193
0.4340
0.0376
0.4942
0.0376
0,8494
0.2132
1.2
0.1523
0.0967
0.1549
0.0898
0.2066
0.1035
0,5468
0.1078
0.5921
0.0890
0.7801
0.1301
1.5
0.3270
0.3257
0.4723
0.5048
1,2103
1.0981315
1.6317
0.8162
1.5914
0.7810
1.7561
0.8064
2.00 Gaussian)
0.6328
0.5954
1.0658
0.9186
3.1391 = 2.2492
2,2096 - 0.9100
2.1675
0.8823
2,8022
1.6839
21

---

The Inverse
Inverse Iqe2.0
Uscemt Frcare Selation
08
8524
0.50
6006
16
66
0600
0666
8632
66
04
5012
000
10
5005
1.0
FEM Solution
8.00
PD2 Solution (.w22.01
Invens Soldiion (ow2
uncemsrry twers Soupion 91-44
60032
6.
48
-0.
0.0
08
6025
60028
6015
43
60312
02
-3
05
3.0
02
04
06
00
1.0
0.000
0.2
04
3.0
04
02
04
05
00 -62
02
06
0.0
0.330
Figure C.3: Solving inverse Darey flow (14) with sparse data using Q-EP (right three in upper row: = 1.0) and GP (right three in lower row: 2.0) respectively. Upper left: true inverse solution a'; lower left: fine-resolution finite element solution u to (14) with a!, Blue crosses are learned inducing points, and red dots indicate locations of observations.
Solution (PINN)
Inverse Solution (PAN)
Febetien La Ery
Febeiinn
hndyet bror
045
-55
-20
-25
-22
Figure C.4: Solving inverse Darcy flow
14 with PINN running for longer (5000) iterations.
Table C.6: Comparing accuracy of various solvers for the inverse problem with FNO-Darcy data in terms of
relative error in I.1 norm (RLE-1),, L2 norm (RLF-22), and Too norm (RIF-oo) respectively. Result each cell are averaged over 10 experiments with different random seeds; values after 1: are standard deviations of these repeated experiments.
Forward PDE Solation
Inverse Solution
Model (g)
RLE-1
RLE-2
RLE
RLE-1
RLE-2
RLE
PINN
0,9019
0.8170
0.8190 0.3859
1.5170 0.1547
0.5
0.3316
0.0677
0,3191 = 0.0402
0.5978 # 0.0535
0.7847 = 0.1394
0.8548
0.1683
1.2072 : 0.2975
1.0
0.3239
0.1297
0.3195
0.0985
0,6041 0.0830
0.7724 : 0.1798
0.8349
0.1841
1.2895
:: 0.2849
1.5
0.2338
0.0197
0,2490
0.0105
0.6149
= 0.1056
0.8466 = 0.1617
0.8744
0.1213
1.2070
= 0.3958
2.0(Gaussian)
0.5341 = 0.5296
0.4904
0.4292
0.6881
0.2356
0.9124 :: 0.1791
0.1502 0.9251
1.1018
+ 0.1219
Nelate
2013
200t
ba un 28410
Eow
02
Figure C.5: Comparing solutions to inverse Darcy flow 14) at refined mesh sizes using Q-EP (left four in upper row: g = 1.0) and GP (left four lower row: 4 = 2.0) respectively, with relative enrors in L1 norm (upper right) and L2 norm (lower right). Blue crosses are learned inducing points. Error bars indicate standard errors based on 10 repeated experiments with different random seeds.
22

---

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

---

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

---

5. Open access to data and code
Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification: Sample codes are included in the supplementary materials and all the codes will be released upon acceptance. Data are simulated and details on how to generate them are included in the section Numerical Experiments and the Technical Appendices.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines () for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" an acceptable answer. Papers cannot be rejected simply for including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Question:
Answer: [Yes]
Justification: All the details of the experiments are disclosed in the section Numerical Experiments and the Technical Appendices.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core the paper to a level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: Please refer to relevant figures for error bars.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
25

---

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

---

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

---

For scraped data from a particular source (e.g., website). the copyright and terms of service of that souree should be provided.
assets are released, the license, copyright information, and terms of use the package should be provided. For popular datasets, has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset.
For existing datasets that are re- re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors encouraged to reach out to the asset's ss creators.
13. New assets
Question: Are assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
The paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/codd/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsoureing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
28

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: The core method development in this research does involve LLMs as any important, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer our LLM policy () for what should should not be described.
29


Reviews:
[{"review_id": "qPxjEyqbMO", "reviewer": "Reviewer_8z5k", "summary": "The authors propose to use q-exponential process to solve and learn PDEs as they are designed to better handle data with abrupt changes and to more accurately model derivative information.", "questions": "1. how do we select the optimal q in real problems? it seems q=1 is nearly always the best and the convergence is faster in theory. do you suggest using q=1 for all tasks\n\n2. it seems that the uncertainty provided by q=1 appears more structurally informative than that of the GP solver q=2. can you elaborate on the reason for this?\n\n3. i feel that it may be a bit misleading to say q-ep is designed to handle abrupt changes in data. does it only refer to q<1? can you comment on that? also, is it possible to adjust the kernel in GP to address the same issue rather than use q-ep? \n\n4. q-ep does not have such trctable conditional and predictive distributions as GP. how do you overcome that?\n\n5. how is q-ep compared with Besov process", "limitations": "no", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n1. The idea of using Q-EP is quite novel\n2. They have great theoretical justification\n3. Uncertainty quantification is provided\n4. They demonstrate their approach in a range of tasks and show superior performance\n\nWeaknesses:\n\nfor highly nonlinear PDEs, the variational approximation to the true posterior distribution might be poor, which could undermine the final solution", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for acknowledging the novelty and strengths of our work. Our **main contribution** is showing that Q-EP (with q=1) is superior to GP (corresponding to q=2) for solving and learning PDEs -- we have demonstrated this point both theoretically and empirically. GP is so popular in (physics-informed) machine learning that it is almost the universal choice for probabilistic PDE solvers. Our argument is that **an easy but better alternative can be offered by Q-EP**.\n\nWeaknesses:\n\nWe adopted sparse variational Bayes for the computational consideration. It works decently well for all the non-linear PDEs included in our paper and it is **an** inference method we used. We admit that extremely non-linear systems, like weather forecasting equations, could pose challenges for our proposed Q-EP based PDE solver. Yet such extreme non-linearity challenges all approximate methods, including GP solvers and neural network based solutions. If the true posterior distribution is the main focus, advanced Monte Carlo methods, such as those involving optimal transportation and normalizing flow, could be used to mitigate the issue. This, however, is not the weakness or limitation of Q-EP itself.\n\n# Questions:\n1. Yes, we advocate Q-EP with $q=1$ as a superior alternative to GP ($q=2$) for solving and learning PDEs. The convergence rate by Theorem 4.2 attains its peak value at $q=1$, which justifies the optimal choice of $q=1$ -- it converges the fastest. This is also verified in our numerical examples and particularly illustrated in the rightmost column of Figure 3, which shows the error reducing as a function of epochs.\n2. We are unsure about whether the uncertainty by Q-EP with $q=1$ is more structurally informative than that of GP as we did not claim so. Similarly as Q-EP, the posterior by GP solver also indicates the uncertainty mainly comes from the boundary in the Eikonal equation, and GP solver also has higher uncertainty in the shock area in Burgers' equation. We will be happy to discuss on specific examples.\n3. We refer to $q<2$ when talking about improved capability of handling inhomogeneous objects compared to GP (which corresponds to $q=2$). Actually, Q-EP (with different qs) has a regularization effect on modeling objects through the parameter $q$ - smaller $q$ imposes stronger regularization and hence reconstructs images sharper (See Figure 1 of Li et al (2023)). \nIt is true that kernel choice impacts the model performance. For example, Matern kernels generally perform better than radius basis function (rbf) in imaging analysis. So it is possible, but not necessary, that GP with Matern kernel yields better results than Q-EP with rbf (Refer to the response to Reviewer MsMP below for an added comparison using rbf kernel). Ideal kernel is chosen such that the resulting process matches the smoothness of modeled data, which is usually unknown. Q-EP, on the other hand, offers a rigorous regularization on function spaces through the parameter $q$ rather than carefully choosing kernels. Note, all the comparisons in the paper were based on the same kernel `Matern` ($\\nu=5/2$) to ensure the fairness in comparing Q-EP with GP.\n4. Q-EP DOES have tractable conditional and predictive distributions as GP. This is exactly Theorem 3.5 of Li et al (2023) and quoted as Theorem 2.2 in our paper. This is one of the theoretical and numerical advantages of Q-EP, when compared with Besov.\n5. This is a very good question! It has been answered in Theorem 3.4 of Li et al (2023). Quoted as Theorem 2.1, Q-EP shares the same series representation as Besov process that is originally defined in terms of series. As also mentioned in the introduction, Q-EP can be viewed as a probabilistic definition of Besov and it has tractable posterior prediction and explicit control of correlation strength (through covariance kernel). Numerically, Q-EP outperforms Besov in low dimensions and behaves similarly in higher dimensions in image reconstruction problems (Li et al, 2023). However, directly applying Besov to solving PDE is impeded by theoretical difficulties, e.g. justifying that the differential of Besov is still Besov, how to obtain conditional and prediction, etc."}, {"review_id": "sPurnEaz8i", "reviewer": "Reviewer_rFWA", "summary": "The paper addresses the limitations of neural network and Gaussian Process (GP)-based methods in solving PDEs, particularly for systems with sharp transitions. It introduces the Q-exponential process (Q-EP), a generalisation of GPs that offers more expressive and flexible function estimators in settings involving non-smooth dynamics and derivative information. Using sparse variational inference, the proposed method achieves scalable and accurate PDE solutions with principled uncertainty quantification. Toy experiments on the Eikonal equation, Burgers’ equation, and inverse Darcy flow show that Q-EP outperforms GPs and neural networks in both accuracy and uncertainty estimation.", "questions": "Please refer to the Weaknesses section for Questions.", "limitations": "yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "I quite like the idea of using the Q-exponential process, which can be seen as a generalised version of GPs for solving PDEs, offering a flexible choice of the exponent $q$. It also makes sense to employ sparse variational inference with inducing points to reduce algorithmic complexity when solving the functions. Overall, the proposed method is well-justified and theoretically sound. However, I have some comments and potential limitations below that I hope the authors can address.\n- However, a natural limitation arises in selecting the optimal value of $q$ for different settings. I’m curious about the performance when $q>2$, and is $q=1$ always optimal, as suggested by the experiments? Please clarify this point.\n- Another limitation I notice is the absence of real-world datasets in the experiments. I suggest the authors incorporate experiments on real datasets -such as physics-informed datasets, time series, or, ideally, image datasets. I would highly appreciate this addition, and it would positively influence my evaluation further. \n- Could the authors mention any potential applications of the Q-exponential process beyond solving PDEs and inverse problems? For example, might it be applicable in modelling processes within diffusion models? I’m curious whether such an extension is feasible.\n\nOverall, I lean toward accepting this work, but only weakly, due to two main limitations: the unclear guidance on selecting the exponent $𝑞$, and the lack of experiments on real-world datasets. I will adjust my score based on the authors’ responses to these concerns.\n\n**[After Rebuttal]**: My concerns have been addressed, and I am happy to raise my confidence score in advocating for the acceptance of this paper (from 2 to 4).", "was_revised": true, "final_justification": "All of my concerns have been addressed by the authors. I am keeping my initial decision for this paper (borderline accept) but raising my confidence score from 2 to 4 for this decision.", "rebuttal": "We are grateful to the reviewer for the appreciation of our work! In fact, we advocate Q-EP with **$\\boxed{q=1}$** over GP (corresponding to $q=2$) and other qs for solving and learning PDEs. The reason is that Q-EP converges the fastest at $q=1$, which is the main take-home message of Theorem 4.2. As explained in Remark 3, the posterior convergence attains the maximal rate at $q=1$. This is also verified in our numerical examples and particularly illustrated in the rightmost column of Figure 3, where Q-EP with $q=1$ has the steepest reduction in error in the process of training.\n\n* As mentioned above, we have the optimal choice of $q=1$ from the posterior contraction perspective. Such choice is justified by Theorem 4.2 and the numerical illustrations. We have also added $q=2.5$ per your request, and it is not surprising that its results are worse as shown in the following tables.\n\n| Model (Eikonal) | MAE | MSE | RLE-1 | RLE-2 | RL-$\\infty$ |\n| -----------------  | ----- | ------ | ------- | ------- | ------------- |\n| **1.0** | **1.68e-3** $\\pm$ 5.41e-4 | **4.42e-6** $\\pm$ 2.26e-6 | **0.0106** $\\pm$ 0.0034 | **0.0110** $\\pm$ 0.0030 | **0.0321** $\\pm$ 0.0071 |\n| 2.0 | 9.64e-3 $\\pm$ 1.33e-3 | 1.30e-4 $\\pm$ 3.61e-5 | 0.0610 $\\pm$ 0.0084 | 0.0612 $\\pm$ 0.0087 | 0.1009 $\\pm$ 0.0124 |\n| 2.5 | 7.71e-2 $\\pm$ 3.18e-2 | 8.36e-3 $\\pm$ 5.29e-3 | 0.4881 $\\pm$ 0.2011 | 0.4668 $\\pm$ 0.1743 | 0.4442 $\\pm$ 0.1390 |\n\n| Model (Burgers) | MAE | MSE | RLE-1 | RLE-2 | RL-$\\infty$ |\n| ------------------ | ----- | ------ | ------- | ------- | ------------- |\n| **1.0** | **9.33e-3** $\\pm$ 1.76e-4 | **1.77e-4** $\\pm$ 2.16e-4 | **0.0242** $\\pm$ 0.0128 | **0.0266** $\\pm$ 0.0140 | **0.0324** $\\pm$ 0.0145 |\n| 2.0 | 7.13e-2 $\\pm$ 8.68e-3 | 1.08e-2 $\\pm$ 1.28e-2 | 0.1848 $\\pm$ 0.0908 | 0.2068 $\\pm$ 0.1118 | 0.2376 $\\pm$ 0.1456 |\n| 2.5 | 2.26e-1 $\\pm$ 3.20e-1 | 2.49e-1 $\\pm$ 6.04e-1 | 0.5879 $\\pm$ 0.8296 | 0.6669 $\\pm$ 0.9576 | 0.7558 $\\pm$ 1.0521 |\n\n* We have now added a more realistic Darcy flow data (publicly available on NVIDIA website, link suppressed due to NIPS policy) used by Fourier Neural Operator (FNO, Google \"nvidia FNO darcy\") and Physics-Informed Neural Operator (PINO) available in **NVIDIA PhysicsNeMo**.\nThis dataset contains thousands of permeability-solution pairs that reflect realistic porous media. Since our method is not to train a surrogate model (see the response to Reviewer MsMP below for elaboration on the difference), we take **only ONE** pair and impose the data on 200x200 mesh to obtain 2000 randomly sampled observations. We then train Q-EP solvers on 60x60 mesh (with collocation points taken from the grid) and with 512 inducing points. This problem has much larger scale (~20 times larger) than all the previous examples. As illustrated in the following table, Q-EP with $q=1$ still achieves remarkable advantage compared with GP ($q=2$) and PINN. Note, Q-EP with $q=1$ has comparable accuracy with $q=0.5$ in forward solution. With only one training pair in the NIVIDA FNO-Darcy data, it is understandably more challenging to obtain the inverse solution, yet for which Q-EP with $q=1$ attains the best accuracy. In both solutions, GP ($q=2.0$) is much worse by most metrics.\n\n| Model ($q$) | RLE-1 (fwd) | RLE-2 (fwd) | RLE-$\\infty$ (fwd) | RLE-1 (inv) | RLE-2 (inv) | RLE-$\\infty$ (inv) |\n| -------------- | ------------- | -------------- | ---------------------| ------------- | ------------ | -------------------- |\n| PINN | 0.9168 $\\pm$ 0.9019 | 0.8170 $\\pm$ 0.7272 | 0.8190 $\\pm$ 0.3859 | 1.6550 $\\pm$ 0.1876 | 1.5170 $\\pm$ 0.1547 | 1.4303 $\\pm$ 0.1190 |\n| 0.5 | 0.3316 $\\pm$ 0.0677 | 0.3191 $\\pm$ 0.0402 | **0.5978** $\\pm$ 0.0535 | 0.7847 $\\pm$ 0.1394 | 0.8548 $\\pm$ 0.1683 | 1.2072 $\\pm$ 0.2975 |\n| 1.0 | 0.3239 $\\pm$ 0.1297 | 0.3195 $\\pm$ 0.0985 | 0.6041 $\\pm$ 0.0830 | **0.7724** $\\pm$ 0.1798 | **0.8349** $\\pm$ 0.1841 | 1.2895 $\\pm$ 0.2849 |\n| 1.5 | **0.2338** $\\pm$ 0.0197 | **0.2490** $\\pm$ 0.0105 | 0.6149 $\\pm$ 0.1056 | 0.8466 $\\pm$ 0.1617 | 0.8744 $\\pm$ 0.1213 | 1.2070 $\\pm$ 0.3958 |\n| 2.0 | 0.5341 $\\pm$ 0.5296 | 0.4904 $\\pm$ 0.4292 | 0.6881 $\\pm$ 0.2356 | 0.9124 $\\pm$ 0.1791 | 0.9251 $\\pm$ 0.1502 | **1.1018** $\\pm$ 0.1219 |\n\n* Since proposed by Li et-al (2023), Q-EP has been applied to:\n  - Regularization of latent representation (ICLR 2025)\n  - Deep probabilistic models (AABI 2025)\n\nWe are glad that the reviewer's suggestion resonates with our ongoing efforts -- yes, we have been working on the diffusion model with Q-exponential noise and achieved favorable preliminary results.\n\nWe hope our added results could convince the reviewer that Q-EP is preferable to GP in solving PDEs and promising as a new tool in physics-informed machine-learning."}, {"review_id": "DuNzdO7ZIN", "reviewer": "Reviewer_MsMP", "summary": "This paper extends recent approaches in solving and learning PDEs based on Gaussian processes to a new class of processes called Q-exponential process. The proposed contribution is the more flexible properties for this new process, especially for PDEs that have stiff solutions, which may be impacted by smoothness assumptions coming from standard GPs.", "questions": "1) Theory: In Section 4, there are theoretical results justifying this approach. However, while they look impressive, the message of this section is unclear in the main text. How do these theorems relate to the results on Gaussian processes, if there are any corresponding results? In terms of handling solutions that are less smooth, is there any noticeable improvement over GPs in theory?\n\n2) UQ: The main claim against a PINN based approach is that there is UQ available from GP and q-exponential based approaches. Can authors elaborate why such tightly concentrated uncertainty would be useful? Is there any realistic example where the UQ is meaningful, e.g,. depending on the observation locations (e.g. in a sparser area)? Can authors construct such an example?\n\n3) I'm not entirely convinced that by tuning and changing GP kernels, similar results you have obtained with the q-exponential process cannot be obtained. You have fixed the kernel. Can you elaborate and provide additional results to argue whether the improvements you observe remain with different (suitable) kernels?\n\n4) The comparisons are really thin. Many recent approaches like DeepONets are not compared to - there are similar probabilistic approaches, see e.g. *Vadeboncoeur, Arnaud, et al. \"Fully probabilistic deep models for forward and inverse problems in parametric PDEs.\" Journal of Computational Physics 491 (2023): 112369* and references therein.\n\nI think the experimental section needs much more work before this paper is ready for NeurIPS or a similar venue.", "limitations": "Yes, very briefly. It could be expanded.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths: This paper is based on a novel idea of using q-exponential processes in the setting of modelling PDEs. This can lead to interesting applications, especially in the settings author mentioned, e.g., PDEs with erratic solutions.\n\nWeaknesses: I list the weaknesses below - please see questions part for further inquiries:\n\n(1) The paper is somewhat an incremental improvement over the previous approaches. The methodology is rather standard, replacing GP with q-exponential process, but deriving a standard variational inference procedure. Therefore, despite this is an interesting extension, I find the whole framework's novelty a bit limited.\n\n(2) Experimental results are insufficient and not very conclusive in my opinion - to justify the whole framework introduced in this paper.\n\n(3) There are some theoretical results - but zero intuition about how to interpret them. Please see my questions.", "was_revised": true, "final_justification": "The authors addressed my concerns -- by including numerical evidence, which tilted the paper from Borderline reject to borderline accept.\n\nIt is not still a \"full accept\" for me. I find this general direction not very innovative - since the inference techniques are very costly, and scalability of these methods (GPs) in general are problematic, thus I doubt this will become a usual way to solve these PDEs.", "rebuttal": "Thanks to the reviewer for acknowledging the novelty of our work as an interesting extension to GP in solving PDEs. Such extension by Q-EP has theoretical and numerical advantages in modeling inhomogeneous objects and derivative information. The novelty lies in the first application of Q-EP to solving and learning PDEs and successful demonstration of its advantages over GP.\n\nWeaknesses:\n\n(1) We admit that we did not invent new inference methods (and it is not our focus). The adoption of sparse variational inference is not standard or straightforward: the unique format of Q-EP density as in Equation (1) prevents many tractable results in ELBO of GP; therefore it requires significant work to overcome these obstacles in deriving computable ELBO for Q-EP. Additionally, as a Bayesian generalization to Chen et~al (2021), which only provides MAP solution, we creatively propose *distribution propagation* (as in Equation (5)) to properly propagate uncertainty (posterior standard deviation) for appropriate UQ, which is not seen in the existing probabilistic PDE solvers. What's more, the posterior contraction theorems applied to solving and learning PDEs represent a tremendous theoretical development. Therefore, we believe our work makes a significant contribution to the physics-informed machine learning and should not be seen merely as a plug-in extension of existing methods.\n\n(2) We have added new experimental results per your request (see below).\n\n(3) We explained the implication of theorems, particularly theorem 4.2, in remark 3. We will elaborate the points in our revision (and below).\n\n# Questions:\n\n1. Long story short. The take-home message of Theorem 4.2 is that the posterior contraction rate of our Bayesian Q-EP solvers attains maximal value at $q=1$ so Q-EP with $q=1$ converges the fastest and should be preferred to GP which is a special case of Q-EP with $q=2$. More specifically, as long as the true PDE solution is at least $L_1$ regular, Q_EP solver for $q=1$ should be adopted for the fastest convergence. This is illustrated by our numerical examples and particularly in the rightmost column of Figure 3 where Q-EP ($q=1$) reduces the error the fastest in the training process.\n\n2. In forward (solving) problems, uncertainty reflects the quality of variational inference—specifically, low uncertainty indicates a well-configured set of inducing points. In inverse (learning) problems, uncertainty reflects the informativeness of the observations—low uncertainty suggests the presence of rich and informative data. One can refer to Figure C.3 and Figure C.4 for a comparison on the uncertainty of inverse solution with different amount of observations. The high uncertainty region usually signals interesting phenomena or data insufficiency, such as the middle region of shocks at $x=0$ in Burgers' equation illustrated in the rightmost columns of Figure C.1. This could be informative and useful for downstream applications such as design of experiment (DoE).\n\n3. A similar question was asked by Reviewer 8z5k as above. Yes, you can possibly find a GP with fine-tuned kernel to beat Q-EP with a poor kernel. The meaning of using Q-EP is to rigorously regularize on function spaces -- this is especially important when it comes to solving PDEs where large derivatives can easily blow up and hurdle the training. Q-EP, compared to GP, imposes stronger regularization with $q<2$ on the derivatives when solving PDE by modeling the function and its derivatives as Q-EPs and hence enjoys a smoother training and faster convergence. To ease your doubts, we have added another comparison with a radius basis function (rbf) kernel, which is kept the same for different models in the following table summarizing the results of solving Burgers' equation. Note, compared to Table 1 with Matern($\\nu=5/2$) kernel, the overall results with rbf are improved (maybe due to the solution smoothness being better matched). **Nevertheless, the Q-EP ($q=1$) with Matern($\\nu=5/2$) kernel yields higher accuracy than GP ($q=2$) with rbf kernel!** Note, the hyper-parameters of both kernels were automatically optimized by `GPyTorch` in all experiments.\n\n| Model (Burgers-**rbf**) | MAE | MSE | RLE-1 | RLE-2 | RL-$\\infty$ |\n| ------------------ | ----- | ------ | ------- | ------- | ------------- |\n| **1.0** | **2.51e-3** $\\pm$ 1.30e-3 | **1.37e-5** $\\pm$ 1.59e-5 | **0.0065** $\\pm$ 0.0033 | **0.0074** $\\pm$ 0.0040 | **0.0150** $\\pm$ 0.0088 |\n| 1.5 | 3.13e-2 $\\pm$ 3.55e-3 | 1.62e-3 $\\pm$ 3.42e-4 | 0.0813 $\\pm$ 0.0092 | 0.0900 $\\pm$ 0.0094 | 0.1268 $\\pm$ 0.0115 |\n| 2.0 | 2.49e-2 $\\pm$ 3.72e-3 | 1.06e-3 $\\pm$ 2.63e-4 | 0.0646 $\\pm$ 0.0097 | 0.0726 $\\pm$ 0.0094 | 0.1144 $\\pm$ 0.0073 |\n| 2.5 | 3.27e-2 $\\pm$ 7.98e-3 | 1.74e-3 $\\pm$ 9.48e-4 | 0.0850 $\\pm$ 0.0207 | 0.0917 $\\pm$ 0.0243 | 0.1220 $\\pm$ 0.0089 |\n\n4. Thanks for the references. First, we want to contend that the numerical results are not thin. We have compared various Q-EP with PINN, B-PINN, GP, and Chen et al (2021) using 4 examples with different focuses including accuracy (tables), convergence (Figure 3), UQ (Figures 2, C.1 and more), influence of data amount (Figures C.3 vs C.4) and resolution (Figure C.6) (refer to the appendix). These examples are standard benchmarks in similar works in the area, including (but not limited to) Chen et~al (2021), Meng and Yang (2023), and more recently Hamelijnck et al (2024). Second, we have also added another more realistic inverse problem from NVIDIA FNO-Darcy dataset (See above response to Reviewer rFWA). Last but not the least, our method is **NOT comparable** to neural operator methods including DeepONets, FNO, PINO and many more:\n\n    * Neural operators require thousands of pairs of inputs (functions of interest such as coefficients, boundaries, etc) and outputs (PDE solutions) to train the neural networks. This usually involves a numerical solver to obtain solutions to the PDE with different inputs. On the contrary, our method does **not require** to prepare solutions as training samples. Instead, it models PDE solution as Q-EP and solves the PDE using variational inference --  the posterior mean is the estimate to the PDE solution and the posterior standard deviation quantifies the uncertainty of the solver.\n    * Neural operators are trained to provide **surrogate** models for quick forward analysis or inverse UQ. While Q-EP probabilistic solver can also get inverse solution (with forward solution as a byproduct, see Section 3.3), it does NOT rely on any surrogate model.\n\nTherefore, it would be an apple-to-orange comparison if comparing Q-EP solver with DeepONet. In checking with our motivating work, Chen et al (2021), and the follow-ups, Meng and Yang (2023) and Hamelijnck et al (2024), none of them compared with neural operator methods. Therefore, we will not include any neural operator method such as DeepONet in our comparison though we appreciate the reviewer's suggestion."}, {"review_id": "d7ffldYqM7", "reviewer": "Reviewer_C3Mw", "summary": "The paper proposes q exponential distributions as an alternative to GPs for solving PDEs with uncertainty quantification. This can be seen as a generalization (q=2 is Gaussian). Specifically for derivative information (which is very relevant for PDEs), the method is claimed as more accurate. The method is well presented within existing literature, and modelling PDEs with uncertainty is a relevant problem for the community. The method is validated on some simple toy problems.", "questions": "Q1: The q-exponential distribution seems similar to the generalized Gaussian distribution $\\propto \\exp(-\\frac{1}{2}(|x-\\mu|/\\alpha)^\\beta)$, can you please comment on the differences in the context of your paper?\n\nQ2: What is matern52? (line 136)\n\nQ3: Apart from replacing the Gaussian with the q-exponential, in what ways does your work differ from Hamelijnck, Solin et al.?\n\nMinor questions and remarks:\n- line 56: preferable\n- line 61: An emerging challenge\n- lines 85-88. Very long sentence that is not totally understandable\n- Fig. 1: labels are too small. Add labels indicating the rows (q=1, q=2, ..). Same for Fig. 2\n- (very minor remark) take a look at your references: **B**esov, **G**aussian, \"amp;\", ..\n- line 508: satisfies", "limitations": "In the checklist the authors refer to the scientific contributions at the end of the introduction. I could not find a discussion on societal impacts, however, this should be trivial for this kind of theoretical work.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n- Extensive exposition on the methods.\n- Learning stochastic PDEs from data is a relevant problem, and not trivial.\n- The authors focus on more accurate gradients, and show improvements, both theoretically and empirically.\n\nWeaknesses:\n- The method is validated on a limited number of toy problems.\n- Writing quality is a bit on the low side. This makes the paper a bit harder to read. Some concrete (friendly) suggestions:\nline 99: The negative log density of q-ED (eq. 1) yields ..\nline 102: Li et al. [24] prove ..\nline 126: Li et al. [24, Theorem 3.5] show that ..\nline 131: remove the sentence \"We often receive ..\"", "was_revised": true, "final_justification": "My main concern is the incremental improvement over Hamelijnck et al (2024), Physics-Informed Variational State-Space Gaussian Processes (PHYSS-GP): replacing the Gaussian with Q-EP.", "rebuttal": "We thank the reviewer for the precise summary of the strengths of our paper. Regarding the numerical experiments, these nonlinear PDEs serve as standard benchmarks in similar/related works such as Li et al (FNO, 2021) and Chen et al (2021). We have also added a large scale inverse problem with NVIDIA FNO-Darcy dataset (See the above response to Reviewer rFWA) to strengthen the numerical evidences. The adoption of Q-EP to solve and learn PDEs and the successful demonstration of its superiority over GP consist of the major contributions to the field of physics-informed machine learning. We will also correct the typos and improve the wording of certain sentences according to the reviewer's suggestions.\n\n# Questions:\n\nQ1. You are absolutely right. Univariate Q-EP coincides with the exponential power distribution (a.k.a. generalized Gaussian distribution), as shown in Equation (3) of Li et al (2023). The challenge of developing a valid stochastic process out of it lies in how to generalize the univariate random variable to a multivariate random vector based on which a process can be defined (similar to using MVN to define GP). The Kolmogorov extension theorem determines the requirements which are not trivially met by arbitrary multivariate generalization. For example, the multivariate extension by Gomez et al (1998) fails to satisfy the consistency condition. Q-EP with density in Equation (1) satisfies Kolmogorov's conditions. One can refer to Section 3 of Li et al (2023) for more details.\n\nQ2. `matern52` refers to the Matern kernel with degree $\\nu=5/2$. That jargon came from `GPyTorch` and slipped in. We will clarify it in the revision.\n\nQ3. This is a very good question! The work by Hamelijnck et al (2024), Physics-Informed Variational State-Space Gaussian Processes (PHYSS-GP) is also a probabilistic PDE solver based on GP and variational inference. We cited it in the *Connection to the literature* paragraph in the introduction and pointed out it is a closely-related work. This work differs from ours in the following aspects:\n\n  * PHYSS-GP builds on GP; while our work uses more general Q-EP that rigorously regularizes derivatives and outperforms GP in solving PDEs - that's the key contribution of our work.\n  * PHYSS-GP is based on spatiotemporal state-space model which is tailored to solving **time-dependent** PDEs; our framework is for general PDEs and can also be adapted using spatiotemporal Q-EP priors. However, it is unclear how PHYSS-GP solves **time-independent** PDEs and what its benefit is.\n  * PHYSS-GP linearizes a non-linear transformation regarding collocation points in its EKS prior; our Q-EP linearizes the non-linear PDE mapping directly to properly propagate variational distribution. These linearizations are different.\n\nMinor questions and remarks:\n\nThanks for your careful reading and we will fix them in the revision.\n\n# Limitations:\n\nWe discussed the limitation of variational inference for UQ in the concluding Section 6. The societal impacts will be realized through advancements in physics-informed machine learning."}]

Author rebuttal:

Meta-review: The reviewed work revisits the use of Gaussian processes for the solution of PDEs by replacing the Gaussian with the q-Gaussian. The authors observe that choosing $q=1$ outperforms the Gaussian (q=2) case on multiple different PDEs. The reviewers appreciate the idea of using the q-Gaussian process, although they are concerned about the additional algorithmic complication due to the use of variational approximations. They are also concerned that the work is somewhat incremental compared to prior approaches. Overall this is a borderline case, but I tend to think that the improvements due to $q=1$ are surprising enough to make this work relevant to the NeurIPS community.
Decision: Accept (poster)