Toward Efficient Kernel-Based Solvers for Nonlinear PDEs
Authors: Zhitong Xu
Da Long
Yiming Xu
Guang Yang
Shandian Zhe
Houman Owhadi
Venue: ICML 2025 poster (2025)
Keywords: PDE Solving
Nonlinear PDEs
Kernel methods
Gaussian processes
Machine Learning based Solvers
Primary area: general_machine_learning->kernel_methods
TL;DR: We developed an efficient kernel-based solvers for nonlinear PDEs.

Abstract: We introduce a novel kernel learning framework toward efficiently solving nonlinear partial differential equations (PDEs). In contrast to the state-of-the-art kernel solver that embeds differential operators within kernels, posing challenges with a large number of collocation points, our approach eliminates these operators from the kernel. We model the solution using a standard kernel interpolation form and differentiate the interpolant to compute the derivatives. Our framework obviates the need for complex Gram matrix construction between solutions and their derivatives, allowing for a straightforward implementation and scalable computation. As an instance, we allocate the collocation points on a grid and adopt a product kernel, which yields a Kronecker product structure in the interpolation. This structure enables us to avoid computing the full Gram matrix, reducing costs and scaling efficiently to a large number of collocation points. We provide a proof of the convergence and rate analysis of our method under appropriate regularity assumptions. In numerical experiments, we demonstrate the advantages of our method in solving several benchmark PDEs.

