Solving High Frequency and Multi-Scale PDEs with Gaussian Processes
Authors: Shikai Fang
Madison Cooley
Da Long
Shibo Li
Mike Kirby
Shandian Zhe
Venue: ICLR 2024 poster (2024)
Keywords: ML PDE solver
gaussian process
PINN
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)

Abstract: Machine learning based solvers have garnered much attention in physical simulation and scientific computing, with a prominent example, physics-informed neural networks (PINNs). However, PINNs often struggle to solve high-frequency and multi-scale PDEs, which can be due to spectral bias during neural network training. To address this problem, we resort to the Gaussian process (GP) framework. To flexibly capture the dominant frequencies, we model the power spectrum of the PDE solution with a student $t$ mixture or Gaussian mixture. We apply the inverse Fourier transform to obtain the covariance function (by  Wiener-Khinchin theorem). The covariance derived from the Gaussian mixture spectrum corresponds to the known spectral mixture kernel. Next,  we estimate the mixture weights in the log domain, which we show is equivalent to placing a Jeffreys prior. It automatically induces sparsity, prunes excessive frequencies, and adjusts the remaining toward the ground truth. Third, to enable efficient and scalable computation on massive collocation points, which are critical to capture high frequencies, we place the collocation points on a grid, and multiply our covariance function at each input dimension. We use the GP conditional mean to predict the solution and its derivatives so as to fit the boundary condition and the equation itself. 
As a result, we can derive a Kronecker product structure in the covariance matrix. We use Kronecker product properties and multilinear algebra to promote computational efficiency and scalability, without low-rank approximations. We show the advantage of our method in systematic experiments. The code is released at {https://github.com/xuangu-fang/Gaussian-Process-Slover-for-High-Freq-PDE}.

Published as a conference paper at ICLR 2024
SOLVING HIGH FREQUENCY AND MULTI-SCALE PDES WITH GAUSSIAN PROCESSES
Shikai Fang,, Madison Cooley,, Da Long", Shibo Li, Robert M. Kirby, Shandian Zhei University of Utah, Salt Lake City, UT 84112, USA (shi kai, wwmolcolyy..oo eyy d1 932, sh bo, k kirby, zhe êcs. .uwh.eduuuu
ABSTRACT
Machine learning based solvers have garnered much attention in physical simula- tion and scientific computing, with a prominent example, physics-informed neural networks (PINNs). However, PINNs often struggle to solve multi-scale PDEs. which can due to the spectral bias during neural network train- ing. To address this problem, we resort to the Gaussian process (GP) framework. To flexibly capture the dominant frequencies, we model the power spectrum of the PDE solution with student mixture or Gaussian mixture. We apply inverse Fourier transform to obtain the covariance function (by Wiener-Khinchin theorem). The covariance derived from the Gaussian mixture spectrum corresponds to the known spectral mixture kernel. Next, we estimate the mixture weights in the log domain, which we show is equivalent to placing a Jeffreys prior. It automatically induces sparsity, prunes excessive frequencies, and adjusts the remaining toward the ground truth. Third, to enable efficient and scalable computation collocation points, which are critical to capture high frequencies, we place the collocation points on grid, and multiply our covariance function at each input di- mension. We use the GP conditional mean to predict the solution and its derivatives as to fit the boundary condition and the equation itself. As result, we can derive a Kronecker product structure in the covariance matrix. We use Kronecker product properties and multilinear algebra to promote computational efficiency and scala- bility, without low-rank approximations. We show the advantage of our method in systematic experiments. The code is released at xuangu- fang/ Gaus wawuussian-PProesssee Slover- for High Preq
and
massive
com/
1 Introduction
Scientific and engineering problems often demand we solve set of partial differential equations (PDEs). Recently, machine learning (ML) solvers have attracted much attention. Compared to traditional numerical methods, ML solvers do not require complex mesh designs and sophisticated numerical tricks, are simple to implement, and can solve inverse problems efficiently and conveniently. The most popular ML solver is the thysics-informed neural network (PINN) (Raissi al. 2019). Consider a PDE of the following general form,
F(u](x)= f(x) (x E 2), =(x)=9(x) (x E a?), =
where F is the differential operator, Ω is the domain, and 0Ω is the boundary of the domain. To solve the PDE, the PINN uses a deep neural network (NN) (x()) to model the solution u. samples No collocation points [
No i from Z and N, points (x) No j=1 from an, and minimizes loss,
0' argming (,(0) L,),
(2)
where L,) 1 No
N, E (x())) g(x))) is the boundary term fit the boundary condition, No M (x f(x!)) is the residual term to fit the equation.
and T,,( N.
Despite many success stories, the PINN often struggles to solve PDEs with high-frequenc, and multi-scal components in the solutions. This is consistent with the "spectrum bias" observed in
"Equal contribution 1 "orresponning author

---

Published as a conference paper at ICLR 2024
NN training Rahaman et al. (019). That is, NNs typically can learn the low- now-freuueccy information efficiently but grasping the high-frequency knowledge is much harder. To alleviate this problem, the recent work Wang et al. (202Ib) proposes to construct set of random Fourier features from zero-mean Gaussian distributions. The random features are then fed into the PINN layers for training (see (2). While effective, the performance of this method is unstable, and is highly sensitive the number and scales of the Gaussian variances, which are difficult to choose beforehand.
In this paper, we resort to an alternative arising ML solver framework, Gaussian processes (GP) (Chen et al. 2021 Long et al. 2022a). We propose GP-HM, a GP solver for High frequency and Multi-scale PDEs. By leveraging the Wiener- Khinchnn theorem, we can directly model the solution in the frequency domain and estimate the target frequencies from the covariance function. We then develop an efficient learning algorithm to scale up to massive collocation points, which are critical to capture high frequencies. The major contributions of our work are as follows.
Model. To flexibly capture the dominant frequencies, we use a mixture of student t or Gaussian distributions to model the power spectrum of the solution. According to the Wiener-Khinchin theorem, we can derive the GP covariance function via inverse Fourier transform, which contains the component weights and frequency parameters. We show that by estimating the weights in the log domain, is equivalent to assigning each weight a Jeffreys prior, which induces strong sparsity, automatically removes excessive frequency components, and drives the remaining toward the ground-truth. This way our GP can effectively extract the solution frequencies. Our covariance function derived from the Gaussian mixture power spectrum corresponds to the known spectral mixture kernel. We therefore are the first realize its rationale and benefit for solving high-frequency and multi-scale PDEs.
Algorithm. To enable efficient computation, we place all the collocation points and the boundary (and/or initial) points on a grid, and model the solution values at the grid with the GP finite projection. To obtain the derivative values in the equation, we compute the GP conditional mean via kernel differentiation. Next, we multiply our covariance function at each input dimension to obtain a product covariance. We then derive a Kronecker product form for the covariance and cross-covariance matrices. We use the properties of the Kronecker product and multilinear algebra to restrict the covariance matrix calculation to substantially reduce the cost and handle massive
each input dimension. In this way, collocation points, without any low rank approximations.
Result. We evaluated GP-HM with several benchmark PDEs that have high-frequency and multi-scale solutions. We compared with the standard PINN and several state-of-the-art variants. We compared with spectral methods (Boyd) 2001) that linearly combine a set of trigonometric bases to estimate the solution. We also compared with several other traditional numerical solvers. In all the cases, GP-HM consistently achieves relative L2 errors at ~ 10-3 or ~ 10-4 or even smaller. By contrast, the competing ML based approaches often failed and gave much larger errors. The visualization of the element-wise prediction error shows that GP-HM also recovers the local solution values much better. We examined the learned frequency parameters, which match the ground-truth. Our ablation study as in Section C of Appendix also shows enough collocation points is critical to the success, implying the importance of our efficient learning method.
2 Gaussian Process
Gaussian processes (GPs) provide an expressive framework for function estimation. Suppose given a training dataset D = {(xn, yn)]1 S n S N), we aim to estimate target function f RRdd - IR. We can assign GP prior,
((-) ~ (P(()))), cov(())),
where m(.) is the mean function and cov(::) :) is the covariance function. In practice, one often sets m(-) 0, and adopts kernel function as the covariance function, namely cov (f(x), fxx)) x(x,x') A nice property of the GP prior is that if f is sampled from a GP, then any derivative (if existent) of f is also a GP, and the covariance between the derivative and the function f is the derivative of the kernel function w.r.t the same input variable(s). For example,
Cov f(x), fxx),((x)) 02 xxx,x'),
(3)
2

---