Toward Efficient Kernel-Based Solvers for Nonlinear PDEs
Zhitong Xu
Da Long",
Yiming Xu
Guang Yang' Shandian Zhe Houman Owhadi
Abstract
designs and intricate numerical techniques, enabling sim- pler, faster, and more convenient implementation and use. Among these solvers, kernel methods or Gaussian processes (GPs) (Williams & Rasmussen (006) hold promise due to their solid mathematical foundations, offering high expres- siveness, robustness, and the ability to quantify and reason under uncertainty. Recently, Chen et al. (2021) introduced a general kernel method to approximate solutions of nonlinear PDEs. They augmented the representation by incorporat- ing differential operators (more generally, linear operators) into the kernels, and jointly estimate the solution values and their derivatives of the PDE on a set of collocation points. The approach has shown promising performance in solving several benchmark nonlinear PDEs, backed by a rigorous error analysis, including both convergence and convergence rates (Chen al. 202T Batlle et al. 2023).
We introduce novel kernel learning framework toward efficiently solving nonlinear partial differential equations (PDEs). In contrast to the state-of-the-art kernel solver that embeds differential operators within kernels, posing challenges with a large number of collocation points, our approach eliminates these operators from the kernel. We model the solution using standard kernel interpolation form and differen- tiate the interpolant to compute the derivatives. Our framework obviates the need for complex Gram matrix construction between solutions and their derivatives, allowing for a straightforward implementation and scalable computation. As an instance, we allocate the collocation points grid and adopt product kernel, which yields a Kronecker product structure in the interpolation. This structure enables us to avoid computing the full Gram matrix, reducing costs and scaling efficiently
Despite its success, the methodology requires manual con- struction of Gram matrix between the solution and its derivatives that show up in the PDEs. It enumerates com- binations between the operators on the kernel to compute different sub-blocks. The process enlarges the size of the Gram matrix (as compared the number of collocation points), and the computation becomes challenging with large number of collocation points, crucial factor for cap- turing complex PDE solutions such as those with potential roughness (or even non-smoothness) and higher frequencies.
large number of collocation points. We provide proof of the convergence and rate analysis of our method under appro- priate regularity assumptions. In numerical experiments, we demonstrate the advantages of our method in solving several benchmark PDEs. Our implementation is released at afttpss
In response, this work proposes an alternative kernel-based framework for solving nonlinear PDEs with several key contributions:
weaithub.comm/aaysliaaa.... Efficient Kernel PDE Solver
Framework: We remove the basis functions associ- ated with the differential evaluation functionals in the approximation and use a standard kernel interpolation for solution modeling. To approximate solution deriva- tives, we directly differentiate the interpolant. By min- imizing the RKHS norm along with a boundary and residual loss, estimate the solution without the need for manually constructing complex Gram matrix. Im- plementation is straightforward and convenient with the aid of modern automatic differential libraries.
Introduction
Solving partial differential equations (PDEs) stands as a central task in scientific and engineering domains. Recently, machine learning (ML)-based solvers have garnered sig- nificant attention. Unlike traditional numerical methods, MI.-based solvers eliminate the need for complex mesh
'Equal contribution 'Uiiversity of Utah University of Ken- tucky "California Institute Technology. Correspondence to: Shandian Zhe <zhe utah.
Computational Method: Our framework allows an immediate use of many existing efficient Gram matrix computation and/or approximation techniques. As an instance, we propose to place collocation points on grid and employ product kernel over each input
Proceedings of the nd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Toward Efficient Kernel- Solvers for Nonlinear PDEs
dimension. This choice induces a Kronecker product structure within both the kernel interpolation and its differentiation, bypassing the need for computing the entire Gram matrix. Such modification results in a substantial reduction in computational costs, enabling efficient processing of tens of thousands or even mil- lions of collocation points. This achieved without any sparse and/or low-rank approximations.
of linear operators, namely,
P())(x) P
xE ?
B(w)(x) = B
(u((u)(x)) E 0
where P(.) and B()) are nonlinear functions, each L3 (1 - - Q) is a linear operator, such as dm U, d2r2z U. a Orgzz u + b u, etc.
Theorem: We provide rigorous analysis of our frame- work. We show the convergence and convergence rate of our method under appropriate PDE stability and reg- ularity assumptions that are similar to the assumptions used in the prior work (Batle et al. 2023). However, our results are not trivial extension of the prior work in that while our framework uses a reduced model space for efficient computation, our convergence re- sults are comparably strong those in the prior work (Chen et al. 2021 Batlle et al. (023) that em- ploys a richer model space. This is achieved through more sophisticated proof. We construct an interpolate of the true solution as an intermediate bridge connect- ing the true solution and our approximation, via which we prove our learning objective is not only feasible, the learned approximation also has a bounded RKHS norm. Next, via using domain decomposition, sam- pling inequality and mean inequality, we are able to bound the L2 norm of the error w.r.t the PDE opera- tors. Combined with the bounded RKHS norm of our approximation, we establish convergence and conver- gence rate. The results theoretically affirm the efficacy of our method in yielding accurate solutions.
To solve the PDE ( Chen et al. (2021) proposed sam- ple a set of collocation points, M ( Ω2 X XMr E aQ], and estimated the solution val- ues and all the relevant linear operators over the solution, (+)(s)())}, evaluated at the collocation points. Specifi- cally, a nested optimization problem was formulated as
minimize uEll s.t. Lg (x)((Xxm) = 2)
||ulle
minimize
1 1 m Ma for j= =1..QQ. Qn.
Mn + ++= S m S M for j = Qn + 1...Q
s.t. P(z
2, Qn) f(x,) 0 (1 s S Mn),
B(2 90+1
29) - 9(xx)= = (Mo + +1 9 S M).
where 2 is the Reproducing Kernel Hilbert Space (RKHS) associated with kernel k((;-). Let us denote all (=// z. With the RKHS U, the similarity (or covariance) between any zí1 and z/2 in Z is
c(2 c : )
Ljx (K)(Xm)) Xmz);
where Lj, is applied along the first argument of k(-,-) and Lie is applied along the second argument, and then we evaluate at Xm, and Xmy- For example, suppose L = O, u and Lạ(u) = 0zgu, then L1 o [2(x)) x') 2 Px(x, xx)/0r,0rr where Ij and x are the first and sec- ond elements of the inputs and x', respectively.
Experiments: Evaluation on Burgers', nonlinear ellip- tical. Eikonal, and Allen-Cahn equations have demon- strated our method's practical efficacy. For less chal- lenging scenarios where a small number of collocation points, e.g., 1000, sufficient, our method achieves comparable or sometimes smaller errors than the ex- isting methods. In more challenging scenarios, such as Burgers' equation with viscosity of 0.001, our method seamlessly scales to tens of thousands of col- location points, yielding low errors on the order of 10-3 to 10-6, underscoring its robustness and accu- racy across a range of problem complexities.
The true solution u" is assumed to reside in U4, and || 1lr is the RKHS norm associated with U. From the representation theorem (Owhadi & Scovel, 2019). it is straightforward to show that the optimum of the aforementioned nested optimizatin problem, denoted by UMs takes the form
uar (x) c(vM((x), C),, 2)
(2)
where ((UM(()) 2) is the similarity between (x)(x) and each element in z, namely, for every 2/ E z,
2
Background
(x)) 2/h) = Ly (x),xxxmm
where L, is applied along the second argument of k((,)) )), and xm is the collocation point corresponding 2/1 Here
Consider a PDE of the general form,
C(z, z) is the Gram matrix of z, constructed from Q x Q sub-blocks,
P(u) (x) (x ()), B(u)= 9(x (x 82), (1)
Cn
CiQ
where x = (ZI,..., x.) T and P and B are nonlinear dif- ferential operators in the interior Ω and boundary 02, re- spectively. We assume P and B are composed by a series
C=
:
(3)
Cqu
Coo
2

---

Toward Efficient Kernel- Solvers for Nonlinear PDEs
where each Cty is the similarity matrix associated with pair of linear operators,
We then apply each linear operator Ly in the PDE over "(x; 17) to approximate [yyx"). Following this way, we do not need to explicitly estimate the values of Lj(u*) at the collocation points as in (Chen et al. ] et al. 2022) (namely z in z). Correspondingly, the Gram matrix K MM is substantially smaller (of size M x M) and it is more convenient to compute - there no need to enumerate pairs of linear operators and apply them to the kernel function to compute different sub-blocks.
City - - h1) - Li (x), hy),
(4)
where S i S Q. L, and Li are applied to the first and second arguments of R, respectively, h, (em)n Im and hj {ah Ju- Therefore, the size of the Gram matrix C is (QuMa + (Q- Qa)(M - M))) x (QnMa + (Q - Qa) - Ma)).
The learning of our model is carried out by addressing the following constrained optimization problem:
This approach can be explained from a probabilistic per- spective (Long et al. 2022). That is, we assign a GP prior over u, and given a sufficiently smooth kernel K, all the linear operators over u, namely, (u) also follow a GP prior, and their projection the collocation points M,, namely, 2, follow multi-variate Gaussian prior distribu- tion, p(z) ~ [zz C). Softening the outer constraints in the aforementioned nested optimization by maximizing a likelihood, this method essentially seeks for an MAP esti- mation of z, and the prediction (2) is the posterior mean conditioned on z.
minimize 1/22 HEll
s.t. Mn Ema (P(u)(xm) - fxm))
+ M-M Em=Ma+1 (B(u)(xm) - g(xm)) ))² Se (6)
ΣM
u takes the kernel interpolation form
where E 2 0 is a given relaxation parameter. Note that since our formulation (5) uses reduced model space (there are no basis functions associated with operators {L,) }) as opposed enable feasibility of
to (Chen et al. 2021), we introduce the optimization and to establish the convergence; see our convergence analysis in Appendix Section A In practice, addressing (6) directly can be cumbersome. We may opt to
3 Our Framework
Despite the success of (Chen et al. 2021). requires compu- tation of a Gram (covariance) matrix (see (3) with dimen- sions typically exceeding the number of collocation points. This can exacerbate computational challenges, particularly when addressing complex PDEs that demand a consider- able number of collocation points (Cho et al. 2024| Florido let al. 2024). Moreover, the construction of the Gram matrix relies on the particular set of linear operators present in the PDE, rendering it cumbersome for implementation and the adoption of efficient approximations, if needed.
optimize an unconstrained objective with soft regularization instead,
minimize L(u(a;)): a, = l|ulie
1
+ αK
Mn Σ
Ma
(P(u)(xm) - [(xm)) f/2
(7)
m-
M
+ 3
m= (B(u)(xm)) - -m0x)=))))(())))))) 6/2
M - Ma
We therefore propose alternative kernel learning frame- work for nonlinear PDE solving, which simplifies the Gram matrix construction and computation. Specifically, we are inspired by the standard kerne//GP regression. Suppose the solution values at the collocation points are known, denoted ))² T. The optimal solution estimate within the framework of standard kernel regression takes the interpolation form: t(x) = x(x, .M) M U'M where K M x(M, M) denotes the kernel matrix com- puted on the collocation points (of size M M). This form is derived by minimizing the RKHS norm while aligning u² at M. In GP regression, t(x) serves as the mean function of the posterior process.
where a, 3 >> 0 are the regularization strengths, and can be simply set to zero
Efficient Computation. In scenarios where PDEs are com- plex and challenging, capturing the solution details might necessitate using a vast array of collocation points. Since our framework uses the standard Gram/kernel matrix con- struct the solution estimate (as illustrated in range of existing kernel approximation and computation methods Guinonero-Candele & Rasmussen 2005 Rahimi & Recht 2007: Farahaat et al. 2011: Lindgren et al. 2011)
as u'u- = (u* (xi(x),...,, 44 "(xx))
wide
can be readily employed to accelerate computation involving K WM. This facilitates the reduction of computational costs and enables scalability to accommodate massive collocation points.
For general PDE solving, as depicted in I.. one often lacks knowledge of the solution values at arbitrary collocation points. Therefore, we regard them as unknown, free vari- ables denoted by n7. We model the solution estimate
As an instance, we propose to induce a Kronecker prod- uet structure to accelcrate computation and scale to large number of collocation points. Specifically, we place the col- location points on a grid, namely, M = s x...x.x sd, where each s includes a collection of locations at input dimension
u(x: n) x(x, MMK,MM M.
3

---

Toward Efficient Kernel- Solvers for Nonlinear PDEs
k, ie., s = (of, :........ mm I These locations can be regular-spaced or randomly sampled. Accordingly, M is an d-dimensional array of size my we employ a product kernel that is decomposed as along the input dimensions, x(x,x' x') =I 2z)) where each Ki is a kernel function of two scalar variables at input dimension j. An example is the widely-used Square Expo- nential (SE) kernel. The product kernel corresponds to a tensor product structure in the latent feature space, offering rich representational power to approximate PDE solutions; see the theoretical justification in (Wang al. 2021).
assumption about PDE stability and the regularity of the domain and boundary (Batlle et al. 2023 Assumption 3.7)
x ... X m.- Next,
with a slight modification. Assumption 4.1. The following conditions hold:
(CI) (Regularity of the domain and its boundary) Q C IRd with d > is a compact set and 0n is a smooth connected Riemannian manifold of dimension d endowed with a geodesic distance paw.
(C2) (Stability of the PDE) 3k, t E N with A > d/2 and t >> (d - 1)/2, and Es, E R such that for any r > 0, it holds that Vun, u2 E (2)).
With the product kernel, the kernel matrix on the collocation points M becomes a Kronecker product, KMM = K, @ @ Ka. where each K = (s², s') is local kernel matrix for dimension j (1 S j I d), of size my x my. We then leverage the Kronecker product properties to efficiently compute the solution estimate (5) as
I - w2 weller S ((I((u)) - -00--PP-ee P(22))l + [B((n) - ((22) Ho((a))
(10)
and Vun, u2 B,
[m (Ii,s') s') @ ... @ Ka(Ea, s*)
KK 2 ... 2 Kil'' n
[|P(un) - P(u2))l( + [BB(un) B(u2)]l
= [K1( KT S... & Ka( sd) K
(8)
= =A x1 (11, " X2 xa [Nd 8*) K
I C [|] - 1/2 (+))
(11)
where A is the tensor view of 77, namely reshaping 77 as my x... X ma array, and Xx is the mode-k tensor-matrix multiplication (Kolda, 2006). In this way, we avoid ex- plicitly computing the full kernel matrix KMr and its in- verse. We only need to invert each local kernel matrix Kj, and hence the cost is substantially reduced. For example, considering 100 x 100 x 100 grid, the full kernel ma- trix is 10 x 108, rending it computationally prohibitive and impracticable for most hardware. By using (8. we only need to invert three 100 x 100 local kernel matrices, which is cheap and fast. Furthermore, since the kernel is decomposed across individual dimensions, taking deriva- tives over the solution estimate u will maintain the structure, e.g., Oeina u(x: 17) = A X1 K2 s' )KT X2... Xd Ka st) '')
where C= CCr) > 0 is a constant independent of d1 and u2. B(r) is an open ball with radius r, H' = W is a Sobolev space where each element and its weak derivatives up to the order of j have a finite I norm.
. (C3) The RKHS U is continuously embedded in HS+ where T > 0.
Lemma 4.2. Let u² E l denote the unique strong solution
of []. Suppose Assumption 4.7 is satisfied, and set of
collocation points M C S given, where Ma C M denotes the collocation points in the interior Q and Man C M the collocation points on the boundary An. Assume the Voronoi diagram based on the collocation points has a uniformly bounded aspect ratio across all the cells 2 Define the fill-distances
We then leverage the structure g efficiently minimize the objective 66) or (7). The computation of each P (x))(x)) and (xu)(x... is straightforward application of the operators P and B to (8) and then evaluate them at the collocation points. This can be done by automatic differential libraries, such as JAX (Frostig et al. [0188.. The RKHS norm in (6) and 7 can be efficiently computed by
hn := sup inf |x - -x'I, xEn x E
han := sup inf Xcon x'E
(x,,,xx'),
(12)
where |- is the Euclidean distance, and pan is geodesic distance defined on On. Set h = max(by, hag). There is
|||| - = 17 T " =n T K2 &... & KJ) n = = T vec x K X2 X )
'Note that this assumption, along with its minor variants, is considered mild and widely used in convergence analysis. For many examples of nonlinear PDEs that satisfy this assumption, see (Batlle et al. (023)..
We can apply any gradient-based optimization algorithm.
"Specifically, for each cell T. that includes X the aspect ratio is defined as the ratio between the smallest radius of a ball centered at X... containing the closure of T., and the maximum radius of the ball centered at contained in Tm. The aspect ratio is assumed to be uniformly bounded across and the total number of collocation points.
4 Convergence Analysis
We now show the convergence of our method. We inherit the road-map of (Batile et al. 2023) and maintain the same
4

---

Toward Efficient Kernel-Basch Solvers for Nonlinear PDEs
always a minimizer of O with the set of collocation points M and c Coh where Co > O is a sufficiently large constant independent of h. Let ul denote such minimizer. When h is sufficiently small, at least h S &CMM- à where C1 > O is a constant, then
relatively inconvenient for implementation. In our work, the solution is approximated by standard kernel interpola- tion, and the Gram matrix is therefore just the kernel matrix over the collocation points. The size of the Gram matrix is smaller. More important, the existent sparse approximation methods can be readily applied to our model, without the need for complex adjustments or novel development.
|/u| - u² (()) S Chelluu lec:
(13)
where p min(k, 4,,, and C > 0 is independent of ut
The computational efficiency of Kronecker product struc- tures has been recognized in various works (Saateci 2012
and h.
Proposition 4.3. Given the set of collocation points M and c c= =Chh where Co > O is a sufficiently large constant, there exists a ,ar >> O such that the minimizer of (7) with aM and 3 BM is also the minimizer of (6) That means, with proper choices of the regularization strengths, the minimizer of 7 enjoys the same convergence result as in 13
Wilson & Nickisch/ 2015 Izmailov et al. 2018 Zhe et al.
2019). Wilson et al. (2015) highlighted that utilizing a regu- lar (evenly-spaced) grid results in Toeplitz-structun
kernel matrices, facilitating On log n.) computation. However, in typical machine learning applications, data is not observed on a grid, limiting the utility of the Kronecker product. In contrast, for PDE solving, estimating solution values grid is natural, making Kronecker products combined with kernels a promising avenue for efficient computation. The recent work (Fang et al. 2024) uses similar computational method to solve high-freet frequency and multi-scale PDEs. The major contribution is to introduce spectral mixture kernel in each dimension to capture the dominant frequencies in the kernel space. This work can be viewed as instance of our proposed framework. We in addition give theoretical analysis about the convergence of our framework. Exten- sive discussions on Bayesian learning and PDE problems are given in (Owhadi 2015).
We can see the convergence results of our framework are as comparably strong as the results for the method of Chen eet al. (202T); see (Batile et al. 2023] Theorem 3.8), though the latter employs a richer model space. We leave the proof in Section A and of the Appendix.
5 Related Work
The prior works of Graepel (2003); Raissi et al. (0117) propose Gaussian Process (GP) models for solving linear PDEs in the presence of noisy measurements of source terms. Chen et al. (202T) introduced a kernel method capa- ble of solving both linear and nonlinear PDEs. The solution approximation is constructed by both kernels and kernel derivatives (more generally, the linear operators of the PDE over the kernels). Hence, the differentiation operators need to be embedded into the kernels to construct the Gram ma- trix whose dimension is typically greater than the number of collocation points. Long et al. (2022) proposed a general framework for incorporating differential equations into GP learning. which is conceptually equivalent to (Chen et al. 2021). However, Long et al. (2022) formulated the method within a probabilistic modeling framework and developed variational posterior inference algorithm. In (Batlle al.
Our work is related to radial basis function (RBF) meth- ods (Hardy 1971] Kansa 1990; Powell 1992 Buhmann
2000 Tolstykk & Shirobokov 2003), which approximate the PDE solution as a linear combination of RBF basis functions, such as Multiquadrics and Gaussian: u(x) =
Σ, (x), where (x) _ <[1x - x,ll). {x,) denote collocation points, and || || is usually chosen as the L2 norm. To estimate the coefficients a = {a,) RBF methods typically convert the PDE into a linear system Aa b. There are also methods that estimate the coefficients via residual minimization, e.g., (Böhmer & Schaback 2020). A key distinction between our method and RBF methods lies in the solution modeling. Our method uses a kernel interpo- lation form (5) to approximate the solution. We can draw an implicit connection via setting a 1= n, where we directly estimate n- - the solution values at the collocation points. One might question why not optimize instead. Empirically, we found that, however, doing significantly degrades performance - often by one order of magnitude in our test benchmarks. This degradation might stem from the fact that each coefficient in a has a strong global influ- ence on the entire solution approximation, making the opti- mization process more challenging, especially for nonlinear PDEs. Additional details and ablation studies supporting this observation are provided in Appendix Section C.5 It is worth noting that while local RBF methods (Fornberg)
2023). systematic theoretical framework is established to analyze the convergence of the method Chen et al. (2021 To alleviate the computational challenge for massive collo- cation points, Chen et al. (2023) adapted the sparse inverse Cholesky factorization [Schafer et al.) 202T) to approxi-
MM.
mate the Gram matrix of (Chen et al., 202 An alternative approach was proposed by Meng & Yang (0023), which ad- justs the Nyström method (Jin et al. 2013) to obtain a sparse approximation of the Gram matrix. Despite the success of these methods, the construction of the sparse approximation needs to carefully handle different sub-blocks in the Gram matrix, where each block corresponds to a pair of lin- ear operators over the kernels, and hence it is complex and
5

---

Toward Efficient Kernel- Solvers for Nonlinear PDEs
& Flyer 2015 Bayona et al. 2010) also approximate the solution and its derivatives via interpolation of nodal val- ues, they primarily rely on neighborhood points through stencils and require repeatedly computing the interpolation weights across different neighborhoods. In contrast, our method employs a single kernel interpolation, offering both convenience and computational efficiency.
method as SKS (Simple Kernel-based Solver) We com- pared with (Chen et al,, 2021) that uses kernel and kernel derivatives (more generally, linear operators) to approximate the solution, which we denote as DAKS (Derivative-Augmented Kernel-based Solver). We used the original implementation from the author In addition, we compared with physics-informed neural net- work (PINN) (Raissi et al.) 2019). a mainstream machine learning based PDE solver. The PINN is implemented with PyTorch (Paszke et al.) 2019). For SKS, we minimize (7) (with = 0), and used ADAM optimization with learning rate 10-3 The maximum number of epochs was set to IM. We stopped the optimization if the loss stopped improving for 1K updates. For DAKS, we used the relaxed Gauss- Newton optimization developed in the original paper. The PINN was first trained by 10K ADAM epochs with learn- ing rate 10 and then by L-BFGS with learning rate 10-1 with a maximum of 50K iterations. The tolerance level for L-BFGS was set to 10 10-9. To identify the architecture for the PINN, we varied the number of layers from (2, 3, 5, 8, 10), and the width of each layer from {10, 20, 30, 100). We used tanh as the activation function. For DAKS and SKS, we used Square Exponential (SE) kernel with dif- ferent length-scales across the input dimensions. We se- lected the nugget term from (5B-5, 1E-5, 5E-6, 1E-6,... 1E-13). However, for solving the nonlinear elliptic PDE with DAKS, we used its default approach that assigns an adaptive nugget for the two sub-blocks in the Gram matrix. This gives the best performance for DAKS. The length- scales were selected from grid search, from [0.1, 0.2] for the nonlinear elliptic and Eikonal PDEs, [0.05, 0.01 for Allen-Cahn equation, and [0.003, 0.05] x [0.02, 0.3] for Burgers' equation. For SKS, we selected and 3 im 77) from the range (10-2, 10010011...... 1010 1012, 10t4 10t5, 1020). jointly with other hyperparameters, including the kernel length scales and nugget terms. To efficiently tune the hyperparameters, we first performed random search to identify a promising group of hyperparameters. We then fixed all other hyperparameters and conducted a grid search over a and B. Finally, we fixed a and B and performed grid search over the remaining hyperparameters. We re- ported the best solution error of each method throughout the running. In Section C.4 of Appendix, we further examined the sensitivity of our method to the kernel parameters.
Numerical Experiments
To evaluate our method, we considered four commonly- used benchmark PDE families in the literature of machine learning based solvers (Raissi et al. 2019; Chen et al. 2021
The Burgers' Equation. We first tested with viscous Burgers' equation,
Ut + utl yuys = O, V(z, ( E (-1,1) x (0,1), u(2; 0) = sin(TT)) u(-1,t) - (11 = 0. (14)
The solution is computed from the Cole Cole-Hopf transforma- tion with numerical quadrature (Chen et al. (202T). We considered two cases: V/ 0.02, and v= 0.001.
Nonlinear elliptic PDE. We next tested with the instance of nonlinear elliptic PDE used in (Chen et al. 2021).
Au(x) + u² (x) f(x). Vx Q u(x) = 0, Vx E an,
(15)
where Ω [0, 1j. the solution is crafted as u(x) = sin(Tæ1) sin( (Tx2) + 4 sin(4 (4xz1) sin(4 and fxx is
correspondingly computed via the equation.
Eikonal PDE. Third, we tested with regularized Eikonal equation as used in (Chen et al. 202T,
[u(x))= f(x)² + eAu(x), VX Ω, m(x) 0, Vx E an,
(16)
where Ω [0, 1J2. f(x) = 1, and 0.1. The solution is computed from highlvere finite difference solver as provided by (Chen et al. 2021).
Allen-Cahn Equation. Fourth, we considered a 2D sta- tionary Allen-Cahn equation with a source function and Dirichlet boundary conditions.
6.1 Solution Accuracy
Ugz + Byy + T - -u) =)= f(x
(17)
Simpler Cases. We first tested all the methods on less chal- lenging benchmarks, for which a small number of colloca- tion points is sufficient. Specifically, we tested with Burgers' equation with viscosity 0.02, the nonlinear elliptic PDE, and Eikonal PDE. These are the same test cases employed
10, We crafted the solution in the form of u = sin(2rar)) cos(2 + sin 1in(2))1) cos(22))))) and f is computed through the equa- tion. We tested with a = 15 and a = 20.
where y = 1, m = 3. and (x.
Method and Settings. We implemented our method
httpp:://itthbl wttphhu...ooyyy/ com/ycmm//ooooo fanc96/ f No -PPpDDaa-PPapa GP solver
with JAX (Frostig et al. 2018). We denote our
6

---

Toward Efficient Kernel- Solvers for Nonlinear PDEs
in CChen et al.
20 Following (Chen et al.
2021), we
Table I: L2 error of solving less challenging PDEs, with a small
varied the number of collocation points from (600, 1200,
number of collocation points. Inside the parenthesis of each top row indicates the grid used by SKS, which takes approximately the same number of collocation used by DAKS. Note that the Gram matrix of DAKS is larger than SKS.
2400, 4800) for Burgers' equation, and (300, 600, 1200,
2400) for nonlinear elliptic PDE and Eikonal PDE. In (Chen et al. ) the collocation points are randomly sampled, and hence we used the average L of DAKS from ten runs on different sets of randomly sampled collocation points for comparison. For SKS, we used a regularly-spacee square grid, for which the total number of grid points is close to that used for DAKS. Note that the size of the Gram matrix of DAKS is larger than SKS. We also examined running DAKS with the same set of grid points as used by SKS, but the performance is worse than randomly selected collocation points. See Appendix Section Ceetion/C.6fio for details. We ran PINN on the same set of grid points used by SKS. The L error is reported in Table [ It can be seen that in most cases, SKS achieves smaller solution error than DAKS, with the excep- tion for Burger-ant and nonlinear elliptic PDE using 4800 and 2400 collocation points, respectively. This might be because DAKS needs to explicitly estimate more variables, includ- ing all kinds of linear operators (e.g., derivatives) over the solution at the collocation points, which increases the opti- mization workload. In addition, the influence of the nugget term might vary across different heterogeneous blocks in the Gram matrix, which complicates the optimization. In most cases, the PINN shows worse performance than SKS, except on Burgers' equation with 600, 1200 and 2400 collocation points. The results have shown that in regimes where the computation of the Gram matrix is not a bottleneck, our method SKS, though adopting a simpler model design, can still achieve comparable or even better solution accuracy.
(a) The Burgers' equation
14 with viscosity v 0.02
600
1200
2400
4800
Method
(25 x 25)
(35 x 35)
(49 x 49)
(70 x 70)
DAKS
1.75E-02
7.90F-03
8.65E-04
9.76E-05
PINN
2.68E-03
6.72E-04
3.60E-04
3.73E-04
SKS
1.44E-02
5,40E-03
7.83E-04
3.21E-04
(b) Nonlinear elliptic PDE [15
300
600
1200
2400
Method
(18 18)
(25 x 25)
(35 x 35)
(49 x 49)
DAKS
1.15F-01
1.15E-04
8.65F-04
1.68E-07
PINN SKS
3.39F-01
1.93E-02
1.28F-03
3.20F-04
1.26E-02
6.93E-05
6.80E-06
1.83E-06
(c) Eikonal PDE [16
Method
300
600
1200
2400
DAKS
L.01-01
1.64E-02
2.27/:-04
1.78E:05
PINN SKS
2.95F-02
1.26F-02
4.53E-03
3.50FF03
6,23E-04
2.68F-04
1.91E-04
2.51F-05
at 10 error level or even bigger. Together this indicates that much larger number of collocation points is needed.
We then ran SKS and PINN with greatly increased colloca- tion points, i.e., dense grids, varying from 6400 to 120K. In such scenarios, running DAKS becomes extremely costly or even infeasible We therefore only report the results of SKS and PINN, as shown in Table One can see that the solution error of SKS is substantially reduced, achieving 10-3 for Burgers' and 10-4 to 10-6 for Allen-Cahn. It is worth noting that PINN using the same set of collocation points also arrives at the 10-3 3 level 1,2 error for Burgers' but the error on Allen-Cahn is still very large, with nearly no improvement upon using much fewer collocation points. This might be due to that the relatively high frequencies in the solution (see (17) are difficult be captured by neu- ral networks, due to their known "spectral bias" (Rahaman) et al. 2019). Thanks to our model design ( can induce Kronecker product structure in the Gram matrix to scale to massive collocation points, without the need for design- ing complex approximations. Note that if a product kernel is used for DAKS, each block Ct in 3 will form a Kro- necker product as well. However, the entire Gram matrix will not exhibit this structure. Consequently, the efficient computation method of SKS cannot be applied to DAKS. To directly examine the influence of the fill distance h (i.e., mesh norm), Appendix Fig. la and Fig. Ib show how the
Next, we tested with more challenging cases, for which massive collocation points are necessary. These cases include Burgers' equation with v = 0.001, and Allen-Cahn equation with a = 15 and a = 20. For Burgers' equation, we found empirically that the spatial resolution is more important than the time resolution, so we set the ratio between the spatial and time resolutions to 3:1. For Allen- Cahn, we still used square-shaped grid. We provide more detailed ablation study in Appendix Section C.2 To verify the necessity of using massive collocation points, we first ran all the methods with the same number of collocation points as adopted in the simpler PDEs, namely, a few hundreds and/or thousands. As we can see from Table 2 the solution errors of all the methods are large, typically around the level of 10-1, indicating failures. Note that, however, when the number of collocation points increases to 4800, SKS can achieve an L error at level 10-4 for solving the 2D Allen-Cahn equation, while the other methods still struggle
Difficult Cases.
'However, if we use a 96 x 50 grid (still including 4800 collo- cation points), SKS gives L error 7.54E-05, which can surpass DAKS. See the ablation study in Appendix Section C.2 for more details.
'oor 6400 collocation points, the size of the Gram matrix of DAKS is 19200 x 19200 for Burgers' and 2D Allen-Cahn, hecause there are three linear operators in each equation.
7

---

Toward Efficient Kernel- Solvers for Nonlinear PDEs
Table 2: L error of solving more challenging PDEs with small number of collocation points. The grids used in SKS are the same in Table
Table 3: L error of solving more challenging PDEs with large number of collocation points.
(a) The Burgers' equation
with viscosity v = 0.001.
(a) The Burgers' equation
14) with viscosity v = 0.001
43200
67500
97200
120000
Method
Method
600
1200
2400
4800
(360x 1 20))
(450x 150)
(540x 180)
(600x
DAKS
3.87E-01
3.12B-01
3.60E-01
2.37E-01
PINN
4.05E-03
6.01E-03
3,94E-03
4.13E-03
PINN
1.27E-01
2.59E-01
3.18E-01
2.65E-01
SKS
3.90E-03
3,50E-03
2.60E-03
2.28E-03
SKS
1,34E-01
1.11E-01
8.04E-02
1.89E-02
(b) The 2D Allen-Cahn equation 17
with
15.0
(b) The 2D Allen-Cahn equation (17 with
15
6400
8100
22500
40000
Method
Method
600
1200
2400
4800
(80x80)
(90x90)
(150x 150)
x200x200)
DAKS
6.84E-01
6.62E-01
6.28E-01
5.74E-01
PINN
5,03E0
5,30E0
4.21EO
5,86E0
PINN
4.02E0
6,20E0
4,26E0
5.39E0
SKS
8.27E-05
3.41E-05
4.34E-06
4.44E-06
SKS
6.80E-01
2.1E-01
5.15E-03
9.20E-05
(c) The 2D Allen- Cahn equation
17 with a
20.0
(c) The 2D Allen-Cahn equation
177 with a
20
Method
6400
8100
22500
40000
Method
600
1200
2400
4800
PINN
4,18E0
4,45ED
5,86E0
5.93E0
DAKS
6.81E-01
6.57E-01
6.19E-01
5.64E-01
SKS
3,98E-04
1.82E-04
4.00E-05
2.98E-05
PINN
4,98E0
5,78E0
5,87E0
3.04E0
SKS
7.07E-01
6.91E-01
1881E-01
9,83E-04
Table 4: I.² Error of finite difference solver and SKS according to the ground -ruth solution.
(a) Nonlinear Elliptic PDE 15
error of SKS varies with /2- hー determined by the collocation points - when solving Burgers' equation with v = 0.02 and v = 0.001, respectively.
Method
18 x 18
25 x 25
35 x 35
49 x 49
FD
3.36E-02
1.78E-02
9.25E-03
4.78E-03
SKS
1.26E-02
6.93E-05
6.80F-06
1.83E-06
Point-wise Error. For fine-graines comparison, we show- case the point-wise error of each method in solving Burgers' (u 0.001) and 2D Allen-Cahn equations. The results and discussion are given by Appendix Section
(b) The 2D Allen-Cahn equation
17) with
15.
Method
80 x 80
90 90
150 x 150
200 200
FD
8,57E-02
6.68E-02
2.33E-02
I.30E-02
SKS
8.27E-05
3,41E-05
4.34E-06
4.44E-06
Ablation Study on Grid Shape. We further examined the influence of the grid shape on the solution accuracy. We compared different choices on Burgers' equation with 0.02 and 0.001. We leave the details in Appendix
(c) The 2D Allen-Cahn equation (17 with
20.
Method
80 x 80
90 x 90
150 x 150
200 x 200
FD
1.62E-01
1.24E-01 1.82-04
4.22E-02
2.34E-02
Section
SKS
3.98E-04
4.00E-05
2.98E-05
Comparison with Conventional Numerical Methods. In addition to comparing SKS with ML-based solvers, we also compared it to a finite difference solver - a widely used conventional numerical approach. We discretized the PDE using numerical differences, specifically employing cen- tered second-order numerical difference approximate the derivatives. The equation was then solved using a Newton- Krylov solver, which computes the inverse of the Jacobian through an iterative Krylov method. We tested on solving the nonlinear elliptic PDE (5)) and the Allen-cahn equa- tion (17). since the ground-truth solutions of these PDEs are known and we can conduct a fair comparison. The re- sults are presented in Table 4 Note that the L2 errors of PINN and SKS have already been reported in Table Tb] (3b] and 3c Our method (SKS) consistently outperforms finite difference. In most cases, the error of SKS is several orders of magnitudes smaller. It implies that using the same grid, SKS is much more efficient in approximating the solution. In addition, with the growth of the grid size, the relative
improvement of our method is often more significant, in particular when solving the nonlinear elliptic PDE with the grid 18 18 increasing to 25 25, and Allen- (a 15) with the grid 90 90 increasing to 150 x 150. Since SKS is efficient in handling a large number of collocation points (i.e., dense grid), it shows the potential of SKS.
Irregular-shaped Domains. While our efficient computa- tion is performed on grids, our method can be applied to irregular-shaped domains by introducing (minimal) virtual grid that encompasses such domains. This allows SKS to be used without any modifications. To validate the effective- ness of this strategy, conducted additional tests, solving the nonlinear elliptic PDE (15) on a circular domain and the Allen-Cahn equation
(17) on a triangular domain. In both cases, our method achieved reasonably good accuracy. Detailed results and discussions are provided in Appendix Section
8

---

Toward Efficient Kernel- Solvers for Nonlinear PDEs
Running Time. Finally, we examined the wall-clock run- time of SKS. We analyzed the runtime of each method when solving Burgers' equation Cahn equation (a = 15) with varying numbers of colloca- tion points. We ran the experiment on a Linux workstation equipped with an Intel(R) Xeo Platinum 8360H Pro- cessor with 24GB memory. The results are presented in SKS is several orders of magnitude faster than both DAKS and PINN per iteration. However, since DAKS employs the Gauss-Newton method, it converges much faster than the ADAM optimizer used by SKS and PINN. Nonetheless, the overall runtime of SKS is still less than 25% of that of DAKS when solving Burger's ' equation, and is close to DAKS when solving Allen-Cahn equation. Additionally, SKS can handle a much larger number of col- location points than DAKS. Overall, the runtime of SKS is significantly less than that of PINN.
Fault Network using a Computational Graph Discovery Approach), and from the Dol Vannevar Bush Faculty Fel- lowship Program. SZ acknolwedges support from NSF CAREER Award IIS US-2204299,, and NSF OAC-231
0.001 and the Allen-
SZ also thanks Bamdad Hosseini for insightful discussions.
Impact Statement
Appendix Table 9
This paper presents work whose goal is to advance Machine Learning for PDE solving. There are many potential societal consequences of our work, none of which we feel must be specifically highlighted here.
References
Adams, R. A. and Fournier, J. Sobolev spaces. Elsevier, 2003.
In addition, we examined the runtime of our model using the naive full matrix computation (i.e., without exploiting the Kronecker product structure). As shown, the per-iteration runtime with naive matrix operations consistently around 100x slower compared to our method leveraging Kronecker product properties. Furthermore, when the number of collo- cation points increases to 22,500 for the Burgers' equation and 43,200 for the Allen-Cahn equation, the naive approach exceeds available memory, resulting in out-of-memory er- rors- - rendering it infcasible to run.
Batlle, P., Chen, Y., Hosseini, B., Owhadi, H., and Stuart, A. M. Error analysis of kernel/GP methods for nonlinear and parametric PDEs. arXiv preprint arXiv:2305.049: 0.4622, 2023.
Bayona, V., Moscoso, M., Carretero, M., and Kindelan, M. RBF-FD formulas and convergence properties. Journal of Computational Physics, 229(22):8281-88))) -295, 2010.
Böhmer, K. and Schaback, R. A nonlinear discretization theory for meshfree collocation methods applied to quasi- linear elliptic equations. ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik, (00((10)e20000)))) 2020.
Conclusion
We have proposed new kernel method for nonlinear PDE solving. We use a standard kernel interpolation to model the solution estimate, which allows more convenient imple- mentation and efficient computation. Our method can easily scale to massive collocation points, which can be important for solving challenging PDEs. The performance of benchmarks is encouraging. In the future, we plan to develop more efficient optimization, e.g., Gaussian-Newton, to further accelerate our method.
Buhmann, M. D. Radial basis functions. Acta numerica, 9: 1-38, 2000.
Chen, Y., Hosseini, B., Owhadi, H., and Stuart, A. M. Solv- ing and learning nonlinear PDEs with Gaussian processes. arXiv preprint arXiv:2 103. 12959, 2021.
series
Chen, Y., Owhadi, H., and Schäfer, F. Sparse Cholesky factorization for solving nonlinear PDEs via Gaussian processes. arXiv preprint
Acknowledgments
2023.
Cho, J., Nam, S., Yang, H., Yun, S.B. Hong, Y., and Park, E. Separable physics-informed neural networks. Advances in Neural Information Processing Systems, 36, 2024.
HO and SZ acknowledge support from the Air Force Office of Scientific Research under MURI award number FA9550- 20-1-0358 (Machine Learning and Physics-Basse -ased Modeling and Simulation). HO acknowledges support from the Air Force Office of Scientific Research under MURI award num- ber FOA-AFRL-AFFSSR-222-- tal Twins), the Department of Energy under award number DE-SC0023163 (SEA-CROGS: Scalable, Efficient, and Ac- celerated Causal Reasoning Operators, Graphs and Spikes for Earth and Embedded Systems), the National Science Foundations under award number 2425909 (Discovering the Law of Stress Transfer and Earthquake Dynamics in
Fang, S., Cooley, M., Long, D., Li, S., Kirby, R., and Zhe, S. Solving high frequency and multi-scale pdes with Gaus- sian processes. In International Conference on Learning Representation, 2024.
(Mathematics of Digi-
Farahat, A., Ghodsi, A., and Kamel, M. A novel greedy algorithm for Nyström approximation. In Proceedings of the Fourteenth International Conference on Artificial

---

Toward Efficient Kernel- Solvers for Nonlinear PDEs
Intelligence and Statistics, pp. 269-277. JMLR Work- shop and Conference Proceedings, 2011.
Meng. R. and Yang. X. Sparse Gaussian processes for solv- ing nonlinear PDEs. Journal of Computational Physics, 490:112340 2023.
Florido, J., Wang, H., Khan, A., and Jimack, P. K. In- vestigating guiding information for adaptive collocation point sampling in pinns. In International Conference Computational Science, pp. 323 -373 Springer, 2024.
Owhadi, H. Multiscale Modeling & Simulation, 13(3):812-823 2015.
Bayesian numerical homogenization.
Forbberg, B. and Flyer, N. Solving PDEs with radial basis functions, Acta Numerica, 14:215-258, 2015.
Owhadi, H. and Scovel, C. Operator-Adapted Wavelets, Fast Solvers, and Numerical Homogenization: From a Game Theoretic Approach to Numerical Approximation and Algorithm Design, volume 35. Cambridge University Press, 2019.
Frostig, R., Johnson, M. J., and Leary, C. Compiling ma- chine learning programs via high-level tracing. Systems for Machine Learning, 4(9), 2018.
Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., al. Pytorch: An imperative style, high- lighpeeffrmmacc deep learning library. Advances in neural information processing systems, 32:8026-8037,, 2019.
Graepel, T. Solving noisy linear operator equations by Gaussian processes: Application to ordinary and partial differential equations. In ICML, volume 3, pp. 231-241, 2003.
Hardy, R. L.. Multiquadric equations of topography and other irregular surfaces. Journal of geophysical research, 7668):1905-1955, 1971.
Powell, M. J. The theory of radial basis function approx- imation in 1990. Advances in numerical analysis, pp. 105-210, 1992.
Izmailov, P., Novikov, A., and Kropotov, D. Scalable Gaus- sian processes with billions of inducing inputs via ten- train decomposition. In International Conference Artificial Intelligence and Statisties, pp. 726-735. PMLR, 2018.
Quinonero-Candell J. and Rasmussen, C. E. A unify- ing view of sparse approximate Gaussian process re- gression. The Journal of Machine Learning Research, 6:1939-1959, 2005.
Rahaman, N., Baratin, A., Arpit, D., Draxler, F., Lin, M., Hamprecht, F., Bengio, Y., and Courville, A. On the spec- tral bias neural networks. In International Conference on Machine Learning, pp. 5301-5310. PMLR, 2019.
Jin, R., Yang, T., Mahdavi, M., Li, Y.-F., and Zhou, Z.-H. Im- proved bounds for the nyström method with application to kernel classification. IEEE Transactions on Information Theory, 59(10):6939-6949 2013.
Rahimi, A. and Recht, B. Random features for large- scale kernel machines. Advances in neural information processing systems, 20, 2007.
Kansa, E. J. Multiquadrics- scattered data approxima- tion scheme with applications to computational fluid- dynamics-15 solutions to parabolic, hyperbolic and elliptic partial differential equations. Computers & mathematics with applications, 1968-1):147-16, 1990.
Raissi, M., Perdikaris, P., and Karniadakis, G. Machine learning of linear differential equations using Gaussian processes. Journal of Computational Physics, 348:683- 693, 2017.
Kolda, Т. G. Multilinear operators for higher-order decom- positions. Technical report, Sandia National Laboratories (SNL), Albuquerque, NM, and Livermore, CA, 2006.
Raissi, M., Perdikaris, P., and Karniadakis, G. F. 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, 2116666-777,0999 2019.
Lindgren, F., Rue, H., and Lindström, An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation ap- proach. Journal of the Royal Statistical Society Series B: Statistical Methodology, 73(4):423-498, 2011.
Saatcci, Y. Scalable inference for structured Gaussian process models. PhD thesis, Citeseer, 2012.
Long, D., Wang. Z., Krishnapriyan, A., Kirby, R., Zhe, S., and Mahoney, M. AutolP: A united framework to inte- grate physics into Gaussian processes. In International Conference on Machine Learning (ICML), pp. 14210- 14222. PMLR, 2022.
Schafer, F., Katzfuss, M., and Owhadi, H. Sparse Cholesky factorization by Kullback-Leibler Minimization. SIAM Journal on scientific computing. 43(3) 2021.
10

---

Toward Efficient Kernel- Solvers for Nonlinear PDEs
Tolstykh, A. I. and Shirobokov, D. On using radial basis functions in "finite difference mode" with applications to elasticity problems. Computational Mechanics, 33(1): 68-79, 2003.
Wang, J., Cockayne, J., Chkrebtii, O., Sullivan, T. J., and Oates, C. J. Bayesian numerical methods for nonlinear partial differential equations. Statistics and Computing,
31:1-20, 2021.
Williams, C. K. and Rasmussen, C. Gaussian processes for machine learning, volume 2. MIT press Cambridge, MA, 2006.
Wilson, A. and Nickisch, H. Kernel interpolation for scalable structured Gaussian processes (KISS-GP). In International conference machine learning, pp. 1775- 1784. PMLR, 2015.
Wilson, A. G., Dann, C., and Nickisch, H. Thoughts on massively scalable Gaussian processes. arXiv preprint arXiv:1511.01870, 2015.
Zhe, S., Xing, W., and Kirby, R. M. Scalable high-order Gaussian process regression. In The 22nd International Conference on Artificial Intelligence and Statistics, pp. 2611-2620. PMLR, 2019.

---

Toward Efficient Kernel-Bass Solvers for Nonlinear PDEs
Appendix
A Proof of Lemma 4.2
Since for any collocation point Xm. P(u*)(xm)) I (xm) when Xm E Ω and B(u") (xu))(xm) = g(Xm) when X, E 00, we can re-write 6
minimize ||u||er weld
s.t. Ma Em= (P(u)(xm) -P(u)))) (Xm)) M-Mo 2m=Ma+1 (Xm))(Xm) - B(u (xm)) S
(18)
u. takes the kernel interpolation form
Step 1. In the first step, we show that for all above some threshold (depending on h), there exists a minimizer uf for (88 and we would like also to bound the RKHS norm of ut, namely || I1x. To this end, we utilize an intermediate optimization problem,
minimize ||u||e well
(19)
s.t. u(xm) = u² (Xm), 1 m M.
Denote the minimizer of (19) by uās- This is standard kernel regression problem. According the representation theorem, "is takes the kernel interpolation form ( and ||u5 l1ee S ||u"
Since uār u² is zero at all the collocations points in S2 according to the sampling inequality (see Proposition A. of (Batlle et al. 2023) ) when the fill-distance fog is sufficiently small (note that hg S f),
|| - u (1)) S h² h =² -='' ())+)
(20)
where S means the inequality holds with a positive constant factor multiplied by the right-hand side, and the constant independent of the terms both sides. Combining with (C2) of Assumption 4.1 we obtain
IP((ur)) - ((u")llara [|B(u))) B(u" (llx S h² [[ur u [1x+++(0)
(21)
Since I is continuously embedded in Ho+7 (12) - (C3) of Assumption 4.1 we have
||u u² (())+ S Il llze-
(22)
Combining 21
22
and the fact ||u liee I ||u, we have
[P(()M) P(=")(())))))))))))) ((")")] (1)) + [[B(urr) - B(u ((("")llen (ln) S hh7ll""lee
(23)
According (C2) of Assumption 4.1 since k > and t 2 d-1 according to Sobolev embedding theorem (Adams &
Fournier 2003 Theorem 4.12), both H (S2) and H² (aN) are continuously embedded into (*) and (aN)) respectively. Therefore,
[P(())) - P(u"))lo & [(P()))) - P(u')((')))))))) P(u"))lN=())) [[(()s)) - B(" B(ux"((loroon )lco) X [[((5s)) B(u") (())))))))))
(24)
At any collocation point, we obviously have
(P(w))((xn) P(u P(u))((mm)) (xm))) [P(ut) P(''p(""))))) (1) (B(u) ( (Xm) B(u" (xm)) < [B(u)) B(u I12
(25)
Combining 23
24 and (25
we can obtain that
Mr
Σ (P(w)))(x))) - (xx))
Mn m=
+ M Ma
Σ (B(u)r(xxm)) (xuu)(xmm))
(26)
m= Ma+1
12

---

Toward Efficient Kernel- Solvers for Nonlinear PDEs
where C 0 is constant independent of h and other terms in the inequality.
The result 26p means that given the collocation points M and Ch²Tllu* l the feasible region of the optimization
problem 8 is nonempty and at least includes "iz. Therefore, the minimizer of (18 must exist and satisfy
|u+ S |us lla S (u="...
(27)
Step Next, we analyzed the error of Plut) and B(u ). For notation convenience, we define two error functions,
{r(x) (x) (x) x
EB (x) [xx))(xx (x)) x 8n.
(28)
We would like to bound the - norm of the error functions, namely, ||Ep (²) and || We first consider the case for Ép. The idea is to decompose S2 into Voronoi diagram based on the collocation points, resulting in Mg regular non-overlapping regions, Ti U... U Tase = Ω2 such that each region T/ only includes one collocation point x,, and its filled-distance h4 S (1 Sis Mn). We therefore can decompose the squared L norm as
Ma Σ}
Ma
/ (x(x) Σ
(29)
Since according to the mean inequality,
{x(x)² (fp(x) - (x(x.) + {x(x,)) s 22(Ep(x)) {x(x:)) + 2xp(xx))
we immediately obtain
IK ( S [1EP - {p(x))lIH((T)) + ((T)))((x))'
(30)
where A(Ti) is the volume of Tg.
The function & - {x(x,) takes zero at x,. Since the aspect ratio of Tj is bounded, we can apply the sampling inequality (which is also called Poincaré inequality this case),
I (x Chtller (x(x)) & h 1lE {
(31)
where C' > is a constant depending on the aspect ratio of Tr Using the mean inequality again,
P {p(x))]l
& 2*
(ILP (TA + [IE(())) 2* (IEEll3ETT) + ((T))
(32)
Since X(Ti)= $ 1rd, combining
29
30
and
32p.
we can obtain
| ( & 2k
||r l Baacer + (hd 122ktd) Er(x))
2k
Ep E + (1e + h22td) Maie,
(33)
where e comes from the constraint of (18). To ensure feasibility and to establish convergence, we set c Ch2 112 as
shown in (66 When h S M and is sufficiently small, we have (hd + 122+d) S (n + 122+d) M - 1 1:2 2
Therefore, we can extend the R.H.S of 33 to
Krll (i $ 1:24 1 11r + 2 1n
(34)
We can follow similar approach to show that
[ (1200 S [EB (oa) + h 222
(35)
Combining 34p and 335).
(IEE)l) HA9 + ())))))))))))) S i.2miiu(t.k)
(IEErll**()) + (IE)ll))()))) 2² [22T|/²1] u
(36)
13