Published as a conference paper at ICLR 2024
where
Id)" and x' (2) ..... 23) Under the GP prior, the function values at any finite input collection, f = [f(x1) ) f(xx)) ).. follow a multi-variate Gaussian distribution, p(f) (f K) where [K], = cf fx fxx)) xxx,,,,). This is called a GP projection. Suppose given f, we want to compute the distribution of the function value at any input x, namely p(f()))f) Since and f(x) also follow multi-variate Gaussian distribution, we obtain conditional Gaussian, (x)(x) N (f(x)) o² (x))), where the conditional mean
u(x)= = (f),
(4)
and of (x) cov( (x). f(x)) [k(x,x1),... xi)..... kxx,, xx) and X XN]]
f) f(x)), cov(f(x), k(x, X)
3 Gaussian Process PDE Solvers
Covariance Design. When the PDE solution includes high frequencies or multi-scale information, one naturally wants to estimate these target frequencies outright in the frequency domain. To this end, consider the solution's 's power spectrum, S(s) [0) where ū(s) is the Fourier transform of u, and s denotes the frequency. The power spectrum characterizes the strength of every possible frequency within the solution. To flexibly capture the dominant high and/or multi-scale frequencies, we use a mixture of student t distributions to model the power spectrum,
S(s) = E g=l St(s; Hgg Pgy 2).
(5)
where Wi > O is the weight of component g, St stands for student / distribution, Pg is the mean, P is the inverse variance, and v is the degree of freedom. Note that W, does not need to be normalized (their summation is not necessary to be one). Each student & distribution characterizes one principle frequency Pag and also robustly models the (potentially many) minor frequencies with fat tailed density (Bishop 2007). An alternative choice is a mixture of Gaussian, S(s) Lg (s:: Myg P²) But the Gaussian distribution has thin tails, hence sensitive to long-tail outliers and can be less robust (in capturing minor frequencies).
Next, we convert the spectrum model into a covariance function to enable our GP solver to flexibly estimate the target frequencies. According to the Wiener-Khinchin theorem (Wiener) 1930 Khint chine 1934). for a wide-sense stationary random process, under mild conditions, its power spectrum and the auto-correlation form a Fourier pair. We model the solution u as drawn from stationary GP, and the auto-correlation is the covariance function, denoted by k(x,, = k(x - x'). We then have
S(s) - k(())e 12 dz, k(=) = f S(s) (2))e
(6)
where = x - and i indicates complex numbers. Therefore, we can obtain the covariance function by applying the inverse Fourier transform over Ss)). However, the straightforward mixture in (5) will lead to complex-valued covariance function. To obtain real-valued covariance, inside each component we add another student distribution with mean -uy so as to cancel out the imaginary part after integration. In addition, to make the derivation convenient, scale the inverse variance and degree of freedom by a constant. We use the following power spectrum model,
S(s) Σ.e
9-1 wy (St(s; Mq: 4T² Pg: 2v) + St(s; Myy 4T Pgr 2v))
Applying inverse Fourier transform in (6). we can derive the following covariance function,
KSMM(I,
9=1 Tu.p, (x,, x²) cos (2TT - z²));
(
where Ne.,, (x, 2² " r(v)
(V2ylk-x)) Py. (v 1=-x11 is the Matérn kernel with degree of Po
freedom 1 and length scale Pq. and K, is the modified Bessel function of the second kind. The details of the derivation is left in Appendix. We now can see that the frequency information Hlg and
'To be pell-posed, the power spectrum for random process is defined in a slightly different way (taking the limit of windowed signal), but it reflects the same insight; see (1 Lathi 1998 Grimmett and Stirzaker 2020) for details.
3

---

Published as a conference paper at ICLR 2024
component weights Wi are embedded into the covariance function. By learning GP model, expect to capture the true frequencies of the solution. One can also construct a symmetric Gaussian mixture in the same way, and via inverse Fourier transform obtain
(1,
Wi exp (-p - z²)9) cos(22((2 z²) Ply)
(9)
q=1
This is known as the spectral mixture kernel (Wilson and Adams 2013) which was originally proposed to construct an expressive stationary kernel according its Fourier decomposition, because in principle the Gaussian mixture can well approximate any distribution, as long as using enough many components. Wilson and Adams (2013) showed that the spectral mixture kernel can well recover many popular kernels, such as rational quadratic and periodic kerne.. In this paper, we take a different motivation and viewpoint. We argue that a similar design can be very effective extracting dominant frequencies in PDE solving.
How to Determine the Component Number? Since the number of dominant frequencies un- known apriori, the solution accuracy can be sensitive to the choice of the component number Q. A too small Q can miss important (high) frequencies while too big Q can bring in excessive noisy frequencies. To address this problem, we set large Q (c.g.. 50), initialize the frequency parameters Hy across a wide range, and then optimize the component weights in the log domain. This turns out to be equivalent to assigning each w, Jefferys prior. Specifically, define w, log(u, (,,). Since do not place an additional prior over Wap we can view p(m,)
x 1. Then we have
d duy
p(w),) = p(Wa)
x
(10)
Wlg
The Jeffreys prior has a very high density near zero, and hence induces strong sparsity during the learning of , (Figueiredo) 200T). Accordingly, the excessive frequency components can be automatically pruned, and the learning drives the remaining My 's toward the target frequencies. This
have been verified by our experiments; see Fig. 4 in Section G
GP Solver Model to Enable Fast Computation. To fulfill efficient and scalable calculation, we multiply our covariance function at each input dimension to construct a product kernel,
(ox))x), = n(x, x110) = II j=1 21,(
(1I)
where 0, {wg Mg: Pa) and θ =00 1=1 are the kernel parameters. Note that the product kernel is equivalent to performing a (high-dimensional) feature mapping for each input dimension and then computing the tensor-po passort across the features. It is a highly expressive structure and commonly used in finite element design (ARNOLD et al. (012).. Next, we create a grid on the domain Ω. We can randomly sample or specially design the locations at each input dimension, and then construct the grid through a Cartesian product. Denote the locations at each input dimension j by h2 = [hy,,-..,
have an M1 X .... x Ma grid,
G= h1 x ... x h;= (x= (==xxxx=.......
E h,, si,
(12)
We will use the grid points on the boundary 02 to fit the boundary conditions and all the grid points as the collocation points to fit the equation.
Denote the solution values at G by u {u(x)]x E 9). which is an Mi x...x Ma array. According to the GP prior over u(-), we have a multi-variate Gaussian prior distribution, p(2f) = N (vec(2)) 0. C). where vee(-) is to flatten U/ into a vector, Cs the covariance matrix computed from G with kernel r(..)). Denote the grid points on the boundary by = G G AQ. To fit the boundary condition, we use Gaussian likelihood, p(g (us) = NV [p. T1 II). where g = vec ([g(x)[x E B)), up are the values of l on B (flatten into a vector), and Tj > 0 is the inverse variance. Next, we want to fit the equation at g. To this end, we need to first obtain the prediction of all the relevant derivatives of u in the PDE, e.g.. 0z, and 0z, u, the grid G. Since u's derivatives also follow the GP prior, the kernel derivative to obtain their cross covariance function (see conditional mean (conditioned on 2/) the prediction. Take 0z, and Day examples. We have
33), with which to compute the GP
oz u(x) = 0n k(x,
Öz u(x) = Ôn, k(x, vec(4))
(13)
where k(x, 9) [k(x, x(),... .., k(x, xhr)] where M = II, Mj and all x constitute G. We can accordingly predict the values of the all the relevant u derivatives G, and combine them to obtain
4

---

Published as a conference paper at ICLR 2024
the PDE (see (T) evaluation at G. which we denote by H. To fit the GP model to the equation, another Gaussian likelihood, p(0 =
f 72 'I, where 0 is an virtual observation, and T2 0. Note that use the same framework as in (Chen et al. 2021 Long al. 2022b). However, there are two critical differences. First, rather than randomly sample the collocation points, we place all the collocation points on grid. Second, rather than assign multivariate Gaussian distribution over the function values and all of its derivatives, we only model the distribution of the function values (at the grid). We then use the GP conditional mean to predict the derivative values. As we will discuss in Section 4, these modeling strategies, coupled with the product covariance (II) enable highly efficient and scalable computation, yet do not need any low rank approximations.
4 Algorithm
We maximize the log joint probability to estimate U, the kernel parameters θ and the likelihood inverse variances T₁ and 72. To ficxibly adjust the influence of the boundary likelihood so as to balance the competition between the boundary and equation likelihoods (Wang et al., (2020a c). we introduce a free pyper-parameter A > O, and maximize the weighted log joint probability,
L(U, O 72) log.No(vee (vec(u))o,, C) + À (gguN( T1
I) + ( (H),, 11)
1 log CI 12-22eed vec(()) CC C (-vvvc(()) + A
1
No 2
T1 ||u -gll=
log TI
2
2
M
T2 ||| (vvv + const.
+ log T2
(14)
2
2
Naive computation of C is extremely expensive when the grid is dense, namely, M is large. That is because the covariance matrix C is between all the grid points, of size M x M (M III M,).
Also, to obtain H. we need to compute the cross-covariance between every needed derivative in the PDE and u across all the grid points. Consequently, the naive computation of the log determinant and inverse of C (see (14) ) and the required cross covaricovee take the time and space complexity (MM3) and (MM2), respectively, which can be infeasible even when each M, is relatively small. For example, when d = 3, and M1 M2 M3 100, we have M 10 and the computation of C will be too costly to be practical (on most computing platforms).
Thanks to that (1) our prior distribution is only on all the function values the grid, and (2) our covariance function is a product over each input dimension (see (TI). We can derive a Kronecker product structure in C, namely, C C₁ X ... & Ca. where C, = (hy, h,) is the kernel matrix on hj the locations at input dimension j. of size M, x Mj. Note that we can also use kaM in (9) Using the Kronecker product properties (Minka 2000), we obtain
M
log [C==> -
log Cy,,
j=1 M, = -=)voc(41=-C) .... @ C vec(U)) = vee (4 X X2..2...X. X2 C77).
(15)
where Xy the tensor-matrix product mode j. Accordingly, we can first compute the local log determinant and inverse at each input dimension (i.e., for each C2), which reduces the time and space complexity to O j=1 Mt) and O(E Mf), respectively. Then we perform the multilinear operation in the last line of (15) i.e., sequentially multiplying the array / with each C- which
takes the time complexity O (EEg M2)M) The computational cost is substantially reduced.
Furthermore, since our product covariance function is factorized over each input dimension, the cross covariance between any derivative of and u itself still maintains a product form- because only the kernel(s) at the corresponding input dimension(s) need to be differentiated. For example,
R(x, cov(O u(x), u(x')) ,x,,, x') Ozyez II, x(Iz,2²))
= Oz K0
j I 21),,
(16)
We found that performing posterior inference over 24 and other parameters, e.g., via variational inference. will degrade the solution accuracy, which partly he hecause the inference and optimization is much more complicated and challenging.
5

---

Published as a conference paper at ICLR 2024
Accordingly, we can also obtain Kronecker product structures in predicting each derivative of u. Take 22, u as an example. According to (33) we can derive that
Der ux) = (02, k( h1) @ o2, k 22. h2) ..... @ k(za, ha)) (CT's &.. vec(())
vec(()
(0,, k(1,, h1 C oz, kk ha) C @ @ k(ra, ha) Cz vec(U)
vec (u x1 Om, k(r1, h1 CT X2 aay k(r2, b2) X3 k(23, ha) X4... X k(rg, C7
Denote the values of Oz u at the grid M by oz dzjag _ (0 "(x) M). Then it is straight- forward to obtain 02,2,44 == =1 u 11 V X2 VIC2 where 1 means taking the derivative w.r.t the first input variable, and we have V1C1 = [0 k(h11, hı): Ohin, k(hIM, h1)l and V,C₂ [o k(h21, h2); Ohaar, k(h2Mg, h2)]- Hence, we just need to perform two tensor- matrix products, which takes (((M) + M2)M) operations, and is efficient and convenient. Similarly, we can compute the prediction of all the associated u derivatives in the PDE operator, with which we can obtain H - the PDE evaluation at the grid in (4. We can then use automatic differentiation to calculate the gradient to maximize (14).
Algorithm Complexity. The time complexity of our algorithm is O M3 + (E, M))M). The space complexity is O M3 + M), including the storage of the covariance matrix at each input dimension, and the solution estimate at grid G, namely U.
5 Related Work
Although the PINN has many success stories, e.g., (Raissi et al. 2020 Chen et al., 2020 Jin et al.
2021 Sirignano and Spiliopoulos) 2018: Zhu et al. 2019:] Geneva and Zabaras 2020 Sahli Costabal et al. 2020), the training is known to be challenging. which is partly due to that applying differential
operators over the NN can complicate the loss landscape (Krishnapriyan et al. 2021). Recent works have analyzed common failure modes of PINNs which include modeling problems exhibiting high- frequency, multi-scale, chaotic, turbulent behaviors (Wang et al. 2020c b a 2022), when the governing PDEs are stiff (Krishnapriyan et al., 2021; Mojgani et al. 2022). One class of approaches to mitigate the training challenge to set different weights for the boundary and residual loss terms. For example, Wight and Zhao (2020) suggested to set a large weight for the boundary loss to prevent the dominance of the residual loss. [Wang et al. (2020a) proposed dynamic weighting scheme based on the gradient statistics of the loss terms. Wang et al. (2020c) developed an adaptive weighting approach based on the eigen-valuee of NTK. Liu and Wang (2021) employed mini-max optimization and updated the loss weights via stochastic ascent. McClenny and Braga-Neto (2020) used a multiplicative soft attention mask to dynamically re- re-weight the loss term on each data point and collocation point. Another strategy is to modify the NN architecture so as to exactly satisfy the boundary conditions, e.g., (Lu et al. [021] Lyu et al.] 2020: Lagaris et al., 1998). However, these methods are restricted to particular types of boundary conditions, and are less flexible than the original PINN framework. Tancii et al. (2020); Wang et al. (2021b) used Gaussian distributions to construct random Fourier features to improve the learning of the figh-freueecy and multi-scale information. The number of Gaussian variances and their scales are critical to the success of these methods. But these hyperparameters are quite difficult choose.
Earlier works (Graepel 2003) have used GP for solving linear PDEs with noisy measurement of source terms. In (Wang et al.) 2021a), the rationale and guarantees of using GP as a prior for PDE solutions are discussed. The work also justifies the usage of the product kernel in terms of sample path properties. The recent work (Chen et al. 2 develops general approach for solving both linear and nonlinear PDEs. Long et al. (202bb) proposed a GP framework to integrate various differential equations. The recent work (Chen et al. (023) uses sparse inverse Cholesky factorization to approximate the kernel matrix so as to handle a large number of collocation points. These methods use SE and Matérn kernels and are challenging to capture high- gegh-freueecc and multi-scale solutions. The recent work (Pööttner et al., 2022) proposes a infoomed GP solver for linear PDEs that generalizes weighted residuals. In (Härkönen et al. 2022). a GP kernel is constructed via the Ehrenpreis- Pirenpreis-Palamodoo fundamental principle and nonlinear Fourier transform to solve linear PDEs with constant coefficients. This work also derives the spectral mixture kernel as an instance of its own kernel design. The computational advantage of using Kronecker product structures have been realized in [Saateci] [0022), and applied in other tasks, such as nonparametric tensor decomposition (Xu et al. 2012). sparse approximation with massive inducing points (Wilson and Nickisch 2015 Izmailov

---

Published as a conference paper at ICLR 2024
et al., 2018), and high- high-dimensiooal output regression (he et al. 2019) In Wilson et al. 2015) it further points out that if one uses a regular (evenly-spaced), each kernel matrix will has a Toeplitz structure, which can lead to (n log n)) computation. However, in machine learning applications, data is typically not observed at a grid and the Kronecker product has a limited usage. By contrast, for PDE solving, it is natural to estimate the solution values on a grid, which opens the possibility of using Kronecker products combined with GP for efficient computation. More general discussions about Bayesian learning and PDE problems are given in (Owhadi 2015 Cockayne et al. 2017). Tensor methods used in numerical computation are discussed in (Gavrilyuk and Khoromskij 2019
6 Experiment
To evaluate GP-HM, we considered three commonly-used benchmark PDE families in the literature of machine learning solvers (Raissi et al., 2019; Wang et al.) 2021b Krishnapriyan et al. [021]: Poisson, Allen-Calle and Advection. Following the prior works, we fabricated series of solutions to thoroughly examine the performance. The details are given in Section B Appendix.
We compared with the following state-of-the-ar art ML solvers: (1) standard PINN, (2) Weighted PINN (W-PINN) that up-weight the boundary loss to reduce the dominance of the residual loss, and to more effectively propagate the boundary information, (3) Rowdy (Jagtap et al. 2022), PINN with an adaptive activation function, which combines a standard activation with several sin or cos activations.
(4) RFF-PINN, feeding Random Fourier Features to the PINN (Wang al. 2021b). To ensure RFF- PINN to achieve the best performance, we followed (Wang et al. 2020c) to dynamically re- re-weight the loss terms based NTK eigenvalues (Wang et al. 2020c). (5) Spectral Method (Boyd 200T), which approximates the solution with a linear combination of trigonometric bases, and estimates the basis coefficients via least mean squares. In addition, we also tested (6) GP-SE and (7) GP-GGtémn GP solvers with the square exponential (SE) and the Matérn kernel. The details about the hyperparameter setting and tuning is provided in Section B of Appendix. We denote our method using the covariance function based on and (9) by GP-HM-StM and GP-HM-GM, respectively. We compared with several traditional numerical solvers: (8) Chebfur that solves PDEs based on Chebyshev interpolants, (9) Finite Difference (FD), which solves the PDEs via discretization based on finite difference. We used PyPDE library to solve ID/2D Poisson equations, and ID advection (using methods of lines). Note that PyPDE does not support solving nonlinear stationary PDEs, namely 1D/2D Allen Cahn Equation in (8). and so we implemented the finite difference with Seipy and Krylov method for root finding. Note also that the Chebfun library does not support 2D Poisson and nonlinear stationary PDEs, namely. DD2D Allen-Cahn equation, and so it has very limited usage. We employed the default settings in Chebfun library. When using PyPDE, we set spacial discretization 400 and 400 time steps (if needed). For ID Allen-cahn, the spatial discretization is set to 400, For 2D Allen-cahn, we used a 45 x 45 grid; otherwise, the root finding either ran forever or failed due to numerical instability. We have also tested the Spectral Galerking method implemented by the Shenfun library F However, we found failed every test case (the relative L22 error is at several thousands). Hence, we did not report the results.
Solution Accuracy. We report the relative L2 error (normalized root-mean- square error) of each method in [an The best result and the smaller error between GP- HM-M and GP-HM- GM are made bold. We can see that, among all the ML solvers, our method achieves the smallest solution error in all the cases except that for the ID Poisson equation with solution U2 RFF-PINN is better. However, in all the cases, the solution error of GP-HM achieves at least 1e-3 level. In quite few cases, our method even reaches error around 1 and 1e-7. It shows that GP-HM can successfully solve all these equations. By contrast, GP solvers using the plain SE and Matérn kernel result in several orders of the magnitude bigger errors. The standard PINN and W-PINN basically failed to solve every equation. While Rowdy improved upon PINN and W-PINN in most cases, the error is still quite large. The inferior performance of the spectral method implies that only using trigonometric bases is not sufficient. With the usage of the random Fourier features, RFF-PINN can greatly boost the performance of PINN and W-PINN in many cases. However, in most cases, is still much inferior to GM-HM. The performance of RFF-PINN is very sensitive to the number and scales of the Gaussian variance, and these hyper-parametes are not easy to choose. We have tried 20
https://ww..chbffftt
7

---

Published as a conference paper at ICLR 2024
ID
2D
Method
M1
U2
Ug
04
Ug
Uk
Ut
PINN
1.36e0
1.40e0
1.00e0
1.42e1
6.03e-1
1.6e0
9.99e-1
W-PINN
1.31e0
2.65e-1
1.86e0
2.60e1
6.94e-1
1.63e0
6.75e-1
RFF-PINN
4.97c-4
2.00c-5
7.29c-2
2.80c-1
5.74c-1
1.69c0
7.99 c-1
Rowdy
1.70c0
1.00c0
1.00e0
1.01e0
1.03e0
2.24e1
7.36c-1
Spectral method
2.36e-2
3,47c0
1.02c0
1.02c0
9,98e-1
1.58c-2
1.0c00
Chehfun
3.05e-11
1.17e-11
5.81e-11
1.14c-10
8.95c-10
N/A
N/A
Finite Difference
5.58e-1
4.78e-2
2.34e-1
1,47c0
1.40e0
2.33e-1
1.75e-2
GP-SE
2.70e-2
9.99e-1
9.99e-1
3.19e-1
9.75e-1
9,99e-1
9.53e-1
GP-Matém
3.32e-2
9.8e-1
5.15e-1
1.83e-2
6.27e-1
6,28e-1
3.54e-2
GP-HM-GM
3.99e-7
2.73e-3
3,92c-6
1.55c-6
1.82c-3
6.46e-5
1.06c-3
GP-HM-SIM
6.53e-7
2.71e-3
3.17e-6
8.97e-7
4.22e-4
6.87e-5
1.02e-3
Table Relative L2 error in solving ID and 2D Poisson equations, where uy's are different high- frequency and multi-seale solutions: U1 sin(100x), uy sin(x) + 0.1.1 (20() + 0.05 co05coss0000 (100z))) us s6r( cos(100z), u4 = sin(200z)) us (500r) 2(x 0.5)² ue = sin(100z) sin(100y) and
u7 sin(6x) sin()0z) sin(6y) sin (20y).
ID Allen-cahn
Method
2D Allen-cahn
1D Advection
U1
U2
PINN
1.41e0
1.14el
I,9e1
1.00e0
W-PINN
1.34e0
1.45el
2.03e1
1.01e0
RFF-PINN
1.24e-3
2.46e-1
7.17e-1
9.96e-1
Rowdy
1.30c0
1.31c0
1.18c0
1.03c0
Spectral method
2.34c-2
2.45c1
2.45c1
2.67c0
Chebfun
1.39e-08
2.94e-10
N/A
1.39c0
Finite Difference
2.32e-01
2.36e-1
3.23e0
1.29e-1
GP-SE
2.74c-2
1.06e-2
3,48c-1
9.99c-1
GP-Matém
3.32c-2
5.160-2
2.960-1
9,99c-1
SP--HM-SM
7.71c-6
4.76c-6
2.99e-3
9.08e-4
GP-HM-GM
4.91c-6
4.24c-6
5.78e-3
3.59e-3
Table Relative L2 error solving ID, Allen-cahn equations and advection equation, where un and
U2 are two test solutions for ID Allen-cahn: i1 = sin(100r), up = sin(6z) cos(100x). The test solution for 2D Allen-cahn is (x)) + 0.1 s20() + cos(100z)) ((iin))))(((())) + 0.1 sin(20y) + cos(1 (os0)) and for ID advection equation is sin (x 200r).
settings and report the best performance (see Section B in Appendix). Compared with traditional solvers, we can see Chebfun performs very well, and achieves the highest solution accuracy except for the ID advection problem. However, Chebfun is limited to 1D problems and temporal PDEs. cannot handle 2D stationary PDEs, no matter linear or nonlinear. Finite Difference can provide reasonable accuracy, but the performance is consistently much worse than GP-HM. This might be due to the challenge in solving the root finding problem, caused by the information implied in the source term. Overall, can that method is general enough to solve different types of PDEs (DD/2D, linear/nonlinear, stationary and non-stationary): to achieve satisfactory accuracy,
scale
do NOT need to change the computation framework re-develop the solver. By contrast, it is known that the success of numerical solvers tightly binds to the specific problem, domain knowledge, skillful implementation, and numerous numerical tricks. Any change of these aspects can cause failures of the solvers and demand for a re-design and re-implementation. therefore brings significant challenges in usage.
We then show the point-wiis wise solution error in Fig. I and in Appendix Fig. 7 We can see that GP-SE difficult capture high frequencies. While GP-Matérn is better, is unable to grasp all the scale information. RFF-PINN successfully captured multi-scale frequencies in Fig. M but it failed in more challenging cases as in Fig. 2 and In 2D Poisson and 1D Advection, the point-wise error of both GP-HM-StM and GP-HM-GNM is quite uniform across the domain and close to zero (dark blue); see Fig. 3 and in Appendix Fig. 6 7 By contrast, the other methods exhibit large errors in few local regions. These results show that GP-HM not only gives superior global accuracy, but locally recovers individual solution values.
Point-wise Error.
Frequency Learning. Third, investigated the learned component weights Wy and frequencies H4 of GP-HM. In Fig. 4 we show the results for two Poisson equations. As we can see, although the number of components is set to be much larger than the number of true frequencies, the estimation
8

---

Published as a conference paper at ICLR 2024
Reound buth
Prediction Boundary
www
u00
GP-HM-StM
GP-HM-GM
GP-SE
GP-Matern
RFF-PINN
Figure Prediction for the 1D Poisson equation with solution sin(x) (r + 0.1 sin (20x)) + 0.05 cos(100zz).
03 GP-HM-StM
0,5 GP-HM-GM
0.5 GP-SE
05 GP-Matern
05 RFF-PINN
Figure 2: Prediction for the 1D) Poisson equation with solution s5000)) - 2(x 0.5)
of most weights Wy is very small (less than 10-10). That means, excessive frequency components have been automatically pruned. The remaining components with significant weights completely match the number of true frequencies in the solution. The frequency estimation (4g is very close to the ground-truth. This demonstrates that the implicit Jefferys prior (by optimizing wy in the log space) can indeed implement sparsity, select the right frequency number, and recover the ground-truth frequency values. Finally, we show additional results in Section C of Appendix.
1.0
Ground-Truth
GP-HM-StM
GP-HM-GM
GP-SE
GP-Matern
RFF-PINN
Figure Point-wise solution error for Allen-cahn equation,
the
solution
(s) + 0.1 sin(20z) + cos(100z)) (sin(y) - 0.1 sin()0g) + cos( 10s0))
u(x, = si.xx/issn0000x///i11100////
uix,y) - sim(100xssin11000)
0.3
106.002 (106)
99.957(100)
202
100.000(100)
Weight (wgg)
93.999 (94)
01:
163
Index of freq components
Index of freq components
Index of freq components
(a) Poisson-1D
(b) Poisson-2D z-dim
(c) Poisson-2D y-dim
Figure 4: The learned component weights and frequency values. For each number pair a(b) the figure, "a" the learned frequency by GP-HM, and is the ground-truth. The expressions on the top are the solutions.
7 Conclusion
We have presented GP-HM, GP solver specifically designed for high-frequency and multi-scale PDEs. On set of benchmark tasks, GP-HM shows promising performance. This might motivate alternative directions of developing machine learning solvers. In the future. plan to develop more powerful optimization algorithms to further accelerate the convergence and to investigate GP-HM in variety of practical applications.

---

Published as a conference paper at ICLR 2024
Acknowledgments
This work has been supported by MURI AFOSR grant FA9550-20-1-0358, NSF CAREER Award IIS IIS-2046295, and and NSF DAC-2311685.
References
ARNOLD, D. N., BOFFI, D., and BONIZZONI. F. (2012). Tensor product finite element differential forms and their approximation properties. arXiv preprint arXiv: 121.6622
Bishop, C. M. (2007). Pattern Recognition and Machine Learning. Springer.
Boyd, J. P. (2001). Chebyshev and Fourier spectral methods. Courier Corporation.
Chen, Y., Hosseini, B., Owhadi, H., and Stuart, A. M. (2021). Solving and learning nonlinear PDEs with Gaussian processes. arXiv preprint arXiv:2103.12959
Chen, Y., Lu, L., Karniadakis, G. E., and Dal Negro, L. (2020). Physics-informen neural networks for inverse problems in nano-optics and metamaterials. Optics express, 28(8): 1618-11633.
Chen, Y., Owhadi, H., and Schäfer, (2023). Sparse cholesky factorization for solving nonlinear pdes via gaussian processes. arXiv preprint arXivv2304.01294.
Cockayne, J., Oates, C., Sullivan, T., and Girolami, M. (2017). Probabilistic numerical methods for pde-constrained bayesian inverse problems. In AIP Conference Proceedings, volume 1853, AIP Publishing.
Figueiredo, M. (2001). Adaptive sparseness using Jeffreys prior. Advances in neural information processing systems, 14.
Frostig, R., Johnson, M. J., and Leary, C. (2018). Compiling machine learning programs via high-level tracing. Systems for Machine Learning. 4(9).
Gavrilyuk, I. and Khoromskij. B. N. (2019). Tensor numerical methods: actual theory and recent applications. Computational Methods Applied Mathematics, 19(1):1-4.
Geneva, N. and Zabaras, N. (2020). Modeling the dynamics of pde systems with physics-constrained deep auto-regressive networks. Journal of Computational Physics, 403:109056.
Graepel, T. (2003). Solving noisy linear operator equations by gaussian processes: Application to ordinary and partial differential equations. ICML, volume 3, pages 234-241.
Grimmett, G. and Stirzaker, D. (2020). Probability and random processes. Oxford university press.
Härkönen, M., Lange-Hegermann M., and Raictua, B. (2022). Gaussian process priors for systems of linear partial differential equations with constant coefficients. arXiv preprint 2rXiv:2212.14331
Izmailov, P., Novikov, A., and Kropotov, D. (2018). Scalable gaussian processes with billions of inducing inputs via tensor train decomposition. In International Conference Artificial Intelligence and Statistics, pages 726-735.
Jagtap, A. D., Shin, Y., Kawaguchi, K., and Karniadakis, G. E. (2022). Deep kronecker neu- ral networks: A general framework for neural networks with adaptive activation functions. Neurocomputing, 468: 165-180.
Jin, X., Cai, S., Li, H., and Kariiadakis, G. E. (2021). Nsfnets (navier-stokes flow nets): Physics- informed neural networks for the incompressible navier-stokes equations. Journal of Computational Physics, 426: 10991..
Khintchine, A. (1934). Korrelationstheorie der stationären stochastischen prozesse. Mathematische Annalen, 109(1 (19(1):00 6164
Krishnapriyan, A., Gholami. A., Zhe, S., Kirby, R., and Mahoney, M. W. (2021). Characterizing possible failure modes in physics-infomee informed neural networks. Advances in Neural Information Processing Systems, 34:26548-2653
10

---

Published as a conference paper at ICLR 2024
Lagaris, I. E.,, Likas, A., and Fotiadis, D. I. (1998). Artificial neural networks for solving ordinary and partial differential equations. IEEE transactions on neural networks, 955:987-11000.
Lathi, B. P. (1998). Modern digital and analog communication systems. Oxford University Press, Inc.
Liu, D. and Wang. Y. (2021). A dual-dimer method for training physics-constrained neural networks with minimax architecture. Neural Networks, 136:112-125.
Long, D., Wang, Z.,, Krishnapriyan, A., Kirby, R., Zhe, S., and Mahoney, M. (2022a). AutolP: A united framework integrate physics into Gaussian processes. In International Conference on Machine Learning, pages 14210-14222. PMLR.
Long, D., Wang, Z., Krishnapriyan, A., Kirby, R., Zhe, S., and Mahoney, M. (2022b). AutolP: A united framework to integrate physics into Gaussian processes. In Chaudhuri, K., Jegelka, S., Song, L.,, Szepesvari, C., Niu, G., and Sabato, S., editors, Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 14210-14222, PMLR.
Lu, L., Pestourie, R., Yao, W., Wang, Z., Verdugo, F., and Johnson, S. G. (2021). neural networks with hard constraints for inverse design. SIAM Journal on Scientific Computing, 43(6):B1 105-B1 1332.
Lyu, L.. Wu, K., Du, R., and Chen, J. (2020). Enforcing exact boundary and initial conditions in the deep mixed residual method. arXiv preprint arXiv:2008.0149.
McClenny, L. and Braga-Neto, U. (2020). Self-adaptive physics-informed informed neural networks using soft attention mechanism. arXiv preprint arXiv:2009.04544 004544.
Minka, T. P. (2000), Old and new matrix algebra useful for statistics. See www. stat. cmu. edu/minkaa/ppper//maarix html, 4.
Mojgani, R., Balajewicz, M., and Hassanzadeh, P. (2022). Lagrangian pinns: causality-conforminn
solution to failure modes of physics-informed neural networks, arXiv preprint arX2v:2205.02902.
Owhadi, H. (2015). Bayesian numerical homogenization. Multiscale Modeling & Simulation, 13(3):812-828 828.
Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al. (2019). Pytorch: An imperative style, high-performance deep learning library. Advances in neural information processing systems, 32.
Pförtner, M., Steinwart, I., Hennig, P., and Wenger, J. (2022). Physics-informed gaussian process regression generalizes linear pde solvers. arXiv preprint arXiv:2212.12474
Rahaman, N., Baratin, A., Arpit, D., Draxler, F., Lin, M., Hamprecht, F., Bengio, Y., and Courville, A. (2019). On the spectral bias of neural networks, In International Conference on Machine Learning. pages 5301-5310. PMLR.
Raissi, M., Perdikaris, P., and Karniadakis, G. E. (2017). Machine learning of linear differential equations using gaussian processes. Journal of Computational Physics, 348:683-693.
Raissi, M., Perdikaris, P., and Karniadakis, G. E. (2019). deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707.
neural networks: A
Raissi, M., Yazdani, A., and Karniadakis, G. E. (2020). Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations. Science, (67(6481)) 1026 1030.
Rasmussen, C. E. and Williams, C. K. (2006). Gaussian Processes for Machine Learning, MIT Press.
Saateci, Y. (2012). Scalable inference for structured Gaussian process models. PhD thesis, Citeseer.
11

---

Published as a conference paper at ICLR 2024
Sahli Costabal, F., Yang, Y., Perdikaris, P., Hurtado, D. E., and Kuhl, E. (2020). Physics-informed neural networks for cardiac activation mapping. Frontiers in Physics, 8:42.
Sirignano, J. and Spiliopoulos, K. (2018). Dgm: A deep learning algorithm for solving partial differential equations. Journal of computational physics, 375:1339-136..
Tancik, M., Srinivasan, P., Mildenhall, B., Fridovich-Keil, S., Raghavan, N., Singhal, U., Ramamoor- thi, R., Barron, J., and Ng. R. (2020). Fourier features let networks learn high frequency functions in low dimensional domains. Advances Neural Information Processing Systems, 33:7537-7547.
Wang, J., Cockayne, J., Chkrebtii, O., Sullivan, T. J., and Oates, C. J. (2021a). Bayesian numerical methods for nonlinear partial differential equations. Statistics and Computing, 31:1-20.
Wang, S., Sankaran, S., and Perdikaris, P. (2022). Respecting causality is all you need for training physics-informed neural networks. arXiv preprint arXiv:2203.07404
Wang, S., Teng, Y., and Perdikaris, P. (2020a). Understanding and mitigating gradient pathologies in physics-informed neural networks. arXiv preprint arXiv:2001.044 .04536.
Wang, S., Wang, H., and Perdikaris, (2020b). On the eigenvector bias of fourier feature networks: From regression to solving multi-scale pdes with physics- -infries neural networks. arXiv preprint arXiv:2012.10047.
Wang, S., Wang. H., and Perdikaris, P. (2021b). On the eigenvector bias of fourier feature networks: From regression to solving multi-scale pdes with physics-informed neural networks. Computer Methods in Applied Mechanics and Engineering, 384:1 13938.
Wang, S., Yu, X., and Perdikaris, P. (2020c). When and why pinns fail to train: A neural tangent kernel perspective. arXiv preprint arXiv:2007.14527
Wiener, N. (1930). Generalized harmonic analysis. Acta mathematica, 55(1):117-258.
Wight, C. L. and Zhao, J. (2020). Solving allen-cahn and cahn-hilliard equations using the adaptive physics informed neural networks. arXiv preprint arXiv:2007. 244444..
Wilson, A. and Adams, R. (2013). Gaussian process kernels for pattern discovery and extrapolation. In International conference on machine learning, pages 1067-1075. PMLR.
Wilson, A. and Nickisch, H. (2015). Kernel interpolation for scalable structured gaussian processes (kiss-gp). In International Conference on Machine Learning, pages 1775-1784.
Wilson, A. G., Dann, C., and Nickisch, H. (2015). Thoughts on massively scalable gaussian processes. arXiv preprint 151v:1511.015
Xu, Z.. Yan, F., and Qi. (2012). Infinite tucker decomposition: nonparametric bayesian models for multiway data analysis. In Proceedings of the 29th International Coference on International Conference on Machine Learning, pages 1675-1682.
Zhe, S., Xing, W., and Kirby, R. M. (2019). Scalable high-order gaussian process regression. In The 22nd International Conference Artificial Intelligence and Statistics, pages 2611-2620.
Zhu, Y., Zabaras, N., Koutsourelakis, P.-S., and Perdikaris, P. (2019).
deep learning for high-dimensional surrogate modeling and uncertainty quantification without labeled data. Journal of Computational Physics, 394:56-81.
12

---

Published as a conference paper at ICLR 2024
Appendix
A Covariance Function Derivation
In this section, we show how to obtain our covariance function in (8) of the main paper. We leverage the fact that the student density a scale mixture of Gaussians with Gamma prior over the inverse variance,
(, = =f6 N(olu, T (aln,TT)GGam G)am (ra,((r,,))))) bjdr
ba
1/2
b (x- a- 1/2
I(a + 1/2)).
(17)
T(a)
2
The key to obtain this is to leverage the form of the normalizer of the Gamma distribution. When merging terms in the Gaussian and Gamma prior in the integration, one can construct another unnormalized Gamma distribution. Accordingly, the integration w.r. gives rises to the normalizer.
we set v = 2a and À = a/b, we immediately obtain the standard student t density,
Staal X, v) T(v/2+ 1/2)
A
1/2
X(x
p)²² 1/21/2
(18)
T(v/2)
TV
where j the mean, the precision (inverse variance) parameters, and is the degree of freedom.
Next, we observe that the spectral density of a Matérn covariance function is a student & density (Ras mussen and Williams 2006). Given the Matérn covariance
21-1 [(v)
222 (V2y'=>r)) 1
K. (v r
x'
Topp, (x,
(19)
Pg
Pg
the spectral density is St(s; 0, p²² 2v). 2v)). That means,
Te.p (A)
St(s; 0, 4 2²²22,2v) exp (i2Ts A)ds,
(20)
where A = 1x - x'.. obtain
From the scale-mixture form
we can set â = v and b = 0/(4x2 p²), and
St(s; O, 4n 4x²22,22)) = / N(s|0, - V)sl0, Gam(rr 8)dr.
(21)
Substituting (1 into 20)
we have
Te.p A) [ Gam(rla, b) / 7(s00,,=))xxp(e A Addddr
(22)
Consider the inverse Fourier transform,
St(s; 4x 2v) exp(i2n s)ds
J Gam(rla, b) I N(ssu (1) (i2ms dsdr,
(23)
we observe that
F- W 1 = / N(slu, (0)),,,--))xpp (i2ms A)
=exp 2 exp (12mp A)
=F-'(W --1 exp (12mp A)
= =/<(s/0,7>>) ( (12TS A) ds exp (12mp- A)
(24)
13

---

Published as a conference paper at ICLR 2024
where F-1 is the inverse Fourier transform, and indicates complex numbers. Note that when set μ = 0, from the second line, we see FF "F = exp (-24 A²). That means, the
inverse transform just moves out a Fourier basis with frequency μ..
Substitute
into
wc obtain
Mt 4n p², 22) eexp(1274 s)ds
Gam(rla, / ((s|0,==)) exp (i2Ts A) dsdr exp(i2mm A) (A) exp(i2 A).
Therefore, when we model the spectral density S(s) as a mixture of student-t distribution,
Q
S(s) We (St(s; Mqp 4n² pg. 2v) + St(s; MMy 4m² Py: 2v))
(25)
It is straightforward to obtain the following covariance function,
ASSI (a, x.)
Σ Wy Te,p, (x, cos((2t - -x')) x²))-
(26)
B Experimental Settings
The Poisson Equation. We considered 1D and 2D Poisson equations with different source functions that lead to various scale information in the solution. We used Dirichlet boundary conditions.
Ugr [(x), x E [0,
Uga + Myy = f(x,)), (c, E [0,2T] x [0,2m].
(27)
For the ID Poisson equation, we created source functions f that give the following and multi-frequency solutions, #1 = sin 100z). U2 sin(x) + 0.1 sin + 0.05 cos( 100x), M4 = sin (6n) cos(1 100x),, and U4 = sin sin (200z). addition, we tested with a challenging hybrid solution that mixes high-frequee freeueerruecyyy with quadratic function, Us = sin(500x) - 2(x 0.5)² where
we set [0, 1]. For the 2D Poisson equation, we tested with the following multi-scale solutions, sin(62) (20xx) + sin(Gy) sin(20y).
Mg sin( 110zz) sin(100y) and
Allen- Allen-Cahn Equation. We considered 1D and 2D Allen-Cahn (nonlinear diffusion-reaction equa- tions with different source functions and Dirichlet boundary conditions.
Ugz + u(uu - 1) = f(x),
[0, 2m],
Map + U My + uuuu 1) = f(x (x,y) E [0, 1] x [0, 1].
(28)
For the ID equation, we tested with solutions ii sin (100r) and U2 sin(6a) cos c1000)).). For the 2D equation, we created the source f that gives the following mixed-scale solution, u = (sin(x:) 0.1 sin(20)) + cos010)z) (sin(y) + 0.1 sin(20y) + cos(100y))
Advection Equation. Third, we evaluated with a ID advection (one-way) equation,
4 + 200μ, 0, x E [0, 2T), t E [0.1].
(29)
We used the Dirichlet boundary conditions, and the solution has an analytical form, u(x, () 1(( 2007) where (xx) is the initial condition for which we chose as l((x)) = sin
Method Implementation. We implemented our method with JAX (Frostig et al. 2018)) while all the competing ML based solvers with Pytorch (Paszke et al. 2019). For all the kernels, we initialized the length-scale to 1. For the Matérn kernel (component), we chose set the number of components Q = 30, and initialized each W. = 1 1/Q. For 1D Poisson and ID Allen-cahn equations, we varied the ID mesh points from 400, 600 and 900. For 2D Poisson, 2D Allen-cahn and ID advection, we varied the mesh from 200 x 200, 400 x 400 and 600 600. We chose ending frequency F from (20, 40, 100), and initialize u, 's with 1 nspace (0, F, Q).
5/2. For our method, we
14

---

Published as a conference paper at ICLR 2024
2m
R GP-HM-GM
2n
" GP-Matern
2n
GP- HM-StM
GP-SE
RFF-PINN
Figure 5: Prediction for the ID Poisson equation with solution x (200xr....
3,000
Ground-Trutee
GP-HM-StM
GP-HM-GM
GP-SE
GP-Matern
RFF-PINN
Figure 6: Point-wise solution error 2D Poisson equation the solution is u(x:) sin sin + s6y( s20y)
We used ADAM for optimization, and the learning rate was set to 10-2, The maximum number of iterations was
1M, and we used the summation of the boundary loss and residual loss less than 10-6 as the stopping condition. The solution estimate II was initialized as zero. We set the A = 500. For W-PINN, we varied the weight of the residual loss from (10, 10 10?). For Rowdy, we combined tanh with sin activation, 6(x) = tanh(x) (x) + Σ sin )) We followed the original Rowdy paper (Jagtap et al. 2022) to set the scaling factor = 10 and varied K from 3, and 9. For the spectral method, we used 200 Trigonometric bases, including cos(n.r and sin(nx) where n 1....... 100. We used the tensor-product for the 2D problems and ID advection. We used the least mean square method to estimate the basis weights. To run RFF-PINNS, we need to specify the number and scales of the Gaussian variances to construct the random features. To ensure broad coverage, we varied the number of variances from {1, 2. 3, 5). For each number, set the variances to be the commonly used values suggested by authors, {1, 20, 50, 100}, combined with randomly sampled ones. The detailed specification is given by Table 3 There are in total 20 settings. We report the best result of RFF-PINN from these settings. For all the PINN based methods, we varied the number of collocation points from 10K and 12K.
Number
Scales
1, 20,0, 50, 100, rand(1. [1, K])
3 rand (2, {1, 20, 30, 100, rand(1, [1, K))H) (0),,,, 2 rand [1, K])
3 x rand (3, {1, 20, 50, 100, r1nd [1, K)}),, x rand (3, [1, K])
2 x {1, 20, 50, 100, rand(1, [1, K]) rand [1, K])
Table 3: The number and scales of the Gaussian variances used in RFF-PINN, where rand(k, A) means randomly selecting elements from the set A without replacement, /x means repeating the sampling to generate I configurations, and the maximum candidate frequency for which we set K 200.
C More Results
C.1 Learning Behavior and Computational Efficiency
We examined the training behavior of our method. As shown in Fig. 8 with the covariance based on the student mixture, GP-HM can converge faster or behave more robustly during the training. Overall, in most cases, GP-HM with covariance based the student mixture performs better than with Gaussian mixture.
The computation efficiency of GP-HM is comparable to PINN- approaches. For example, solving ID Poisson and Allen-cahn equations, the average per-iteration time of GP-HM (mesh 200), PINN and RFF-PINN are 0.006, 0.004 and 0.004 seconds. For 2D Poisson and Allen-cahn equations
15

---

Published as a conference paper at ICLR 2024
Ground Truth
GP-HM-StM
GP-HM-GM
GP-SE
GP-Matern
RFF-PINN
Figure 7:
Point-wise solution error for ID Advection equation and the solution is sin(2 2001).
ID
2D
Method
122
Uy
Us
MG
M
PINN
622
688
624
619
275
RFF-PINN
562
546
576
555
544
3, 3,394
5, 493
Spectral method
502
495
600
480
517
5, 778
7. 060662
Chebfun
1.05
1.22
1.19
1.38
3.90
N/A
N/A
Finite Difference
1.25c-02
1.27c-2
1.22c-2
1.22c-2
1.22c-2
N/A
N/A
GP-HM-GM 536
1,858
775
703
3,510
173
561
GP-HM-SIM 683
164
914
852
4,263
5, 263
6, 435
Table 4: Running time in seconds solving ID Poisson equations, where uy 's are different high- frequency and multi-scale solutions: u1 sin(100x), u2 sin(z) + 0.1 sin(20x) 0.05 co00x))))) uy sin(6x) cos( 100z), u4 sin(200.2)
sin sin(500z)
2(x 0.5)². Sim(100z) sin (100y) and
uy s6n()x) (20r) sin(6y) sin (20y)-
GP HM StM GP-4M-GN
12 error log scale
3
60K Training epoch
40K
100K
200K
400K
690K
B00K
Training epach
(a) 1D Poisson with solution ug.
(b) 2D Poisson with solution ur
Figure 8: The learning curve.
and ID advection, the average per-iteration time of GP-HM (mesh 200 200) is 0.022 seconds while PINN and RFF-PINN (with two scales) took 0.006 and 0.02 seconds, respectively. We examined the running time on Linux workstation with NVIDIA GeForce RTX 3090 GPU. Thanks to the usage of the grid structure and the product covariance, our GP solver can scale to a large number of collocation points, without need for additional low rank approximations.
We also reported the total running time for every test case in Table We can see that the running time of GP-HM is comparable PINN and RFF-PINN in most cases. However, the ML based solvers are slower than traditional methods. This might be because the ML solvers use optimization to find the solution approximation while the numerical methods often use interpolation and fixed point iterations, which are usually more efficient.
C.2 Influence of Collocation Point Quantity
examined how the number of collocation points influences the solution accuracy. To this end, we tested with ID Poisson and 2D Poisson equation, whose solutions include high frequencies. In Fig. we show the solution accuracy with different grid sizes (resolutions). We can see that both PDEs, using low resolutions gives much worse accuracy, e.g., less than 200 in ID and 200 x 200 in 2D Poisson. The decent performance obtained only when resolutions high enough, e.g., 300 ID and 400 400 in 2D Poisson. That means, the number of collocation points is large (particularly for
9
16

---

Published as a conference paper at ICLR 2024
ID Allen-cahn
Method
2D Allen-cahn
ID Advection
1/
02
PINN
509
828
2,509
2,496
RFF-PINN
1,227
1,172
4,421
2,495
Spectral method
504
552
3,840
2,188
Chebfun
6.57
6.0
N/A
1.39
Finite Difference
2.32c-1
2.36c-1
1,130
12.6
GP-HM-SIM
735
2,291
7,447
2,574
GP-HM-GM 612
2,013
6,238
2,239
Table 5: Running time in seconds solving ID. 2D Allen-cahn equations and ID advection equation, where un and Ug are two test solutions for ID Allen-cahn: μ1 sin(100z), U2 sin(6x) cos( 100z). The test solution for 2D Allen-cahn (s))(x) + 0.1 sin(20z) advection equation is sin (x 200t).
(00100000(ss(((((( + 0.1 sin(20g) cos(100y)), and for ID
Log (Relative L, Emror)
Logs: (Nelative La Error)
GP+HM-N GM
GP-M-GM GP-MM-SUM
Resalution
Resolution
(a) ID Poisson with Mg
(b) 2D Poisson with
Figure The solution error using different grid resolutions.
2D problems, e.g., 160K collocation points for the resolution 400 x 400). However, it is extremely costly or practically infeasible for the existent GP solvers incorporate massive collocation points, due to the huge covariance matrix. Our GP solver (defined on grid) and computational method can avoid computing the full covariance matrix, and highly efficiently scale to high resolutions. The results have demonstrated the importance and value of our model and computation method.
D Limitation and Discussion
The learning of GP-HM can automatically prune useless frequencies and meanwhile adjusts jog for the preserved components, namely, those with nontrivial values of wy. align with the true frequencies in the solution. However, the selection and adjustment of the covariance components often require many iterations, like tens of thousands, see Fig. Sa More interestingly, the first-order optimization approaches, like ADAM, perform well, yet the second-order optimization, which in theory converges much faster, such as L-BFGS, performs badly. This might be because the component selection and adjustment is a challenging optimization task, and might easily encounter inferior local optimums. To overcome this limitation and challenge, we plan to try with alternative sparse prior distribution over the weights Wg. such as the horse-shoe prior and the spike-and-slab prior, to accelerate the pruning and frequency learning. We also plan to try other optimization strategies, such alternating updates of the component weights and frequencies, to see if we can accelerate the convergence and if we can embed and take advantage the second-order optimization algorithms.
found that
17


Reviews:
[{"review_id": "V2Nf468K96", "reviewer": "Reviewer_GK3t", "summary": "This work uses Gaussian processes to solve partial differential equations. The authors propose to use a spectral mixture kernel and learn the mixture weights from data with a sparsity-inducing prior. To achieve scalability, they place collocation points on a grid and assume a product kernel which induces Kronecker structure in the resulting covariance matrices. The approach is then evaluated on three common model PDEs in one and two dimensions.", "questions": "- What about the marginal uncertainty output by the Gaussian process? Does it capture the approximation error of the mean or is it miscalibrated? \n- What's the impact of the number of components of the mixture kernel on the approximated solution? An ablation such as a plot of number of components vs error would be informative. Also how does it affect the difficulty of the optimization problem? I could imagine choosing too many components makes the optimization problem significantly harder to solve in practice.\n- Do the kernel matrices you define have Toeplitz structure? Stationary kernels on a 1D grid would have. This could further accelerate the required computations (see Section 3.2. of Wilson et al. 2015).", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Overall, the paper presents a promising idea for how to solve certain classes of PDEs with Gaussian processes. The idea of using a sparsity-inducing prior over the frequencies of a spectral mixture kernel seems very useful and the paper demonstrates the robustness and applicability of the approach.", "weaknesses": "There are some weaknesses in the paper regarding the claimed novelty of the approach, the experimental evaluation and the proper attribution of existing approaches.\n\n## Originality\nIn my opinion the following two statements from the paper are incorrect or at least too strong for what is presented in the paper:\n\n1. Quote from the Abstract (also mentioned in contributions): \"The covariance derived from the Gaussian mixture spectrum corresponds to the known spectral mixture kernel. We are the ﬁrst to discover its rationale and effectiveness for PDE solving.\" In my opinion, this is incorrect. Härkönen et al. (2023) construct a specific kernel for linear PDEs, which as they write explicitly at the end of Section 4.1 recovers the spectral mixture kernel (Wilson et al. 2013) as a special case.\n2. \"By contrast, for PDE solving, it is natural to estimate the solution values on a grid, which opens the possibility of using Kronecker products for efﬁcient computation. To our knowledge, our work is the ﬁrst to realize this beneﬁt and use the Kronecker product structure to efﬁciently solve PDEs.\"I disagree with this being the first instance of Kronecker product structure being used to efficiently solved PDEs. First of all, Kronecker product structure (on regular grids) is being used to solve PDEs (either via preconditioners or via the connection of separation of variables and tensor numerical methods (e.g. Gavrilyuk et al. 2019)). Second, as the authors mention themselves in the manuscript, the computational efficiency of product kernels on regular grids is well-known for GPs (Saatci, 2012). Further, PDE solvers based on GPs which use tensor product kernels (specifically Matern) have also been proposed previously (e.g. Wang et al. 2021).\n\n## Quality of Presentation\nThe paper is generally well-written and easy to understand. The ideas were presented in a coherent way and accessible to readers outside the specific subfield. The plots were largely informative, but I would have preferred larger plots to make them more legible. \n\n## Theoretical Results\nThe theoretical results are mostly standard applications of known results and look correct to me, although I did not check every derivation in detail.\n\n## Experimental Evaluation\nThe experimental evaluation is largely well-done. I have some questions about the baselines that are compared against. Specifically, the vanilla GP baseline with a Matern kernel. Is this using a product kernel? The authors write in the paper that \"[...] it is extremely costly or practically infeasible for the existent GP solvers to incorporate massive collocation points, due to the huge covariance matrix.\". Why should one not be able to leverage the Kronecker product + Toeplitz structure (see Sections 3.1 and 3.2 of Wilson et al. 2015) of a product of Matern kernels for the vanilla GP case as well? If this was not done here, the baseline that is presented here seems rather weak. I can imagine that a scale mixture still outperforms a vanilla Matern kernel, but the performance gap in terms of time and memory should be significantly smaller I would expect.\n\n## Related Work\nThe paper discusses related work from the domain of PINNs well, but misses some of the work on Gaussian-process based PDE solvers. For example, the following papers used GPs to solve PDEs:\n\n- J. Cockayne, C. Oates, T. Sullivan, M. Girolami, *Probabilistic numerical methods for PDE-constrained Bayesian inverse problems*, AIP Conf. Proc. 1853 (2017)\n- Wang, Junyang, et al. *Bayesian numerical methods for nonlinear partial differential equations.* Statistics and Computing 31, 2021 URL: https://arxiv.org/abs/2104.12587\n- Pförtner, Marvin, et al. *Physics-informed Gaussian process regression generalizes linear PDE solvers.* arXiv preprint arXiv:2212.12474 (2022).\n- Chen, Yifan, Houman Owhadi, and Florian Schäfer. *Sparse Cholesky factorization for solving nonlinear PDEs via Gaussian processes.* arXiv preprint arXiv:2304.01294 (2023).\n\n**References**\nPapers cited in the review for reference:\n\n- Härkönen, Marc, Markus Lange-Hegermann, and Bogdan Raita. *Gaussian Process Priors for Systems of Linear Partial Differential Equations with Constant Coefficients.* International Conference on Machine Learning (ICML), 2023, URL: https://arxiv.org/abs/2212.14319\n- Wilson, Andrew, and Ryan Adams. *Gaussian process kernels for pattern discovery and extrapolation.* International Conference on Machine Learning (ICML), 2013, URL: https://arxiv.org/abs/1302.4245\n- Wang, Junyang, et al. *Bayesian numerical methods for nonlinear partial differential equations.* Statistics and Computing 31, 2021 URL: https://arxiv.org/abs/2104.12587\n- Gavrilyuk, Ivan and Khoromskij, Boris N.. *Tensor Numerical Methods: Actual Theory and Recent Applications* Computational Methods in Applied Mathematics, vol. 19, no. 1, 2019, pp. 1-4. URL: https://doi.org/10.1515/cmam-2018-0014\n- Wilson, Andrew Gordon, Christoph Dann, and Hannes Nickisch. “Thoughts on Massively Scalable Gaussian Processes.” arXiv, November 5, 2015. URL: https://doi.org/10.48550/arXiv.1511.01870.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "IBCqE5simr", "reviewer": "Reviewer_SL7y", "summary": "The authors propose a solver for partial differential equations (PDE), specifically designed for high-fidelity and multi-scale PDEs, using Gaussian Process. Current PDE solver methods can be highly unstable and sensitive to hyperparameters. The authors bypass this problem by modeling the solution in the frequency domain and estimate target frequencies from the covariance function of the Gaussian process. The authors also propose an efficient algorithm to scale up the learning algorithm.", "questions": "See weakness section.", "limitations": "", "rating": 6, "confidence": 2, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The authors described their method in detail, analyzed the runtime complexity, and provided multiple sets of experiments.", "weaknesses": "The paper abstract right now is very hard to follow for interested readers, instead of focussing on what the authors did in their methodologies, it should focus mainly on the contributions from a high level. The introduction clarifies the authors’ motivation well, however certain rearrangement of figures can help make it better. When the authors say their method is focussed on high fidelity and multi-scale PDEs, they can exhibit this with a figure (perhaps by moving figure 2,3 up in the first 2 pages).\n\n“While effective, the performance of this method is unstable, and is highly sensitive to the number and\nscales of the Gaussian variances, which are difficult to choose beforehand” - this statement should be associated with an example figure/experiment/prior work.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "qvXLdLfe8t", "reviewer": "Reviewer_FJgw", "summary": "In this paper, the authors consider the problem of solving high-frequency partial differential equations (PDEs) using a Gaussian Process approach. The aim is to propose an alternative to PINNs, which struggle for high-frequency or multiscale PDE using a GP framework. To do so, the authors consider a prior distribution of the solution to the PDE as a mixture of student t distributions to capture high frequencies. Then, they propose a computationally fast GP method, which consists of constructing a sample grid on the domain through a Cartesian product of points, which is then exploited in the algorithm.", "questions": "- Is there any concrete advantage of using the Student t distribution over the Gaussian distribution? The authors state that the Student t distribution is beneficial for high frequency but the performance seems very similar in the experiments of Table 1.\n- The experiments in Section 6 should provide the source terms (at least in Appendix) to allow for reproducibility.\n- Equation (6) should include the domain on which the integration is performed.\n- Minor comment: 2nd sentence of p.2: \"the performance of this method is unstable\" -> \"this method is unstable\"\n- Section 6: Timings should be reported in main text as it's one of the main advantage of the proposed method.", "limitations": "", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "- The paper is well-structured and the problem is clearly introduced and motivated.\n- The authors perform several numerical experiments on synthetic datasets, which show that their method outperform PINNs in terms of resulting accuracy of the solution.", "weaknesses": "- One of the weaknesses of the paper is that the difference between prior GP-based works in the literature is unclear. The authors should clearly state the differences between their work and the existing ones, in particular the papers by Chen et al (JCP, 2021) and Raissi et al (JCP, 2017), cited in this work. I believe that the two differences are: (1) The use of student-t distribution and (2) using a Kronecker product grid to speed-up computations. In the current version, the title and abstract of the paper might be misleading to readers in the sense that this paper is not the first one to consider solving a PDE with GP.\n- I am very surprised by the Spectral method experiment performed by the authors. A quick experiment with Chebfun (https://www.chebfun.org/) allows me to solve the 1D Poisson equation to 10 digits accuracy in a fraction of a second.\n- The authors state at the end of Section 5 that their work is \"the first to [...] use the Kronecker product structure to efficiently solve PDEs\". I believe this sentence is erroneous. It is completely natural when using spectral methods to exploit the Kronecker product structure of the domain (see e.g. Fortunato & Townsend, 2019).\n- The authors offer very little conclusion and the main limitation section from the Appendix should appear in the main text. One of the key limitation of the method is that it is limited to simple geometry, where one can use a Kronecker product structure. However, this severely limits the applicability of the method and questions its advantages over spectral methods.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "l1IZ5JmNI8", "reviewer": "Reviewer_J8Ew", "summary": "The paper discusses the numerical solution of multi-scale partial differential equations (PDEs) with Gaussian processes (GPs).\nTo this end, the submission proposes to use variations of the spectral mixture kernel as a prior model.\nAdditionally, selecting a covariance function that factorises into the product of dimension-wise covariances implies a Kronecker structure in the covariance matrices.\nAs a result, the algorithm requires less-than-cubic complexity to evaluate the loss function, and a user can optimise the PDE solution and (hyper)parameters efficiently. \nExperiments demonstrate the advantage of this method over physics-informed neural networks and those Gaussian-process-based solvers that use more traditional prior distributions.", "questions": "None", "limitations": "", "rating": 6, "confidence": 5, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The topic of the paper (machine-learning-based PDE solvers) is of growing interest in the ICLR community and contributes some valuable concepts to this area. More specifically, the following contributions are novel, to the best of my knowledge:\n\n* Using a spectral mixture kernel for solving multi-scale PDEs\n* Numerically optimising the joint distribution of PDE solution and PDE information can be efficient if a Kronecker-factorised prior is used (something that is not necessarily the case for more traditional MLE or MAP estimation without further assumptions).\n\nThe manuscript is easy to follow for a reader who is roughly familiar with physics-informed neural networks, PDEs, and GPs.\nOverall, I consider this to be a good paper.", "weaknesses": "The weaknesses of this submission are mostly presentational (with two exceptions that concern the experiments).\nI believe the following issues can _and should_ all be corrected before publication.\n\n\n1. The manuscript misses a range of closely related work about solving partial differential equations with Gaussian processes:\n\n    * Howhadi [1] discusses Bayesian numerical homogenisation (solving multi-scale PDEs with GPs), a central piece of related work that this paper should address prominently. \n    * Accessing the derivatives of the Gaussian process as a function of $\\vec{\\mathcal{U}}$ (the trick suggested in Equation (13)) is discussed in Section 2 and Section 3 in the paper by Krämer et al. [2].\n    * The related work section (Section 5) mentions the work on solving PDEs with GPs by Chen et al. (2021) and Long et al. (2022b) but completely misses all work that belongs to the field of probabilistic numerical algorithms (see, for instance, [1-3] below and the references therein).\n\n    I am aware that large parts of the ICLR community may be more familiar with physics-informed neural networks than with Gaussian-process-based differential-equation solvers and that taking this into account for the exposition therefore makes sense for an ICLR paper. However, the works mentioned above are closely related and must be discussed appropriately.\n\n2. I am surprised the submission never mentions traditional approaches to solving multi-scale problems, such as those based on numerical homogenisation (a starting point could be [4]). I know that the paper targets the ICLR community and that large parts of this community may be more interested in machine-learning-based PDE solvers than traditional PDE solvers. Still, existing approaches should be credited by (at least) mentioning their existence. (This relates to point 4. below)\n\n3. I am also surprised that the experiments only discuss the precision of the approximations rather than the work required to achieve this precision. The experiments in Appendix C suggest that the per-iteration runtime of the proposed algorithm is comparable to that of a physics-informed neural network. However, Figure 9 indicates that tens of thousands of iterations are needed. What is each algorithm's overall runtime (training time) to achieve the precisions in Tables 1 and 2?\n\n4. What would we find if we included a non-machine-learning-based solver in the experiments and compared the runtime to the accuracy of all these methods? I would like to see this comparison. I am not saying that the proposed method must outperform traditional solvers, which have been studied and optimised for a long time, but that the context is critical.\n\n\n5. This could be personal taste, but I think the paper might benefit from a slightly more rigorous notation surrounding Equation (14), for instance, by clarifying which joint probability density Equation (14) is the logarithm of. I am imagining something like mentioning $$\n\\mathcal{L} = \\log p(\\mathcal{U}, \\mathcal{H} \\mid \\Theta, \\tau_1, \\tau_2)$$ somewhere at the beginning of Section 4. Contextualising this term against the maximum-a-posteriori and maximum-likelihood loss for the same estimation problem would be even better. What do the authors think?\n\n\nI consider points 1 and 2 essential and 3 and 5 slightly less important to resolve before publication. \nI expect all four of those issues to be relatively straightforward to fix.\nPoint 4 would be optional, as incorporating it is non-trivial. However, I would like to see those results.\n\nI am open to revising my score once these issues have been addressed.\n\n\n\n\n\n**References:**\n\n[1] Howhadi, Bayesian numerical homogenization. Multiscale Modeling \\& Simulation. 2015.\n\n[2] Krämer et al.: Probabilistic numerical method of lines for time-dependent partial differential equations. AISTATS. 2022.\n\n[3] Hennig et al.: Probabilistic numerics. Cambridge University Press. 2022.\n\n[4] Altmann et al.: Numerical homogenisation beyond scale separation. Acta Numerica. 2021.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (poster)