---

Toward Efficient Kernel-Bace Solvers for Nonlinear PDEs
10
10-1
L Error
2 Error
10:2
10
10-3
0.03
0.62
0.01
0.11
0.08
0.06
0.04
0.01
h
h
(a) Burgers' equation with 0.02
(b) Burgers' equation with 0.001
Figure 1: Fill-distance (mesh norm) h L error of SKS.
Since U He+t (C3) of Assumption 4.1 have U Ha Leveraging and (C2) of Assumption 4.1
immediately obtain
[ + ((sll S 1| u² (1= S llue - u ullle S ||u"
(37)
Combining (36) and 37). and (CI) of Assumption 4.1 we arrive at
|| (1000)) X 1/llu//l
(38)
where p minkk, r). When ut converges u*.
B Proof of Proposition 4.3
The constraint optimization problem 66) is equivalent to the following mini-max optimization problem,
Mg
min max #20
(P(w)(xm. P(u") (xm))
Mn
(B(u)(xm) B(u" (xm)))
(39)
M
M9 m=Ma+1
Suppose the feasible region is non-empty. Denote the optimum of 39 by (w!,w!) w!). Then ul is minimizer of 6 Now if we set 3 = wt in and optimizing (7) will recover the minimizer u.
More Results
C.1 Point-wise Error
For fine-gra -ree-gr evaluation, we examined how the point-wise error of DAKS and SKS varies along with the increase of
collocation points. To this end, we altered the number of collocation points from 600, 4800 and 120K on Burgers' equation
with v 0.001, and from 600, 2400, and 40K on 2D Allen-Cahn equation with both 15 and The results
14

---

Toward Efficient Kernel- Solvers for Nonlinear PDEs
shown in Fig. 2 and It can be seen that across all the three PDEs, the solution error of SKS decreases and more along with the increase of collocation points. Note that for Allen-Cahn with 15, the visual difference between SKS using 2400 and 40K collocation points is little, though numerically the difference is at three orders of magnitudes (5.15E-03 vs. 4.44E-06). For DAKS, the point-wise eπrоr decreases substantially the number of collocation points grows when solving Burgers' equation (see Fig. 2 but not obviously on solving Allen-Cahn equation (see Fig. [ 4) This consistent with the global error shown in Table 2 This might be because the quantities collocation point used are not sufficient to lead to a qualitative boost DAKS. However, scaling up to much more collocation points, such as 400K, incurs substantial increase of the computational cost.
0.00
0.25
0.50
0.75
1.00
1.25
-0.5 0.0 0.5
Ground truth
DASK, 600
DAKS, 4800
SKS, 600
SKS, 4800
SKS, 120000
Figure 2: Point-wise solution error for Burgers' equation
14 with viscosity = 0.001.
0.0
0.5
1.0
1.5
Ground truth
DAKS, 600
DAKS, 2400
SKS, 600
SKS, 2400
SKS, 40000
Figure 3: Point wise solution error 2D Allen-Cahn equation
15.0.
0.00
0.25
0.50
0.75
1.00
1.25
1.50
-1
Ground truth
DAKS, 600
DAKS, 2400
SKS, 600
SKS, 2400
SKS, 40000
Figure 4: Point-wise solution error for 2D Allen-Cahn equation 17 with 20.0.
C.2 Ablation Study on Grid Shape
investigated how the grid shape can influence performance of our method. To this end, we tested on Burgers' equation with 0.02 and 0.001. For the former case, we fixed the number of collocation points be 4800 and varied the time resolution from 10 80, and the spatial resolution obtained by dividing 4800 by the time resolution and rounding up to fixed the number of the collocation points to 120K, and varied the time resolution from 100 to 600. We show the T error of using each grid shape in Table 5 It can be seen that the grid shape does influence the error. In particular, when / is small, i.e., / = 0.001, the higher the spatial resolution, the smaller the error. The smallest enror achieved when we use the space-time resolution 1200 100. The time resolution seems to have much less effect on the solution accuracy. This is reasonable, because on Burgers' equation, a smaller viscosity (v) increases the sharpness the shock wave (spatial function). Naturally, the higher the spatial resolution, the more accurate the sharpness can be captured.
an integer. Similarly, for 0.001,
15

---

Toward Efficient Kernel-Bace Solvers for Nonlinear PDEs
Table 5: L error of SKS using different grid shapes to solve Burgers' equation tine-resolution.
14. The grid shape is depicted as "spatial-resolution X
(a) Viscosity 0.02, with 4800 collocation points.
Grid shape
60 80
69 x 70
96 x 50
160 30 480 10
L error
2.27E-03
4.10E-04
7.54E-05
1.72E-04
5.98B-04
(b) Viscosity v = 0.001 with 120000 collocation points.
Grid shape
200 600
240 500
400 300
600 x 200
1200 x 100
L error
1.93E-02
6.68E-03
3.44E-03
2.28E-03
1.63E-03
Table 6: L error of solving PDEs on inreularly-shaped domains.
L Error
SKS
DAKS
PINN
Nonlinear Elliptic
8,40E-0K 8.30E-02
4.86E-05
4,20E-02
Allen-Can (a = 15)
6.06E-01
1.00E+00
In summary, we believe that in general the grid shape should be viewed influence factor in running method, which needs to be carefully selected. The appropriate choice may also connect to the intrinsic property of the PDE itself.
C.3 Irregularly-Shaped Domains
We tested on solving the nonlinear elliptic PDE 15 and the Allen-Cahn equation 17p with = 15. For the nonlinear elliptic PDE, the domain is an inscribed circle within [0, 1] x [0, 1J. For the Allen-Cahn equation, the domain is a triangle with vertices at at (0, 0), (1, 0) and (0.5, 1). The solution is prescribed as our paper, with boundary conditions derived from the solution. For both PDEs, our method (SKS) used virtual grid on [0, 1] x [0, 1] that covers the domain. For DAKS and PINN, we sampled the same number of collocation points from the domain. For fair comparison, all the methods used the same set of 192 uniformly sampled collocation points the boundary. The error of each method given in Table 6 The point wise error is shown in Fig. 5 As see, irregularlysshapee domains, method SKS still obtains reasonably good accuracy for both cases (note that the Allen-Cahn case is much more challenging).
Allen cahn (a=15)
15
05
15
15
-05
18
10
05
Ground truth
14
PINN
62
DAKS
SKS
10
Nonlinear Elliptic
655
0.55
0.03
632
02
0.02
Ground truth
PINN
DAKS
SKS
Figure 5: Point-wise solution error for Nonlinear Elliptic and 2D Allen-Cahn equation
with 15.0 on irregularly-hhaped domains.
16

---

Toward Efficient Kernel-Basce Solvers for Nonlinear PDEs
Table 7: L error of SKS using different length scales.
(a) Nonlinear Elliptic PDE (TS).
Length-scans
0.05
0.1
0.2
0.3
I error
1.19E-02
6.80E-06
4.62E-05
9.14E-04
(b) Allen-Cahn Equation (17)
where = 15.
Length-scale
0.08
0.06
0.04
0,02
L error
4. 4.98P-01
3.83E-01
5.15E-03
2. 2.19F-01
Table 8: Sensitivity of PINN to the network width and depth.
(a) L error for 1epth=10 and layer width from 10 to 50.
Layer width
10
20
30
40
50
Nonlinear Elliptic
1.75F-01
7.68F-02
7.01
3.00F-02
4.06F-02
Allan-cahn (a 15)
6.31E+00
5.61E+00
7.08F+00
7.65F+00
1.03E+01
(b) L error for depth= 10 and layer width from 60 to 100.
Layer width
60
70
80
90
100
Nonlinear Elliptic 1.66E-02
3.58E-02
2.74E-02
2.32E-02
5.01E-02
Allan-cahn (a = 15)
1.30E+01
1.26E+01
1,20+001
1.34E401
1.09E401
(c) L error for layer width=10 and depth from 3
10.
Depth
3
5
10
Nonlinear Elliptic
3.78E-02
7.33F-02
7.01F-02 7.08F+00
6.57F-02
Allan-cahn
15)
1.13F+01
1.35F+01
1.10F+01
(d) L error for layer width=10 and depth from 15 to 30.
Depth
15
20
25
30
Nonlincar Elliptic
9.16E-02
1.13E-01
I 1.50E-01
7.49E-02
Allan-cahn (a = 15)
1.14E+00
9.30E+00
1.14E+00
9.35E+00
C.4 Sensitivity to Hyper-Parametee
To examine the sensitivity to the choice of kernel parameters, we run SKS to solve nonlinear Elliptic PDE (15) and Allen-Cahn equation (17 with a 15, with varying set of length-scale parameters. For the nonlinear elliptic PDE, we employed the grid of size 35 x 35 while for the Allen-Cahan equation, we used the grid of size 49 x 49. The results given in Table 7 As we can see, different length-scale parameters results in changes of orders of magnitude in the solution error. For example, switching the length-scale from 0.05 to 0.1, the L error for solving the nonlinear elliptic PDE decreases from 1.22 to 6.80E-06. Hence, our method is sensitive to the choice of the kernel parameters.
As comparison, we also examined the sensitivity of PINN to the choice of the architectures. To this end, we fixed the depth 8 and varied the layer width from 10 100, and also fixed the width 30 while varying the depth from 30. We then tested PINN solving the same PDEs. The solution error is reported in Table One can see that when layer width is greater than 50 or the depth is beyond 8, there is no significant improvement in performance. The accuracy remains within the same magnitude with only minor variations. However, larger networks lead to substantially increased computational costs for PDE solving.
C.5 Optimizing Basis Coefficients
To investigate the difference with RBF methods, we tested on optimizing the basis coefficients, namely a = Kn1n, instead of the solution values 77, for our model. We employed the same optimization approach. We tested on solving Burger's s equation with / = 0.02. The results are reported in Table 10 We can see that the relative L error increases significantly as compared to learning 77, worsening by an order of magnitude when using 1200 or 2400 collocation points, and by over 50% when using 4800 collocation points. This suggests that optimizing the coefficients a is more challenging, as each
17

---

Toward Efficient Kernel- Solvers for Nonlinear PDEs
Table 9: Runtime in seconds with respect to the number of collocation points. SKS-Naive means running SKS using naive full kemel
matrix computation, without leveraging the Kroneeker product structure (see Section 4 Note that N/A means the method not able with the corresponding number of collocation points.
(a) The Burgers' equation
14 with v 0.001.
Method
2400
4800
43200
SKS (per-iter)
4.6E-4
9.8E-4
6.8E-3
SKS-Naive (per-iter)
1.4E-02
5.4E-02
N/A
DAKS (per-iter)
7.43
38.5
N/A
PINN (per-iter)
2.7B-1
5.2E-1
4.1E-1
SKS (total)
22.15
94.25
2007.1
DAKS (total)
89.14
462.18
N/A
PINN (total)
706.24
721.94
6454.7
(b) Allen-Cahn equation (17
with
15.
Method
2400
4800
6400
8100
22500
SKS (per-iter)
3.6E-4
9.1E-4
1.2E-3
1.8E-3
5,9E-3
SKS-Naive (per-iter)
1.1E-2
4.3E-2
7.2E-2
1.EE-
N/A
DAKS (per-iter)
2.1
10.5
N/A
N/A
PINN (per-iter)
5.6E-2
1E-1
1,3E-
1,5E-
4.3E-
SKS (total)
27.1
99,56
116.8
132.57
474.34
DAKS (total)
16.44
84.18
N/A
N/A
N/A
PINN (total) 2821
5112
6287
7614
21375
Table 10: L error of learning vs. ar on solving Burgers' equation with v = 0.02.
Method Optimizing 77
1200
2400
4800
5.40E-03
7.83E-04
3.21E-04
Optimizing coefficients --11
2.80E-02
1.56E-03
7.14E-04
coefficient exerts a strong global influence on the entire solution approximation, including both boundary conditions and internal regions. We empirically observed that conflicts often arise between fitting the boundary conditions and minimizing the PDE residuals, resulting in less effective optimization outcomes.
C.6 Effectiveness Regular Grid Points for DAKS
The original DAKS (Chen et al. 2021) utilized randomly sampled collocation points. In this study, we evaluated the performance DAKS using regular grid points as collocation points to solve Burgers' equation with v = 0.02. The evaluation was conducted on exactly the same grids used our method (SKS). As shown in Table regular grid points proved less effective than random collocation points for DAKS.
D Limitation and Discussion
Currently, the most effective training for SKS is fulfilled by stochastic optimization, namely ADAM. We need to run a large number of ADAM epochs to achieve promising solution accuracy. It means that the PDE solving procedure is slow. The second-order optimization methods, such as L-BFGS, neither improve the solution accuracy accelerate the convergence. We have also tried the relaxed Gauss-Newton approach used in DAKS. However, this method can only achieve good performance on the nonlinear elliptic PDE, and easily diverges on the other cases. This might stem from that we take derivatives (or other linear operators) over the kernel interpolation form, which makes the convergence of the fixed point iterations used in DAKS much more difficult. We plan to develop novel Gauss-Newton relaxations to ensure convergence and stableness (at least in practice) so that we can further accelerate the PDE solving.
18

---

Toward Efficient Kernel-Basee Solvers for Nonlinear PDEs
Table 11: L error of different sampling methods of DAKS. Results are based solving Burgers with 0.02.
Sampling method
600
1200
2400
4800
Grid
4.09E-01
3,85E-01
4.27E-02
5,67E-02
Random
1.75E-02
7.90E-03
8.65E-04
9.76E-05
19


Reviews:
[{"review_id": "fbwhaR7rSr", "reviewer": "Reviewer_P51e", "summary": "The authors propose an assymetric RBF collocation method for solving general PDEs from a Gaussian process/RKHS point of view. They parametrize the solution $u(x;\\eta)$ as a (Gaussian RBF) kernel interpolant given function values on collocation points, then differentiate this representation to obtain derivatives values at collocation points, which they then optimize to approximately satisfy the PDE on collocation points using a least squares formulation. \n\nThey then take advantage of the fact that under the assumption that the kernel $k$ factorizes across different dimensions of the input space, then both $\\eta \\to u(x;\\eta)$ aand $\\eta \\to Lu(x;\\eta)$ posses efficient efficient evaluation formulas via Kronecker-factored matrix vector products.  \n\nWith the parametrization and fast evaluation in hand, their algorithm consists of minimizing a regularized squared residual of PDE mismatch on interior collocation points and boundary mismatch using stochastic optimization (ADAM).", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "The claims made in the paper seem well substantiated by evidence. I am only unsure as to why SKS would outperform DAKS--the theory for DAKS seems stronger, but a good implementation is difficult, and especially with the RBF kernel, matrices can get very ill-conditioned, making it hard to disambiguate discretization error from numerical roundoff error.", "methods_and_evaluation_criteria": "The problems and evaluation criteria seem good. I am a little bit skeptical of ADAM being the best solver for the problem, compared to a quasi-Newton or Newton-Krylov method, but their experiments corroborate this choice.", "theoretical_claims": "The proofs look sound. I read the first part of the proof of Lemma 4.2 carefully, and it looks all correct except for a minor typo in equation (22) which seems like its should be $u_m^*$ on the left.", "experimental_designs_or_analyses": "The experiments seem valid", "relation_to_broader_scientific_literature": "The main contributions are as follows:\n1. Convergence analysis of least squares methods for assymetric collocation--a bound on the attainable squared residual norm, and convergence rates of such least squares solutions to a true solution of the PDE. \n2. Taking advantage of the Kronecker factored structure of the kernel matrices when the kernels factorize to obtain fast evaluation formulas and efficient optimization\n3. Interpretation of RBF collocation methods as optimization over function values that are interpolated rather than coefficients of an expansion.", "essential_references_not_discussed": "None that I know of.", "other_strengths_and_weaknesses": "The strengths include a nice analysis of feasibility for the least squares problem, and a simple algorithm for general PDEs. \n\nThe sensitivity to hyperparameters seems to be a weakness (lengthscale and nugget).", "other_comments_or_suggestions": "The term \"Gaussian-Newton\" seems odd to refer to \"Gauss-Newton\" methods. \n\nThe implementation on irregular domains is a bit unclear. Throughout the rest of the paper, it is assumed that the basis points are the same as collocation points, but this seems hard to do with irregular domains (especially for the boundary. This should be clarified. \n\nThe authors should include an error plot as a function of the mesh norm $h$ to show the order of convergence of the algorithm (using Gaussian RBF should theoretically give spectral convergence?). \n\nThe authors should clarify assumption C3 in relation to their choice of Gaussian RBF kernel. Is it true that the RKHS associated to a Gaussian RBF kernel can be continuously embedded in $H^k$ for every $k$?", "questions_for_authors": "What was the batch size and hyperparameters used in the stochastic optimization?\n\nWhat did the results with attempting to apply LBFGS to improve the solutions look like? In the PINN literature, it seems that it often significantly improves the solution after running ADAM to warm up.\n\nDid the gradient norms of the objective go to zero at the end (did you reach a minimum, or stall out due to ill conditioning)?", "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": "Ke5IOFNrNX", "reviewer": "Reviewer_fQAE", "summary": "The paper proposes a new twist to a kernel-based solver for PDEs. It builds on previous work by relaxing the constraints associated to the PDE to facilitate optimization. By placing the collocation points on a grid and using a decomposable kernel, the inversion of a large matrix is broken down into many inversions of small matrices, which improves the performance. A convergence analysis is provided, showing that this solver enjoys similar properties as previous work i.e. the Sobolev norm of the residuals goes to zeros as the grid uniformly covers the input space. Numerical experiments highlight the benefits of the approach and compare to competitors (PINNs and the former kernel-based solver).\n\n## update after rebuttal\n\nThe authors clarified my question about how the constraint was handled using lagrangian duality. This and the rebuttal to other reviews convinced me to increase my score.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your valuable and insightful comments.\n\n>The introduction lacks references. Only the direct relevant work is cited (2 papers)\n\nR1: Thank you for the helpful suggestion. We will include additional references in the introduction to provide broader context and better situate our work within the existing literature from the outset.\n\n>no standard deviations\n\nR2: Great comment. In our experiments, we observed that despite using stochastic optimization (ADAM), our method (SKS) consistently converge to the same solution. So is PINN. We ran multiple trials on each PDE and found that the standard deviation of the error  is negligible. For example, below we show  the standard deviation of the $L^2$ error for our method when solving Burgers' equation ($\\nu = 0.01$) and the nonlinear elliptic equation. Given the extremely low variance, we chose to omit standard deviation values from the main results (Note that the finite difference method is deterministic and does not exhibit variance). We will clarify this point in our paper.\n\nBurgers($\\nu = 0.01$):\n\n| 600                  \t| 1200                 \t| 2400                 \t| 4800            \t|\n|---------------------|----------------------|----------------------|-----------------|\n| 1.44E-02 $\\pm$ 1.52E-07 \t| 5.40E-03 $\\pm$ 1.39E-07 \t| 1.12E-03 $\\pm$ 7.76E-08 \t| 3.21E-04 $\\pm$ 0.0 \t|\n\nNonlinear elliptic:\n\n| 300             \t| 600             \t| 1200            \t| 2400            \t|\n|-----------------|-----------------|-----------------|-----------------|\n| 1.26E-02 $\\pm$ 0.0 \t| 6.93E-05 $\\pm$ 0.0 \t| 6.80E-06 $\\pm$ 0.0 \t| 1.83E-06 $\\pm$ 0.0 \t|\n\n>I think that there is a problem with the formulation of (7). As it is, $\\epsilon$  has no impact on the problem thus I hardly see how it can be equivalent to (6)... A typical way to do this is to square this quantity again.\n\nR2: Thank you for your insightful question. Here is our clarification.\n\nFirst, **$\\epsilon$ plays an important theoretical role** (though in practice, it can often be set to zero). Our convergence proof and rate estimate (see Lemma 4.2 and Appendix Section A) are established by systematically varying $\\epsilon$. Specifically, we set $\\epsilon = C_0 h^{2\\tau}$ and vary $\\epsilon$ by adjusting the fill-in distance $h$. Each choice of $\\epsilon$ defines a distinct instance of problem (6), with its own solution. We prove that as $\\epsilon \\to 0$ (i.e., $h \\to 0$), the solution of (6) converges to the ground-truth solution of the PDE (see Appendix Section A for details).\n\nSecond, to show that solving (7) (with appropriately chosen $\\alpha$ and $\\beta$) is equivalent to solving (6), we use the **Lagrangian formulation** of (6). It is a standard result that constrained optimization problems can be reformulated as **mini-max problems over the Lagrangian**. The Lagrangian is a *linear* combination of the objective and the constraints — it does **not** involve squaring the constraints. From this perspective, it is straightforward to show that when $\\alpha$ and $\\beta$ are selected as the optimal dual variables in the mini-max problem, solving (7) yields the optimal $u$, and thus is equivalent to solving (6). This equivalence holds for **any** $\\epsilon$, not just in the limit as $\\epsilon \\to 0$. The full proof is provided in Appendix Section B.\n\nFinally, we agree that an alternative approach --- minimizing the objective plus the squared constraints (i.e., a penalty method) --- is also a good and viable idea. However, this method is primarily justified for equality constraints (i.e., $\\epsilon = 0$) (While it is possible to design penalty terms for inequality constraints, they typically introduce non-differentiable terms.) Moreover, equivalence to the original constrained problem is only guaranteed as as **the penalty weights tend to infinity**. Therefore, we believe our formulation in (7) offers stronger theoretical guarantees and greater practical flexibility when approximating solutions of (6). We will include a more detailed discussion of this alternative approach in our paper.\n\n>Could you include a discussion about the tuning of two regularization parameters (7)?\n\nR3: Thank you for the great suggestion — we completely agree. In our experiments, we selected $\\alpha$ and $\\beta$ in Equation (7) from a wide range:\n$[10^{-2}, 10^{-1}, 1, 10, 10^2, 10^3, 10^4,10^5, 10^6, 10^7, 10^8, 10^{10}, 10^{12}, 10^{14}, 10^{15}, 10^{20}]$,\njointly with other hyperparameters, including the kernel length scales and nugget terms (see their ranges in Lines 284-297 right). To efficiently tune hyperparameters, we first performed a random search to identify a promising group of hyperparameters. We then fixed all other hyperparameters and conducted a grid search over $\\alpha$ and $\\beta$. Finally, we fixed $\\alpha$ and $\\beta$ and performed a grid search over the remaining hyperparameters. We will include this discussion in the paper.", "claims_and_evidence": "Yes, the claims are valid and supported by evidence, both theoretical and empirical. The performances are not uniformly superior to competitors but I do not think it is an issue. The wording in the conclusion (\"encouraging performances\") might be closer to reality that within the abstract (\"demonstrate the advantage\").", "methods_and_evaluation_criteria": "Yes, the methods make sense. They seem to be standard in this area of research and have been investigated in other papers.", "theoretical_claims": "To the best of my knowedge, the theoretical claims are valid.", "experimental_designs_or_analyses": "I have no issue with the experimental design.", "relation_to_broader_scientific_literature": "The paper proposes a novel way to encode the constraints for solving PDEs, improving previous work. The convergence analysis is similar to what was done for previous work as well. Overall, the paper is very close to their main source of inspiration.", "essential_references_not_discussed": "I would have liked to see discussed [Learning Partiel Differential Equations in Reproducing Kernel Hilbert Spaces, George Stepaniants, JMLR 2023] and references therein. Overall I think that the overview in the introduction lacks references about related approaches beyond the scope of [Chen et al. 2021].", "other_strengths_and_weaknesses": "+: The paper is relevant to the ICML community.\n\n-: The introduction lacks references. Only the direct relevant work is cited (2 papers) which makes it hard to situate the paper at first.\n\n-: There are no standard deviations on the performance metrics", "other_comments_or_suggestions": "I think that there is a problem with the formulation of (7). As it is, $\\epsilon$ has no impact on the problem thus I hardly see how it can be equivalent to (6). The regularization term should grow as $(P(u)(x) - f(x))^2$ is away from $[0, \\epsilon]$. A typical way to do this is to square this quantity again.", "questions_for_authors": "The proposed formulation for handling the constraints relies on two additional regularization parameters. While you demonstrate the equivalence between formulations (6) and (7), this is only true for the perfect regularization strength, which we do not know. Could you include a discussion about the tuning of these parameters ?", "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": "uee8Kucxg9", "reviewer": "Reviewer_fonM", "summary": "This paper introduces a novel kernel learning framework for efficiently solving nonlinear partial differential equations (PDEs). Unlike existing methods that embed differential operators within kernels, this approach eliminates these operators from the kernel, using standard kernel interpolation to model the solution. By differentiating the interpolant, the method avoids the need for complex Gram matrix construction, which simplifies implementation and enables efficient computation. The framework leverages Kronecker product structures for scalable computation, allowing it to handle large numbers of collocation points. The authors provide a rigorous convergence analysis and demonstrate the method's efficacy on several benchmark PDEs, showing competitive performance and scalability, particularly in challenging scenarios requiring dense grids.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the valuable and constructive comments.\n\n>C1: Can you provide more details on the virtual grid on an irregular domain, which is important for the application of your method?\n\nR1: Great question. Our virtual grid is chose as the smallest rectangular grid that fully covers the irregular domain. Specifically,  for the  nonlinear elliptic PDE (see Appendix C.3 and Figure 4), the domain is a circle inscribed within $[0, 1]\\times [0, 1]$, and we hence set the virtual grid on $[0, 1]\\times [0, 1]$ with a grid size of **$50 \\times 50$**. For the Allen-Cahn equation (see Appendix Figure 4), the domain is a triangle with vertices at (0, 0), (1,0) and (0.5, 1). To fully cover this domain, we again place the virtual grid on $[0, 1] \\times [0, 1]$, using a grid size of $50 \\times 50$. We will include these details in our paper.\n\n> C2: Can you provide the running time comparison with the finite difference method?\n\nR2: Great suggestion. Below, we provide the time comparison results in seconds (for the Allen-chan equation with $a=15$).\n\n| \t                           | 2400     \t | 4800     \t | 6400     \t | 8100    | 22500   |\n|-----------------------------|------------|------------|------------|---------|---------|\n| **FD (per-iter)**             \t | 1.6E-02 \t  | 1.4E-02 \t  | 1.4E-02    | 1.4E-02 | 1.5E-02 |\n| SKS (per-iter)              | 3.6E-4     | 9.1E-4     | 1.2E-3     | 1.8E-3  | 5.9E-3  |\n| DAKS (per-iter)             | 2.1        | 10.5       | N/A        | N/A     | N/A     |\n| PINN (per-iter)             | 5.6E-2     | 1E-1       | 1.3E-1     | 1.5E-1  | 4.3E-1  |\n|    |      |        |      |      |      |\n| **FD (total)**\t                 | 0.13 \t     | 0.14 \t     | 0.16 \t     | 0.18    | 1.15    |\n| SKS (total)\t                | 27.1 \t     | 99.56 \t    | 116.8 \t    | 132.57  | 474.34  |\n| DAKS (total)\t               | 16.44 \t    | 84.18 \t    | N/A \t      | N/A     | N/A     |\n| PINN (total)\t               | 2821 \t     | 5112 \t     | 6287 \t     | 7614    | 21375   |\n\nAs we can see, the finite difference (FD) is computationally much more costly per iteration compared to our method (SKS). This might be due to the expensive cost of computing the inverse Jacobian during the root finding procedure. However, FD typically converges within just a few dozen iterations — substantially faster than the stochastic optimization used in both SKS and PINN — leading to a much lower total runtime. We will include this comparison in the paper.\n\n> C3: In future work, can the proposed method be extended to the neural operator setting for further acceleration?\n\nR3: We appreciate the reviewer for bringing up this excellent idea — we completely agree. In fact, we have already begun exploring this direction in our ongoing work. One line of effort involves replacing the trunk network in the Deep Operator Network (DeepONet) with Gaussian process bases. This enables us to adopt a similar computational strategy to accelerate both training and prediction, while also providing a natural framework for uncertainty quantification. We are also extending our approach to deep kernel-based operator learning, where our method can further accelerate function transformations in the latent channel space. We look forward to continuing our exploration in this promising direction.", "claims_and_evidence": "The main claim is that the newly proposed method has better computation efficiency and scalability, which are well supported by the experiments.", "methods_and_evaluation_criteria": "The kernel method and PDE solution accuracy criteria are suitable.", "theoretical_claims": "I roughly checked the correctness of theoretical claims.", "experimental_designs_or_analyses": "The description and analyses of experiments are valid.", "relation_to_broader_scientific_literature": "Kernel Method, Numerical method of PDEs.", "essential_references_not_discussed": "No", "other_strengths_and_weaknesses": "The main strength of the proposed method is the acceleration of kernel-based PDE solvers.", "other_comments_or_suggestions": "Readability can be improved by adding: 1）intuitive and informal description of math derivations in Sections 2 and 3,  2) bullets in Section 5.", "questions_for_authors": "1. Can you provide more details on the virtual grid on an irregular domain, which is important for the application of your method?\n2. Can you provide the running time comparison with the finite difference method?\n3. In future work, can the proposed method be extended to the neural operator setting for further acceleration?", "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": "woBuWo72A8", "reviewer": "Reviewer_DKwk", "summary": "This paper presents a fair contribution to kernel-based PDE solvers for nonlinear PDEs, improving upon prior methods by eliminating the need for differential operator-embedded kernels. The proposed algorithm enhances computational efficiency by leveraging Kronecker product properties and avoiding complex Gram matrices. The paper also provides convergence proofs and error rate analysis under regularity assumptions. The proposed method is evaluated on Burgers’, nonlinear elliptic, Eikonal, and Allen-Cahn equations, showing its comparable accuracy and improved scalability compared to the baselines.\n\n## update after rebuttal\nThank you for addressing the concerns. I am keeping my original score.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your valuable and constructive comments.\n\n>The method assumes a structured grid, which limits its generalization to unstructured meshes and other complex discretization routines.\n\nR1: We appreciate your insightful comments. Indeed, we agree that our method relies on a structured grid and does not directly apply to more complex discretizations. A simple yet effective workaround — which we have validated (see Section 6 and Appendix Section C.3) — is to employ a virtual grid that covers the irregular domain.\n\nIn future work, we plan to explore two directions to better generalize our method. First, we will investigate domain decomposition and hierarchical grid construction to adaptively adjust the grid resolution across local regions while preserving computational efficiency. Second, we aim to learn a mapping from unstructured mesh points to a latent grid, where our efficient kernel solvers can be applied. This approach is inspired by recent work such as GeoFNO[1] in the neural operator literature.\n\nWe will include this discussion in the paper.\n\n[1] Li, Z., Huang, D. Z., Liu, B., & Anandkumar, A. (2023). Fourier neural operator with learned deformations for pdes on general geometries. Journal of Machine Learning Research, 24(388), 1-26.\n\n>The hyperparameter selection process (e.g., kernel length scales) is not discussed in detail.\n\nR2: Thank you for your great suggestion. Given the wide range of hyperparameters --- including $\\alpha$ and $\\beta$ (see our response R3 to Reviewer fQAE), as well as kernel length-scales and nugget values (see Lines 284–297 right column) --- performing a full grid search would be prohibitively expensive. Therefore, we adopt a hybrid strategy. We begin with a random search to identify a promising set of hyperparameters. Then, we perform a grid search over $\\alpha$ and $\\beta$, keeping the other parameters fixed. Once $\\alpha$ and $\\beta$ are selected, we fix them and conduct a grid search over the remaining hyperparameters, including the nugget and length-scales. We will include a detailed description of this procedure in the paper.\n\n>I'm also curious how the runtime of the Kronecker based method compares with the naive full matrix method.\n\nR2: Great question. Below, we provide the runtime of our model using the naive full matrix computation (i.e., without exploiting the Kronecker product structure). As shown, the per-iteration runtime with naive matrix operations is consistently around **100× slower** compared to our method, which leverages Kronecker product properties. Furthermore, when the number of collocation points increases to 22,500 for the Burgers’ equation and 43,200 for the Allen–Cahn equation, the naive approach exceeds available memory, resulting in out-of-memory (OOM) errors — rendering it infeasible to run. We will include these results in our paper.\n\n| \t        Allen-Cahn ($a=15$)              | 2400     \t | 4800     \t | 6400     \t | 8100    | 22500 |\n|-------------------------------------------|------------|------------|-----------|---------|-------|\n| **Naive computation** (per-iter)            \t | 1.1E-02 \t  | 4.3E-02 \t  | 7.2E-02   | 1.1E-01 | OOM   |\n| SKS (per-iter)            \t               | 3.6E-4     | 9.1E-4     | 1.2E-3    | 1.8E-3  | 5.9E-3  |\n| DAKS (per-iter)            \t              | 2.1        | 10.5       | OOM          | OOM     | OOM     |\n| PINN (per-iter)            \t              | 5.6E-2     | 1E-1       | 1.3E-1    | 1.5E-1  | 4.3E-1  |\n\n| \t  Burgers' $\\nu=0.001$                  | 2400     \t | 4800     \t | 43200     \t |\n|------------------------------------------|------------|------------|-------------|\n| **Naive computation** (per-iter)           \t | 1.4E-02 \t  | 5.4E-02 \t  | OOM         |\n| SKS (per-iter)           \t               | 4.6E-04 \t  | 9.8E-04 \t  | 6.8E-03     |\n| DAKS (per-iter)           \t              | 7.43 \t     | 38.5 \t     | OOM         |\n| PINN (per-iter)           \t              | 2.7E-01 \t  | 5.2E-01 \t  | 4.1E-01     |", "claims_and_evidence": "Claims are largely supported by theoretical derivations and experiments.", "methods_and_evaluation_criteria": "- Strengths: The methods and evaluation criteria are in general well-aligned with PDE-solving tasks. The benchmarks involve multiple nonlinear PDE formulations, and the baselines include DAKS, PINNs, and legacy finite difference methods too. The paper also offers scalability analysis, which is appreciated. \n- Weakness: The method assumes a structured grid, which limits its generalization to unstructured meshes and other complex discretization routines.", "theoretical_claims": "This paper provides convergence analysis, proving that the method maintains error bounds similar to prior kernel PDE solvers despite using a smaller model space. Other claims are also briefly reviewed and look good, but not checked in detail.", "experimental_designs_or_analyses": "- Strengths: The experimental setup is in general well-designed. It covers multiple PDE formulations, varying levels of difficulty (small vs. large collocation points), and different kernel settings. \n- Weaknesses: The hyperparameter selection process (e.g., kernel length scales) is not discussed in detail. I'm also curious how the runtime of the Kronecker based method compares with the naive full matrix method.", "relation_to_broader_scientific_literature": "This paper is built upon kernel-based PDE solvers and relates to Gaussian Process models for PDE solving. The Kronecker product approach also aligns with structured kernel methods. There is sufficient technical contribution.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "N/A", "other_comments_or_suggestions": "N/A", "questions_for_authors": "See the \"Weaknesses\" items above.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: The paper introduces a kernel learning framework for solving nonlinear partial differential equations efficiently. The key idea is to pose the problem as a kernel interpolant by relaxing the constraints. By placing the collocation points on a grid and using a decomposable kernel, the inversion of a large matrix is broken down into many inversions of small matrices, which improves the performance. Theoretical and experimental results are presented to demonstrate the efficacy of the proposed approach.

All the reviewers are positive about the paper, and the contribution of a fast kernel learning framework is interesting and useful. The rebuttal to the authors' comments appears satisfactory. I strongly urge the authors to incorporate the reviewer comments and the rebuttal discussion in the camera-ready version. Based on the reviewer comments, I propose to accept the paper.
Decision: Accept (poster)