Toward Efficient Kernel-Based Solvers for Nonlinear PDEs

## Abstract

Abstract This paper introduces 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.

## 1 Introduction

Solving partial differential equations (PDEs) stands as a central task in scientific and engineering domains.
Recently, machine learning (ML)-based solvers have garnered significant attention. Unlike traditional numerical methods, ML-based solvers eliminate the need for complex mesh designs and intricate numerical techniques, enabling simpler, faster, and more convenient implementation and use. Among these solvers, kernel methods or Gaussian processes (GPs) (Williams and Rasmussen,, ) hold great promise due to their solid mathematical foundations, offering high expressiveness, robustness, and the ability to quantify and reason under uncertainty.
Recently, introduced a general kernel method to approximate solutions of nonlinear PDEs. They augment the representation by incorporating differential operators (more generally, linear operators) into the kernels, and jointly fit the solution values and their derivatives of the PDE on a set of collocation points. Their approach has demonstrated promising performance in solving several benchmark nonlinear PDEs, backed by a rigorous error analysis, including both convergence and convergence rates (; Batlle et al.,, ).

Despite its success, the methodology requires manual construction of a Gram matrix between the solution and its derivatives (which show up in the PDEs). It enumerates combinations between the operators on the kernel to compute different sub-blocks. The process enlarges the size of the Gram matrix (as compared to the number of collocation points), and the computation becomes challenging with a large number of collocation points, a crucial factor for capturing complex PDE solutions such as with potential roughness (or even non-smoothness) and higher frequencies.

In response, this work proposes an alternative kernel-based framework for solving nonlinear PDEs with several key contributions:

- •
Framework: We remove the basis functions associated with the differential evaluation functionals in the approximation and use a standard kernel interpolation for solution modeling. To approximate solution derivatives, we directly differentiate the interpolant. By minimizing the RKHS norm along with a boundary and residual loss, we estimate the solution without the need for manually constructing a complex Gram matrix. Implementation is straightforward and convenient with the aid of modern automatic differential libraries.
- •
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 a grid and employ a product kernel over each input dimension. This choice induces a Kronecker product structure within both the kernel interpolation and its differentiation, effectively bypassing the need for computing the entire Gram matrix. Such modification results in a substantial reduction in computational costs, enabling the efficient processing of tens of thousands or even millions of collocation points. This is achieved without resorting to any sparse and/or low-rank approximations.
- •
Theorem: We provide a rigorous analysis of our framework. We show the convergence and convergence rate of our method under appropriate PDE stability and regularity assumptions that are very similar to the assumptions used in the prior work (Batlle et al.,, ). However, our results are not a trivial extension of the prior work in that while our framework uses a reduced model space for efficient computation, our convergence results are as comparably strong as those in the prior work (; Batlle et al.,, ) that employs a richer model space. This is achieved through a much more sophisticated proof. We construct an interpolate of the true solution as an intermediate bridge connecting 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, sampling inequality and mean inequality, we are able to bound the $L_{2}$ norm of the error w.r.t the PDE operators. Combined with the bounded RKHS norm of our approximation, we establish convergence and convergence rate. The results theoretically affirm the efficacy of our method in yielding accurate solutions.
- •
Experiments: Evaluation on Burgers’, nonlinear elliptical, Eikonal, and Allen-Cahn equations have demonstrated our method’s efficacy. For less challenging scenarios where a small number of collocation points, e.g., $1000$, is sufficient, our method achieves comparable or sometimes smaller errors than the existing methods. In more challenging scenarios, such as Burgers’ equation with a viscosity of $0.001$, our method seamlessly scales to tens of thousands of collocation points, yielding low errors on the order of $10^{-3}$ to $10^{-6}$, underscoring its robustness and accuracy across a range of problem complexities.

## 2 Background

We consider a PDE of the general form,

$$ $\displaystyle{\mathcal{P}}(u)$ $\displaystyle=f({\bf x}),\;\;{\bf x}\in\Omega,$ $\displaystyle{\mathcal{B}}(u)$ $\displaystyle=g({\bf x}),\;\;{\bf x}\in\partial\Omega,$ (1) $$

where ${\bf x}=(x_{1},\ldots,x_{d})^{\top}$, and ${\mathcal{P}}$ and ${\mathcal{B}}$ are nonlinear differential operators in the interior $\Omega$ and boundary $\partial\Omega$, respectively. We assume ${\mathcal{P}}$ and ${\mathcal{B}}$ are composed by a series of linear operators, namely,

$$ $\displaystyle{\mathcal{P}}(u)({\bf x})$ $\displaystyle=P\left(L_{1}(u)({\bf x}),\ldots,L_{Q_{\Omega}}(u)({\bf x})\right ),\;{\bf x}\in\Omega,$ $\displaystyle{\mathcal{B}}(u)({\bf x})$ $\displaystyle=B\left(L_{Q_{\Omega}+1}(u)({\bf x}),\ldots,L_{Q}(u)({\bf x}) \right),\;{\bf x}\in\partial\Omega,$ $$

where $P(\cdot)$ and $B(\cdot)$ are nonlinear functions, each $L_{j}$ ($1\leq j\leq Q$) is a linear operator, such as derivatives of $u$ and their linear combinations: $\partial_{x_{1}x_{1}}u$, $\partial_{x_{1}x_{2}}u$, $a\cdot\partial_{x_{2}x_{2}}u+b\cdot u$, etc.

To solve the PDE (1), proposed to sample a set of collocation points, ${\mathcal{M}}=\{{\bf x}_{1},\ldots,{\bf x}_{M_{\Omega}}\in\Omega,{\bf x}_{M_{
\Omega}+1},\ldots,{\bf x}_{M}\in\partial\Omega\}$ in the domain, and estimated the solution values and all the relevant linear operators over the solution, $\{L_{j}(u)(\cdot)\}_{j}$, evaluated at the collocation points. Specifically, a nested optimization problem was formulated as

$$ $\displaystyle\underset{{\bf z}}{\text{minimize}}\begin{cases}\underset{u\in{ \mathcal{U}}}{\text{minimize}}\;\;\;\|u\|_{{\mathcal{U}}}\\ \text{s.t.}\;\;L_{j}(u)({\bf x}_{m})=z^{j}_{m},1\leq j\leq Q,\\ 1\leq m\leq M_{\Omega}\;\;\text{if }j\leq Q_{\Omega}\\ \text{otherwise }1\leq m\leq M-M_{\Omega}\end{cases}$ $\displaystyle\text{s.t.}\;\;\;P(z^{1}_{m},\ldots,z^{Q_{\Omega}}_{m})-f({\bf x} _{m})=0,\;\;1\leq m\leq M_{\Omega},$ (2) $\displaystyle B(z^{Q_{\Omega}+1}_{m^{\prime}},\ldots,z^{Q}_{m^{\prime}})-g({ \bf x}_{m^{\prime}})=0,\;\;1\leq m^{\prime}\leq M-M_{\Omega},$ $$

where ${\bf z}$ includes all $\{z^{j}_{m}\}$, ${\mathcal{U}}$ is the Reproducing Kernel Hilbert Space (RKHS) associated with a kernel $\kappa(\cdot,\cdot)$.With the RKHS ${\mathcal{U}}$, the similarity (or covariance) between any $z^{j_{1}}_{m_{1}}$ and $z^{j_{2}}_{m_{2}}$ in ${\bf z}$ is

$$ $\displaystyle c(z^{j_{1}}({\bf x}_{m_{1}}),z^{j_{2}}({\bf x}_{m_{2}}))=L_{j_{1 }}\circ L_{j_{2}}(\kappa)({\bf x}_{m_{1}},{\bf x}_{m_{2}}),$ $$

where $L_{j_{1}}$ is applied along the first argument of $\kappa(\cdot,\cdot)$ and $L_{j_{2}}$ is applied along the second argument, and then we evaluate at ${\bf x}_{m_{1}}$ and ${\bf x}_{m_{2}}$. For example, suppose we have $L_{1}(u)=\partial_{x_{1}}u$ and $L_{2}(u)=\partial_{x_{2}}u$, then $L_{1}\circ L_{2}(\kappa)({\bf x},{\bf x}^{\prime})={\partial^{2}\kappa({\bf x}
,{\bf x}^{\prime})}/{\partial x_{1}\partial x^{\prime}_{2}}$ where $x_{1}$ and $x_{2}^{\prime}$ are the first and second elements of the inputs ${\bf x}$ and ${\bf x}^{\prime}$, respectively.

The true solution $u^{*}$ is assumed to reside in ${\mathcal{U}}$, and $\|\cdot\|_{\mathcal{U}}$ is the RKHS norm associated with ${\mathcal{U}}$. From the representation theorem (Owhadi and Scovel,, ), it is straightforward to show that the optimum of (2), denoted by $v_{M}$, has the following form

$$ $\displaystyle v_{M}({\bf x})=c(v_{M}({\bf x}),{\bf z}){\bf C}({\bf z},{\bf z}) ^{-1}{\bf z},$ (3) $$

where $c(v_{M}({\bf x}),{\bf z})$ is the similarity between $v_{M}(x)$ and each element in ${\bf z}$, namely, $\forall z^{j}_{m}\in{\bf z}$,

$$ $c(v_{M}({\bf x}),z^{j}_{m})=L_{j}(\kappa)({\bf x},{\bf x}_{m}),$ $$

where $L_{j}$ is applied along the second argument of $\kappa(\cdot,\cdot)$, and ${\bf x}_{m}$ is the collocation point corresponding to $z^{j}_{m}$. Here ${\bf C}({\bf z},{\bf z})$ is the Gram matrix of ${\bf z}$, constructed from $Q\times Q$ sub-blocks,

$$ $\displaystyle{\bf C}$ $\displaystyle=\left(\begin{array}[]{ccc}{\bf C}_{11}&\ldots&{\bf C}_{1Q}\\ \vdots&\ddots&\vdots\\ {\bf C}_{Q1}&\ldots&{\bf C}_{QQ}\\ \end{array}\right),$ (7) $$

where each ${\bf C}_{ij}$ is the similarity matrix associated with a pair of linear operators,

$$ $\displaystyle{\bf C}_{ij}$ $\displaystyle=c({\bf h}_{i},{\bf h}_{j})=L_{i}\circ L_{j}(\kappa)({\bf h}_{i}, {\bf h}_{j}),$ (8) $$

where $1\leq i,j\leq Q$, $L_{i}$ and $L_{j}$ are applied to the first and second arguments of $\kappa$, respectively,
${\bf h}_{i}=\{z_{m}^{i}\}_{m}$ and ${\bf h}_{j}=\{z_{m}^{j}\}_{m}$. Therefore, the size of the Gram matrix ${\bf C}$ is $\left(Q_{\Omega}M_{\Omega}+(Q-Q_{\Omega})(M-M_{\Omega})\right)\times\left(Q_{
\Omega}M_{\Omega}+(Q-Q_{\Omega})(M-M_{\Omega})\right)$.

The method () can be explained from a probabilistic perspective. That is, we assign a GP prior over $u$, and given a sufficiently smooth kernel $\kappa$, all the linear operators over $u$, namely, $L_{j}(u)$ also follow a GP prior, and their projection on the collocation points ${\mathcal{M}}$, namely, ${\bf z}$, follow a multi-variate Gaussian prior distribution, $p({\bf z})\sim\mathcal{N}({\bf z}|{\bf 0},{\bf C})$. Softening the outer constraints in the aforementioned nested optimization by maximizing a likelihood, this method essentially seeks for an MAP estimation of ${\bf z}$, and the prediction (3) is the posterior mean conditioned on ${\bf z}$.

## 3 Our Framework

Despite the success of (), it requires computation of a Gram (covariance) matrix (see (7)) with dimensions typically exceeding the number of collocation points.
This can exacerbate computational challenges, particularly when addressing complex PDEs that demand a considerable number of collocation points (Cho et al.,, ; Florido et al.,, ). 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.

We therefore propose an alternative kernel learning framework for nonlinear PDE solving, which simplifies the Gram matrix construction and computation. Specifically, we are inspired by the standard kernel/GP regression. Suppose the solution values at the collocation points are known, denoted as ${\bf u}^{*}_{\mathcal{M}}=\left(u^{*}({\bf x}_{1}),\ldots,u^{*}({\bf x}_{M})
\right)^{\top}$. The optimal solution estimate within the framework of standard kernel regression takes the interpolation form:
$t({\bf x})=\kappa(\mathbf{x},\mathcal{M})\mathbf{K}_{MM}^{-1}\mathbf{u}^{*}_{
\mathcal{M}}$,
where $\mathbf{K}_{MM}=\kappa(\mathcal{M},\mathcal{M})$ denotes the kernel matrix computed on the collocation points (of size $M\times M$). This form is derived by minimizing the RKHS norm while aligning $u^{*}$ at $\mathcal{M}$. In GP regression, $t(\mathbf{x})$ serves as the mean function of the posterior process.

In the broader context of PDE solving, as depicted in (1), one often lacks knowledge of the solution values at arbitrary collocation points. Therefore, we regard them as unknown, free variables denoted by $\bm{\eta}$. We model the solution estimate as

$$ $\displaystyle u({\bf x};\bm{\eta})=\kappa({\bf x},{\mathcal{M}}){\bf K}_{MM}^{ -1}\bm{\eta}.$ (9) $$

We then apply each linear operator $L_{j}$ in the PDE over $u({\bf x};\bm{\eta})$ to approximate $L_{j}(u^{*})$. Following this way, we do not need to explicitly estimate the values of $L_{j}(u^{*})$ at the collocation points (namely $z^{j}_{m}$ in (2)). Correspondingly, the Gram matrix ${\bf K}_{MM}$ is substantially smaller (of size $M\times M$) and it is more convenient to compute — there is no need to enumerate pairs of linear operators and apply them to the kernel function to compute different sub-blocks.

The learning is carried out by addressing the following constrained optimization problem:

$$ $\displaystyle\begin{cases}\underset{u\in{\mathcal{U}}}{\text{minimize}}\;\;\|u \|_{\mathcal{U}}\\ \text{s.t. }\frac{1}{M_{\Omega}}\sum_{m=1}^{M_{\Omega}}\left({\mathcal{P}}(u)( {\bf x}_{m})-f({\bf x}_{m})\right)^{2}\\ +\frac{1}{M-M_{\Omega}}\sum_{m=M_{\Omega}+1}^{M}\left({\mathcal{B}}(u)({\bf x} _{m})-g({\bf x}_{m})\right)^{2}\leq\epsilon,\\ \text{$u$ takes the kernel interpolation form }\eqref{eq:u-form}.\end{cases}$ (10) $$

where $\epsilon\geq 0$ is a given relaxation parameter. Note that since our formulation (9) uses a reduced model space as opposed to (), we introduce $\epsilon$ to enable feasibility of the optimization and to establish the convergence; see our convergence analysis in Appendix Section A.
In practical scenarios, addressing (10) directly can be cumbersome. We may opt to optimize an unconstrained objective with soft regularization instead,

$$ $\displaystyle\underset{\bm{\eta}}{\text{minimize}}\;\;$ $\displaystyle\mathcal{L}(u(x;\bm{\eta});\alpha,\beta):=\|u\|_{\mathcal{U}}^{2} +\alpha\cdot\left[\frac{1}{M_{\Omega}}\sum\nolimits_{m=1}^{M_{\Omega}}({ \mathcal{P}}(u)({\bf x}_{m})-f({\bf x}_{m}))^{2}-\epsilon/2\right]$ $\displaystyle+\beta\cdot\left[\frac{1}{M-M_{\Omega}}\sum\nolimits_{m={M_{ \Omega}+1}}^{M}({\mathcal{B}}(u)({\bf x}_{m})-g({\bf x}_{m}))^{2}-\epsilon/2 \right],$ (11) $$

where $\alpha,\beta>0$ are the regularization strengths, and $\epsilon$ can be simply set to zero.

Efficient Computation. In scenarios where PDEs are complex and challenging, capturing the solution details might necessitate employing a vast array of collocation
points. Since our framework uses the standard kernel matrix to construct the solution estimate (as illustrated in (9)), a wide range of existing kernel approximation and computation methods (Quinonero-Candela and Rasmussen,, ; Rahimi and Recht,, ; Farahat et al.,, ; Lindgren et al.,, ) can be readily employed to accelerate computation involving $\mathbf{K}_{MM}^{-1}$. This facilitates the reduction of computational costs and enables scalability to accommodate massive collocation points.

As an instance, we propose to induce a Kronecker product structure to accelerate computation and scale to a large number of collocation points. Specifically, we place the collocation points on a grid, namely, ${\mathcal{M}}={\bf s}^{1}\times\ldots\times{\bf s}^{d}$, where each ${\bf s}^{k}$ includes a collection of locations at input dimension $k$, i.e., ${\bf s}^{k}=(s^{k}_{1},\ldots,s^{k}_{m_{k}})^{\top}\in\mathbb{R}^{m_{k}}$. These locations can be regular-spaced or randomly sampled. Accordingly, ${\mathcal{M}}$ is an $d$-dimensional array of size $m_{1}\times\ldots\times m_{d}$. Next, we employ a product kernel that is decomposed as along the input dimensions, $\kappa({\bf x},{\bf x}^{\prime})=\prod_{j=1}^{d}\kappa_{j}(x_{j},x^{\prime}_{j})$, where each $\kappa_{j}$ is a kernel function of two scalar variables at input dimension $j$. As a result, the kernel matrix on the collocation points ${\mathcal{M}}$ becomes a Kronecker product,

$$ ${\bf K}_{MM}={\bf K}_{1}\otimes\ldots\otimes{\bf K}_{d},$ $$

where each ${\bf K}_{j}=\kappa_{j}({\bf s}^{j},{\bf s}^{j})$ is a local kernel matrix for dimension $j$ ($1\leq j\leq d$), of size $m_{j}\times m_{j}$. We then leverage the Kronecker product properties to efficiently compute the solution estimate (9),

$$ $\displaystyle u({\bf x};\bm{\eta})=\left[\kappa_{1}(x_{1},{\bf s}^{1})\otimes \ldots\otimes\kappa_{d}({\bf x}_{d},{\bf s}^{d})\right]\cdot\left[{\bf K}_{1} \otimes\ldots\otimes{\bf K}_{d}\right]^{-1}\bm{\eta}$ $\displaystyle=\left[\kappa_{1}(x_{1},{\bf s}^{1}){\bf K}_{1}^{-1}\otimes\ldots \otimes\kappa_{d}(x_{d},{\bf s}^{d}){\bf K}_{d}^{-1}\right]\bm{\eta}$ $\displaystyle={\mathcal{A}}\times_{1}\left[\kappa_{1}(x_{1},{\bf s}^{1}){\bf K }_{1}^{-1}\right]\times_{2}\ldots\times_{d}\left[\kappa_{d}(x_{d},{\bf s}^{d}) {\bf K}_{d}^{-1}\right],$ (12) $$

where ${\mathcal{A}}$ is the tensor view of $\bm{\eta}$, namely reshaping $\bm{\eta}$ as a $m_{1}\times\ldots\times m_{d}$ array, and $\times_{k}$ is the mode-$k$ tensor-matrix multiplication (Kolda,, ). In this way, we avoid explicitly computing the full kernel matrix ${\bf K}_{MM}$ and its inverse. We only need to invert each local kernel matrix ${\bf K}_{j}$, and hence the cost is substantially reduced. For example, considering a $100\times 100\times 100$ grid, the full kernel matrix is $10^{6}\times 10^{6}$, rending it computationally prohibitive and impracticable for most hardware. By using (12), we only need to invert three $100\times 100$ local kernel matrices, which is cheap and fast. Note that, since the kernel is decomposed across individual dimensions, taking derivatives over the solution estimate $u$ will maintain the structure, e.g., $\partial_{x_{1}x_{d}}u({\bf x};\bm{\eta})={\mathcal{A}}\times_{1}\left[
\partial_{x_{1}}\kappa_{1}(x_{1},{\bf s}^{1}){\bf K}_{1}^{-1}\right]\times_{2}
\ldots\times_{d}\left[\partial_{x_{d}}\kappa_{d}(x_{d},{\bf s}^{d}){\bf K}_{d}
^{-1}\right]$.

We then leverage the structure (12) to efficiently minimize the objective (10) or (11). The computation of each ${\mathcal{P}}(u)({\bf x}_{m})$ and ${\mathcal{B}}(u)({\bf x}_{m})$ is a straightforward application of the operators ${\mathcal{P}}$ and ${\mathcal{B}}$ to (12) and then evaluate them at the collocation points. This can be done by automatic differential libraries, such as JAX (Frostig et al.,, ). The RKHS norm in (10) and (11) can be efficiently computed by

$$ $\displaystyle\|{\bf u}\|_{\mathcal{U}}^{2}$ $\displaystyle=\bm{\eta}^{\top}{\bf K}_{MM}^{-1}\bm{\eta}=\bm{\eta}^{\top}\left [{\bf K}_{1}\otimes\ldots\otimes{\bf K}_{d}\right]^{-1}\bm{\eta}=\bm{\eta}^{ \top}{\rm vec}({\mathcal{A}}\times_{1}{\bf K}_{1}^{-1}\times_{2}\ldots\times_{ d}{\bf K}_{d}^{-1}).$ (13) $$

We can apply any gradient-based optimization algorithm..

## 4 Convergence Analysis

We now show the convergence of our method. We inherit the road-map of (Batlle et al.,, ) and maintain the same assumption about PDE stability and the regularity of the domain and boundary (Batlle et al.,, , Assumption 3.7)(^1^11Note 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.,, ).), with a slight modification.

###### Assumption 4.1 .

The following conditions hold:

- •
(C1) (Regularity of the domain and its boundary) $\Omega\subset\mathbb{R}^{d}$ with $d>1$ is a compact set and $\partial\Omega$ is a smooth connected Riemannian manifold of dimension $d-1$ endowed with a geodesic distance $\rho_{\partial\omega}$.
- •
(C2) (Stability of the PDE) $\exists\gamma>0$ and $\exists k,t\in\mathbb{N}$ with $k>d/2$ and $t>(d-1)/2$,
and $\exists s,l\in\mathbb{R}$ such that for any $r>0$, it holds that $\forall u_{1},u_{2}\in B_{r}(H^{l}(\Omega))$,
$\displaystyle\|u_{1}-u_{2}\|_{H^{l}(\Omega)}\leq C\left(\|{\mathcal{P}}(u_{1})
-{\mathcal{P}}(u_{2})\|_{H^{{0}}(\Omega)}+\|{\mathcal{B}}(u_{1})-{\mathcal{B}}
(u_{2})\|_{H^{{0}}(\partial\Omega)}\right),$
(14)
and $\forall u_{1},u_{2}\in B_{r}(H^{s}(\Omega))$,
$\displaystyle\|{\mathcal{P}}(u_{1})-{\mathcal{P}}(u_{2})\|_{H^{k}(\Omega)}+\|{
\mathcal{B}}(u_{1})-{\mathcal{B}}(u_{2})\|_{H^{t}(\partial\Omega)}\leq C\|u_{1
}-u_{2}\|_{H^{s}(\Omega)},$
(15)
where $C=C(r)>0$ is a constant independent of $u_{1}$ and $u_{2}$, $B(r)$ is an open ball with radius $r$, $H^{j}=W^{j,2}$ is a Sobolev space where each element and its weak derivatives up to the order of $j$ have a finite $L^{2}$ norm.
- •
(C3) The RKHS ${\mathcal{U}}$ is continuously embedded in $H^{s+\tau}(\Omega)$ where $\tau>0$.

###### Lemma 4.2 .

Let $u^{*}$ denote the unique strong solution of (1).
Suppose Assumption 4.1 is satisfied, and a set of collocation points ${\mathcal{M}}\subset\overline{\Omega}$ is given, where ${\mathcal{M}}_{\Omega}\subset{\mathcal{M}}$ denotes the collocation points in the interior of $\Omega$ and ${\mathcal{M}}_{\partial\Omega}\subset{\mathcal{M}}$ the collocation points on the boundary $\partial\Omega$. Define the fill-distances

$$ $\displaystyle h_{\Omega}:=\underset{{\bf x}\in\Omega}{\sup}\;\;\underset{{\bf x }^{\prime}\in{\mathcal{M}}_{\Omega}}{\inf}|{\bf x}-{\bf x}^{\prime}|,\;\;\;h_{ \partial\Omega}:=\underset{{\bf x}\in\partial\Omega}{\sup}\;\;\underset{{\bf x }^{\prime}\in{\mathcal{M}}_{\partial\Omega}}{\inf}\rho_{\partial\Omega}({\bf x },{\bf x}^{\prime}),$ (16) $$

where $|\cdot|$ is the Euclidean distance, and $\rho_{\partial\Omega}$ is a geodesic distance defined on $\partial\Omega$. Set $h=\text{max}(h_{\Omega},h_{\partial\Omega})$.
There is always a minimizer of (10) with the set of collocation points ${\mathcal{M}}$ and $\epsilon=C_{0}h^{2\tau}$ where $C_{0}>0$ is a sufficiently large constant independent of $h$.
Let $u^{\dagger}$ denote such a minimizer.
When ${h}$ is sufficiently small, at least $h\leq C_{1}M^{-\frac{1}{d}}$ where $C_{1}>0$ is a constant, then

$$ $\displaystyle\|u^{\dagger}-u^{*}\|_{H^{l}(\Omega)}\leq Ch^{\rho}\|u^{*}\|_{ \mathcal{U}},$ (17) $$

where $\rho=\min(k,t,\tau)$, and $C>0$ is independent of $u^{\dagger}$ and $h$.

###### Proposition 4.3 .

Given the set of collocation points ${\mathcal{M}}$ and $\epsilon=C_{0}h^{2\tau}$ where $C_{0}>0$ is a sufficiently large constant, there exists $\alpha_{M},\beta_{M}>0$ such that the minimizer of (11) with $\alpha=\alpha_{M}$ and $\beta=\beta_{M}$ is also the minimizer of (10). That means, with proper choices of the regularization strengths, the minimizer of (11) enjoys the same convergence result as in (17).

We can see the convergence results of our framework are as comparably strong as the results for the method of ; see
(Batlle et al.,, , Theorem 3.8), though the latter employs a richer model space.
We leave the proof in Section A and B of the Appendix.

###### Assumption 4.1 .

###### Lemma 4.2 .

###### Proposition 4.3 .

## 5 Related Work

The prior works of Graepel, (); Raissi et al., () propose Gaussian Processes (GPs) models for solving linear PDEs in the presence of noisy measurements of source terms. Recently, Wang et al., () delve into the rationale and guarantees of using GPs as a prior for PDE solutions. The effectiveness of the product kernel is also justified in terms of sample path properties.

introduced a kernel method capable 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 matrix whose dimension is typically much greater than the number of collocation points. In (Batlle et al.,, ), a systematic theoretical framework is established to analyze the convergence and the convergence rate of the method of . To alleviate the computational challenge for massive collocation points, Chen et al., () adapted the sparse inverse Cholesky factorization (Schafer et al.,, )
to approximate the Gram matrix of (). An alternative approach was proposed by Meng and Yang, (), which adjusted the Nyström method (Jin et al.,, ) 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 sub-block corresponds a pair of linear operators over the kernels, and hence it is complex and relatively inconvenient for implementation. In our work, the solution is approximated by a standard kernel interpolation, 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.

The computational efficiency of Kronecker product structures has been recognized in various works (Saatcci,, ; Xu et al.,, ; Wilson and Nickisch,, ; Izmailov et al.,, ; Zhe et al.,, ). Wilson et al., () highlighted that utilizing a regular (evenly-spaced) grid results in Toeplitz-structured kernel matrices, facilitating $O(n\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 on a grid is natural, making Kronecker products combined with kernels a promising avenue for efficient computation. The recent work (Fang et al.,, ) uses a similar computational method to solve high-frequency and multi-scale PDEs. The major contribution is to introduce a spectral mixture kernel in each dimension to capture the dominant frequencies in the kernel space. This work can be viewed as an instance of our proposed framework. We in addition give a theoretical analysis about the convergence of our framework.
For more extensive discussions on Bayesian learning and PDE problems, readers are referred to (Owhadi,, ).

Our work is also connected to the radial basis function (RBF) method for solving PDEs (Hardy,, ; Kansa,, ; Tolstykh and Shirobokov,, ; Shu et al.,, ; Fornberg et al.,, ; Safdari-Vaighani et al.,, ). The RBF method typically approximates the PDE solution as a linear combination of RBF bases, i.e., kernels, $u({\bf x})\approx\sum_{j}\alpha_{j}\kappa({\bf x}-{\bf x}_{j})$ where $\{{\bf x}_{j}$} are collocation points. To estimate the coefficients $\{\alpha_{j}\}$, the RBF method typically converts the PDE into a linear system ${\bf A}\bm{\alpha}={\bf b}$ and accordingly solves $\bm{\alpha}=(\alpha_{1},\alpha_{2},\ldots)^{\top}$. Our method differs mainly in two folds. First, we use kernel regression form (9) to approximate the solution, and do not explicitly estimate the coefficients $\bm{\alpha}$. In other words, we have implicitly $\bm{\alpha}={\bf K}_{MM}^{-1}\bm{\eta}$. We instead estimate $\bm{\eta}$ — the solution values at the collocation points. Second, we do not convert the PDE solving into solving a linear system; instead, we convert it into solving an optimization problem (10) or (11). While we can also use the RBF formulation to approximate the solution in our framework, we found that, optimizing $\bm{\alpha}$ (instead of $\bm{\eta}$) will severely degrade the performance, which might be due to the complexity of the optimization procedure.

## 6 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.,, ; ).

The Burgers’ Equation. We first tested with a viscous Burgers’ equation,

$$ $\displaystyle u_{t}+uu_{x}-\nu u_{xx}=0,\;\;\forall(x,t)\in(-1,1)\times(0,1],$ $\displaystyle u(x,0)=-\sin(\pi x),\;\;u(-1,t)=u(1,t)=0.$ (18) $$

The solution is computed from the Cole–Hopf transformation with numerical quadrature (). We considered two cases: $\nu=0.02$, and $\nu=0.001$.

Nonlinear elliptic PDE. We next tested with the instance of nonlinear elliptic PDE used in (),

$$ $\displaystyle-\Delta u({\bf x})+u^{3}({\bf x})$ $\displaystyle=f({\bf x}),\;\;\;\forall{\bf x}\in\Omega,$ $\displaystyle u({\bf x})$ $\displaystyle=0,\;\;\;\forall{\bf x}\in\partial\Omega,$ (19) $$

where $\Omega=[0,1]^{2}$, the solution is crafted as $u({\bf x})=\sin(\pi x_{1})\sin(\pi x_{2})+4\sin(4\pi x_{1})\sin(4\pi x_{2})$, and $f({\bf x})$ is correspondingly computed via the equation.

Eikonal PDE. Third, we tested with a regularized Eikonal equation as used in (),

$$ $\displaystyle|\nabla u({\bf x})|^{2}=f({\bf x})^{2}+\epsilon\Delta u({\bf x}), \;\;\forall{\bf x}\in\Omega,$ $\displaystyle u({\bf x})=0,\;\;\forall{\bf x}\in\partial\Omega,$ (20) $$

where $\Omega=[0,1]^{2}$, $f({\bf x})=1$, and $\epsilon=0.1$.
The solution is computed from a highly-resolved finite difference solver as provided by ().

Allen-Cahn Equation. Fourth, we considered a 2D stationary Allen-Cahn equation with a source function and Dirichlet boundary conditions.

$$ $\displaystyle u_{xx}+u_{yy}+\gamma(u^{m}-u)$ $\displaystyle=f(x,y),\;\;(x,y)\in[0,1]\times[0,1],$ (21) $$

where $\gamma=1$ and $m=3$. We crafted the solution in the form $u=\sin(2\pi ax_{1})\cos(2\pi ax_{2})+\sin(2\pi x_{1})\cos(2\pi x_{2})$, and $f$ is computed through the equation. We tested with $a=15$ and $a=20$.

Method and Settings. We implemented our method with JAX (Frostig et al.,, ). We denote our method as SKS (Simple Kernel-based Solver)
We compared with () 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 implementation from the original authors(^2^22[https://github.com/yifanc96/NonLinPDEs-GPsolver](https://github.com/yifanc96/NonLinPDEs-GPsolver)). In addition, we compared with physics-informed neural network (PINN) (Raissi et al.,, ), a mainstream machine learning PDE solver. The PINN is implemented with PyTorch (Paszke et al.,, ).
For SKS, we minimize (11) (with $\epsilon=0$), and used ADAM optimization with learning rate $10^{-3}$. The maximum number of epochs was set to 1M. For DAKS, we used the relaxed Gauss-Newton optimization propsoed in the original paper across all the experiments. The PINN was first trained by 10K ADAM epochs with learning rate $10^{-3}$ 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^{-9}$.
To identify the architecture for the PINN, we varied the number of layers from {2, 3, 5, 8, 10}, and the the width of each layer from $\{10,20,30,\ldots,100\}$. We used tanh as the activation function.
For DAKS and SKS, we used Square Exponential (SE) kernel with different length-scales across the input dimensions. We selected the nugget term from {5E-5, 1E-5, 5E-6, 1E-6, $\ldots$, 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 a grid search, from $[0.1,0.2]^{2}$ for the nonlinear elliptic and Eikonal PDEs, $[0.05,0.01]^{2}$ for Allen-Cahn equation, and $[0.003,0.05]\times[0.02,0.3]$ for Burgers’ equation.
We reported 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.

### 6.1 Solution Accuracy

Simpler Cases. We first tested all the methods on less challenging benchmarks, for which a small number of collocation points is sufficient. Specifically, we tested with Burgers’ equation with viscosity $\nu=0.02$, the nonlinear elliptic PDE, and Eikonal PDE. These are the same test cases employed in (). Following  (), we varied the number of collocation points from
{600, 1200, 2400, 4800} for Burgers’ equation, and
{300, 600, 1200, 2400} for nonlinear elliptic PDE and Eikonal PDE. In (), the collocation points are randomly sampled, and hence we used the average $L^{2}$ of DAKS from ten runs on different sets of randomly sampled collocation points for comparison. For SKS, we used a regularly-spaced, 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 ran PINN on the same set of grid points used by SKS. The $L^{2}$ error is reported in Table 1. It can be seen that in most cases, SKS achieves smaller solution error than DAKS, with the exception for Burgers(^3^33However, if we use a $96\times 50$ grid (still including 4800 collocation points), SKS gives $L^{2}$ error 7.54E-05, which can surpass DAKS. See the ablation study in Appendix Section C.2 for more details.) and nonlinear elliptic PDE using 4800 and 2400 collocation points, respectively. This might be because DAKS needs to explicitly estimate more variables, including all kinds of linear operators (e.g., derivatives) over the solution at the collocation points, which increases the optimization 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 (18) with viscosity $\nu=0.02$.**
| Method | 600 ($25\times 25$) | 1200 ($35\times 35$) | 2400 ($49\times 49$) | 4800 ($70\times 70$) |
| --- | --- | --- | --- | --- |
| DAKS | 1.75E-02 | 7.90E-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 |

Difficult Cases. Next, we tested with more challenging cases, for which massive collocation points are necessary. These cases include Burgers’ equation with $\nu=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 a square-shaped grid. We provide a 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^{2}$ error at level $10^{-4}$ for solving the 2D Allen-Cahn equation, while the other methods still struggle at $10^{-1}$ error level or even bigger. Together this indicates that a much larger number of collocation points is needed.

We then ran SKS and PINN with greatly increased collocation points, i.e., dense grids, varying from 6400 to 120K. In such scenarios, running DAKS becomes extremely costly or even infeasible(^4^44For 6400 collocation points, the size of the Gram matrix of DAKS is $19200\times 19200$ for Burgers’ and 2D Allen-Cahn, because there are three linear operators in each equation.). We therefore only report the results of SKS and PINN, as shown in Table 3. 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}$ level $L^{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 (21)) are difficult to be captured by neural networks, due to their known “spectral bias” (Rahaman et al.,, ). Thanks to our model design (3), we can induce a Kronecker product structure in the Gram matrix to scale to massive collocation points, without the need for designing complex approximations.

**(a) The Burgers’ equation (18) with viscosity $\nu=0.001$.**
| Method | 600 ($42\times 14$) | 1200 ($60\times 20$) | 2400 ($84\times 28$) | 4800 ($120\times 40$) |
| --- | --- | --- | --- | --- |
| DAKS | 3.87E-01 | 3.12E-01 | 3.60E-01 | 2.37E-01 |
| PINN | 1.27E-01 | 2.59E-01 | 3.18E-01 | 2.65E-01 |
| SKS | 1.34E-01 | 1.11E-01 | 8.04E-02 | 1.89E-02 |

**(a) The Burgers’ equation (18) with viscosity $\nu=0.001$.**
| Method | 43200 ($360\times 120$) | 67500 ($450\times 150$) | 97200 ($540\times 180$) | 120000 ($600\times 200$) |
| --- | --- | --- | --- | --- |
| PINN | 4.05E-03 | 6.01E-03 | 3.94E-03 | 4.13E-03 |
| SKS | 3.90E-03 | 3.50E-03 | 2.60E-03 | 2.28E-03 |

Point-wise Error. For a fine-grained comparison, we showcase the point-wise error of each method in solving Burgers’ ($\nu=0.001$) and 2D Allen-Cahn equations. The results and discussion are given by Appendix Section C.1.

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 $\nu=0.02$ and $\nu=0.001$. We leave the details in Appendix Section C.2.

**(a) Nonlinear Elliptic PDE (19).**
| Method | $18\times 18$ | $25\times 25$ | $35\times 35$ | $49\times 49$ |
| --- | --- | --- | --- | --- |
| Finite Difference | 3.36E-02 | 1.78E-02 | 9.25E-03 | 4.78E-03 |
| SKS | 1.26E-02 | 6.93E-05 | 6.80E-06 | 1.83E-06 |

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 a centered second-order numerical difference to 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 (19) and the Allen-cahn equation (21), since the ground-truth solutions of these PDEs are known and we can conduct a fair comparison. The results are presented in Table 4. Note that the $L^{2}$ errors of PINN and SKS have already been reported in Table 1(b), 3(b) and 3(c). 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\times 18$ increasing to $25\times 25$, and Allen-Cahn ($a=15$) with the grid $90\times 90$ increasing to $150\times 150$. Since our method is efficient in handling a large number of collocation points (i.e., dense grid), it shows the potential of our method.

Irregular-shaped Domains.
While our efficient computation is performed on grids, our method can be readily applied to irregular-shaped domains by introducing a virtual grid that encompasses such domains. This allows SKS to be used without any modifications. To validate the effectiveness of this strategy, we conducted additional tests, solving the nonlinear elliptic PDE (19) on a circular domain and the Allen-Cahn equation (21) on a triangular domain. In both cases, our method achieved reasonably good accuracy. Detailed results and discussions are provided in Appendix Section C.3.

Running Time.
Finally, we examined the wall-clock runtime of SKS. We analyzed the runtime of each method when solving Burgers’ equation $\nu=0.001$ and the Allen-Cahn equation ($a=15$) with varying numbers of collocation points. We ran the experiment on a Linux workstation equipped with an Intel(R) Xeon(R) Platinum 8360H Processor with 24GB memory.
The results are presented in Appendix Table 9. 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 collocation points than DAKS. Overall, the runtime of SKS is significantly less than that of PINN.

## 7 Conclusion

We have proposed a new kernel method for nonlinear PDE solving. We use a standard kernel interpolation to model the solution estimate, which allows more convenient implementation and efficient computation. Our method can easily scale to massive collocation points, which are necessary for challenging PDEs. The performance on a series of benchmarks is encouraging. In the future, we plan to develop more efficient optimization, e.g., Gaussian-Newton, to further accelerate our method.

## 8 Acknowledgements

This work has been supported by MURI AFOSR grant FA9550-20-1-0358, NSF CAREER Award IIS-2046295 and NSF 2311685. We thank Bamdad Hosseini for valuable discussions and suggestions.

## Appendix

## Appendix A Proof of Lemma 4.2

Since for any collocation point ${\bf x}_{m}$, ${\mathcal{P}}(u^{*})({\bf x}_{m})=f({\bf x}_{m})$ when ${\bf x}_{m}\in\Omega$ and ${\mathcal{B}}(u^{*})({\bf x}_{m})=g({\bf x}_{m})$ when ${\bf x}_{m}\in\partial\Omega$, we can re-write (10) as

$$ $\displaystyle\begin{cases}\underset{u\in{\mathcal{U}}}{\text{minimize}}\;\;\|u \|_{\mathcal{U}}\\ \text{s.t. }\frac{1}{M_{\Omega}}\sum_{m=1}^{M_{\Omega}}\left({\mathcal{P}}(u)( {\bf x}_{m})-{\mathcal{P}}(u^{*})({\bf x}_{m})\right)^{2}\\ +\frac{1}{M-M_{\Omega}}\sum_{m=M_{\Omega}+1}^{M}\left({\mathcal{B}}(u)({\bf x} _{m})-{\mathcal{B}}(u^{*})({\bf x}_{m})\right)^{2}\leq\epsilon,\\ \text{$u$ takes the kernel interpolation form }\eqref{eq:u-form}.\end{cases}$ (22) $$

Step 1. In the first step, we show that for all $\epsilon$ above some threshold (depending on $h$), there exists a minimizer $u^{\dagger}$ for (22), and we would like also to bound the RKHS norm of $u^{\dagger}$, namely $\|u^{\dagger}\|_{\mathcal{U}}$. To this end, we utilize an intermediate optimization problem,

$$ $\displaystyle\begin{cases}\underset{u\in{\mathcal{U}}}{\text{minimize}}\;\;\|u \|_{\mathcal{U}}\\ \text{s.t. }u({\bf x}_{m})=u^{*}({\bf x}_{m}),\;\;1\leq m\leq M.\end{cases}$ (23) $$

Denote the minimizer of (23) by $u^{*}_{M}$.
This is a standard kernel regression problem. According to the representation theorem, $u^{*}_{M}$ takes the kernel interpolation form (9), and $\|u^{*}_{M}\|_{\mathcal{U}}\leq\|u^{*}\|_{\mathcal{U}}$.

Since $u^{*}_{M}-u^{*}$ is zero at all the collocations points, according to the sampling inequality (see Proposition A.1 of (Batlle et al.,, )), when the fill-distance $h$ is sufficiently small (note $h_{\Omega}\leq h$),

$$ $\displaystyle\|u^{*}_{M}-u^{*}\|_{H^{s}(\Omega)}\lesssim h^{\tau}\|u^{*}_{M}-u ^{*}\|_{H^{s+\tau}(\Omega)}$ (24) $$

where $\lesssim$ means the inequality holds with a positive constant factor multiplied by the right-hand side, and the constant is independent of the terms on both sides. Combining with (C2) of Assumption 4.1, we can obtain

$$ $\displaystyle\|{\mathcal{P}}(u^{*}_{M})-{\mathcal{P}}(u^{*})\|_{H^{k}(\Omega)} +\|{\mathcal{B}}(u^{*}_{M})-{\mathcal{B}}(u^{*})\|_{H^{t}(\partial\Omega)} \lesssim h^{\tau}\|u^{*}_{M}-u^{*}\|_{H^{s+\tau}(\Omega)}.$ (25) $$

Since ${\mathcal{U}}$ is continuously embedded in $H^{s+\tau}(\Omega)$ — C(3) of Assumption 4.1, we have

$$ $\displaystyle\|u_{M}-u^{*}\|_{H^{s+\tau}({\mathcal{X}})}\lesssim\|u_{M}-u^{*} \|_{\mathcal{U}}.$ (26) $$

Combining (25), (26) and the fact $\|u^{*}_{M}\|_{\mathcal{U}}\leq\|u^{*}\|_{\mathcal{U}}$, we have

$$ $\displaystyle\|{\mathcal{P}}(u_{M}^{*})-{\mathcal{P}}(u^{*})\|_{H^{k}(\Omega)} +\|{\mathcal{B}}(u_{M}^{*})-{\mathcal{B}}(u^{*})\|_{H^{t}(\partial\Omega)} \lesssim h^{\tau}\|u^{*}\|_{\mathcal{U}}.$ (27) $$

According to (C2) of Assumption 4.1, since $k>\frac{d}{2}$ and $t>\frac{d-1}{2}$ , according to Sobolev embedding theorem (Adams and Fournier,, , Theorem 4.12), both $H^{k}(\Omega)$ and $H^{t}(\partial\Omega)$ are continuously embedded into $C^{0}(\Omega)$ and $C^{0}(\partial\Omega)$, respectively.
Therefore,

$$ $\displaystyle\|{\mathcal{P}}(u^{*}_{M})-{\mathcal{P}}(u^{*})\|_{C^{0}(\Omega)} \lesssim\|{\mathcal{P}}(u^{*}_{M})-{\mathcal{P}}(u^{*})\|_{H^{k}(\Omega)},$ $\displaystyle\|{\mathcal{B}}(u_{M}^{*})-{\mathcal{B}}(u^{*})\|_{C^{0}(\partial \Omega)}\lesssim\|{\mathcal{B}}(u_{M}^{*})-{\mathcal{B}}(u^{*})\|_{H^{t}( \partial\Omega)}.$ (28) $$

At any collocation point, we obviously have

$$ $\displaystyle\left({\mathcal{P}}(u^{*}_{M})-{\mathcal{P}}(u^{*})\right)^{2} \leq\|{\mathcal{P}}(u^{*}_{M})-{\mathcal{P}}(u^{*})\|^{2}_{C^{0}(\Omega)},$ $\displaystyle\left({\mathcal{B}}(u^{*}_{M})({\bf x}_{m})-{\mathcal{B}}(u^{*})( {\bf x}_{m})\right)^{2}\leq\|{\mathcal{B}}(u^{*}_{M})-{\mathcal{B}}(u^{*})\|^{ 2}_{C^{0}(\partial\Omega)}.$ (29) $$

Combining (27), (28) and (29), we can obtain that

$$ $\displaystyle\frac{1}{M_{\Omega}}\sum_{m=1}^{M_{\Omega}}\left({\mathcal{P}}(u^ {*}_{M})({\bf x}_{m})-{\mathcal{P}}(u^{*})({\bf x}_{m})\right)^{2}$ $\displaystyle+\frac{1}{M-M_{\Omega}}\sum_{m=M_{\Omega}+1}^{M}\left({\mathcal{B }}(u^{*}_{M})({\bf x}_{m})-{\mathcal{B}}(u^{*})({\bf x}_{m})\right)^{2}\leq Ch ^{2\tau}\|u^{*}\|^{2}_{\mathcal{U}},$ (30) $$

where $C>0$ is a constant independent of $h$ and other terms in the inequality.

The result (30) means that given the collocation points ${\mathcal{M}}$ and $\epsilon=Ch^{2\tau}\|u^{*}\|^{2}_{\mathcal{U}}$, the feasible region of the optimization problem (22) is nonempty and at least includes $u^{*}_{M}$. Therefore, the minimizer of (22) must exist and satisfy

$$ $\displaystyle\|u^{\dagger}\|_{\mathcal{U}}\leq\|u^{*}_{M}\|_{\mathcal{U}}\leq \|u^{*}\|_{\mathcal{U}}.$ (31) $$

Step 2. Next, we analyzed the error of ${\mathcal{P}}(u^{\dagger})$ and ${\mathcal{B}}(u^{\dagger})$. For notation convenience, we define two error functions,

$$ $\displaystyle\xi_{P}({\bf x})$ $\displaystyle={\mathcal{P}}(u^{\dagger})({\bf x})-{\mathcal{P}}(u^{*})({\bf x} ),\;\;{\bf x}\in\Omega,$ $\displaystyle\xi_{B}({\bf x})$ $\displaystyle={\mathcal{B}}(u^{\dagger})({\bf x})-{\mathcal{B}}(u^{*})({\bf x} ),\;\;{\bf x}\in\partial\Omega.$ (32) $$

We would like to bound the $L^{2}$ norm of the error functions, namely, $\|\xi_{P}\|_{H^{0}(\Omega)}$ and $\|\xi_{B}\|_{H^{0}(\partial\Omega)}$. We first consider the case for $\xi_{P}$. The idea is to decompose $\Omega$ into $M_{\Omega}$ regular non-overlapping regions, ${\mathcal{T}}_{1}\cup\ldots\cup{\mathcal{T}}_{M_{\Omega}}=\Omega$, such that each region ${\mathcal{T}}_{i}$ only includes one collocation point ${\bf x}_{i}$, and its filled-distance $h_{i}\lesssim h$ ($1\leq i\leq M_{\Omega}$). We therefore can decompose the squared $L^{2}$ norm as

$$ $\displaystyle\|\xi_{P}\|^{2}_{H^{0}(\Omega)}=\sum_{i=1}^{M_{\Omega}}\int_{{ \mathcal{T}}_{i}}\xi_{P}({\bf x})^{2}{\rm d}{\bf x}=\sum_{i=1}^{M_{\Omega}}\| \xi_{P}\|_{H^{0}({\mathcal{T}}_{i})}^{2}.$ (33) $$

Since according to the mean inequality,

$$ $\xi_{P}({\bf x})^{2}=\left(\xi_{P}({\bf x})-\xi_{P}({\bf x}_{i})+\xi_{P}({\bf x }_{i})\right)^{2}\leq 2\left(\xi_{P}({\bf x})-\xi_{P}({\bf x}_{i})\right)^{2}+ 2\xi_{P}({\bf x}_{i})^{2},$ $$

we immediately obtain

$$ $\displaystyle\|\xi_{P}\|_{H^{0}({\mathcal{T}}_{i})}^{2}\lesssim\|\xi_{P}-\xi_{ P}({\bf x}_{i})\|_{H^{0}({\mathcal{T}}_{i})}^{2}+\lambda({\mathcal{T}}_{i})\xi _{P}({\bf x}_{i})^{2},$ (34) $$

where $\lambda({\mathcal{T}}_{i})$ is the volume of ${\mathcal{T}}_{i}$.

Since the function $\xi_{P}-\xi_{P}({\bf x}_{i})$ takes zero at ${\bf x}_{i}$, we can apply the sampling inequality again. That is, when the fill-distance $h_{i}$ is sufficiently small, we have

$$ $\displaystyle\|\xi_{P}-\xi_{P}({\bf x}_{i})\|_{H^{0}({\mathcal{T}}_{i})} \lesssim h_{i}^{k}\|\xi_{P}-\xi_{P}({\bf x}_{i})\|_{H^{k}({\mathcal{T}}_{i})}.$ (35) $$

Since $h_{i}\lesssim h$, when $h$ is sufficiently small, we further have

$$ $\displaystyle\|\xi_{P}-\xi_{P}({\bf x}_{i})\|_{H^{0}({\mathcal{T}}_{i})} \lesssim h^{k}\|\xi_{P}-\xi_{P}({\bf x}_{i})\|_{H^{k}({\mathcal{T}}_{i})},$ (36) $$

and then using the mean inequality,

$$ $\displaystyle\|\xi_{P}-\xi_{P}({\bf x}_{i})\|^{2}_{H^{0}({\mathcal{T}}_{i})} \lesssim h^{2k}\left(\|\xi_{P}\|^{2}_{H^{k}({\mathcal{T}}_{i})}+\|\xi_{P}({\bf x }_{i})\|^{2}_{H^{k}({\mathcal{T}}_{i})}\right)=h^{2k}\left(\|\xi_{P}\|^{2}_{H^ {k}({\mathcal{T}}_{i})}+\xi_{P}({\bf x}_{i})^{2}\right).$ (37) $$

Since $\lambda({\mathcal{T}}_{i})\lesssim h^{d}$, combining (33), (34) and (37), we can obtain

$$ $\displaystyle\|\xi_{P}\|^{2}_{H^{0}(\Omega)}$ $\displaystyle\lesssim h^{2k}\sum_{i}\|\xi_{P}\|^{2}_{H^{k}({\mathcal{T}}_{i})} +(h^{d}+h^{2k})\sum_{i}\xi_{P}({\bf x}_{i})^{2}$ $\displaystyle\lesssim h^{2k}\|\xi_{P}\|^{2}_{H^{k}(\Omega)}+(h^{d}+h^{2k}) \cdot M_{\Omega}\cdot\epsilon$ (38) $$

where $\epsilon$ comes from the constraint of (22). To ensure feasibility and to establish convergence, we set $\epsilon=Ch^{2\tau}\|u^{*}\|^{2}_{\mathcal{U}}$ as shown in (30).
When $h\lesssim M^{-\frac{1}{d}}$ and is sufficiently small, we have $(h^{d}+h^{2k})M_{\Omega}\leq(h^{d}+h^{2k})M\leq 1+h^{2k-d}\leq 2$ (since $k>d/2$). Therefore, we can extend the R.H.S of (38) to

$$ $\displaystyle\|\xi_{P}\|^{2}_{H^{0}(\Omega)}\lesssim h^{2k}\|\xi_{P}\|^{2}_{H^ {k}(\Omega)}+h^{2\tau}\|u^{*}\|^{2}_{\mathcal{U}}.$ (39) $$

We can follow a similar approach to show that

$$ $\displaystyle\|\xi_{B}\|^{2}_{H^{0}(\partial\Omega)}\lesssim h^{2t}\|\xi_{B}\| ^{2}_{H^{t}(\partial\Omega)}+h^{2\tau}\|u^{*}\|^{2}_{\mathcal{U}}.$ (40) $$

Combining (39) and (40),

$$ $\displaystyle\left(\|\xi_{P}\|_{H^{0}(\Omega)}+\|\xi_{B}\|_{H^{0}(\partial \Omega)}\right)^{2}\lesssim h^{2\cdot\min(t,k)}\left(\|\xi_{P}\|_{H^{k}(\Omega )}+\|\xi_{B}\|_{H^{t}(\partial\Omega)}\right)^{2}+h^{2\tau}\|u^{*}\|^{2}_{ \mathcal{U}}.$ (41) $$

Since ${\mathcal{U}}\hookrightarrow H^{s+\tau}$ — (C3) of Assumption 4.1, we have ${\mathcal{U}}\hookrightarrow H^{s}$.
Leveraging (31) and (C2) of Assumption 4.1, we immediately obtain

$$ $\displaystyle\|\xi_{P}\|_{H^{k}(\Omega)}+\|\xi_{B}\|_{H^{t}(\partial\Omega)} \lesssim\|u^{\dagger}-u^{*}\|_{H^{s}(\Omega)}\lesssim\|u^{\dagger}-u^{*}\|_{{ \mathcal{U}}}\lesssim\|u^{*}\|_{\mathcal{U}}.$ (42) $$

Combining (41) and (42), and (C1) of Assumption 4.1, we arrive at

$$ $\displaystyle\|u^{\dagger}-u^{*}\|_{H^{l}(\Omega)}\lesssim h^{\rho}\|u^{*}\|_{ \mathcal{U}}$ (43) $$

where $\rho=\min(k,t,\tau)$. When $h\rightarrow 0$, obviously $u^{\dagger}$ converges to $u^{*}$.

## Appendix B Proof of Proposition 4.3

The constraint optimization problem (10) is equivalent to the following mini-max optimization problem,

$$ $\displaystyle\min_{u}\max_{w\geq 0}\|u\|_{\mathcal{U}}+w$ $\displaystyle\Bigg{[}\frac{1}{M_{\Omega}}\sum_{m=1}^{M_{\Omega}}\left({ \mathcal{P}}(u)({\bf x}_{m})-{\mathcal{P}}(u^{*})({\bf x}_{m})\right)^{2}$ $\displaystyle+\frac{1}{M-M_{\Omega}}\sum_{m=M_{\Omega}+1}^{M}\left({\mathcal{B }}(u)({\bf x}_{m})-{\mathcal{B}}(u^{*})({\bf x}_{m})\right)^{2}-\epsilon\Bigg{ ]}.$ (44) $$

Suppose the feasible region is non-empty. Denote the optimum of (44) by $(u^{\dagger},w^{\dagger})$. Then $u^{\dagger}$ is a minimizer of (10). Now if we set $\alpha=\beta=w^{\dagger}$ in (11), and optimizing $\eqref{eq:our-loss}$ will recover the minimizer $u^{\dagger}$.

## Appendix C More Results

### C.1 Point-wise Error

For a fine-grained 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 $\nu=0.001$, and from 600, 2400, and 40K on 2D Allen-Cahn equation with both $a=15$ and $a=20$. The results are shown in Fig. 1, 2 and 3. It can be seen that across all the three PDEs, the solution error of SKS decreases more and more along with the increase of collocation points. Note that for Allen-Cahn with $a=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 error decreases substantially as the number of collocation points grows when solving Burgers’ equation (see Fig. 1), but not obviously on solving Allen-Cahn equation (see Fig. 2 and 3). This is consistent with the global error shown in Table 2. This might be because the quantities of collocation point used are not sufficient to lead to a qualitative boost of DAKS. However, scaling up to much more collocation points, such as 400K, incurs a substantial increase of the computational cost.

Figure: Figure 1: Point-wise solution error for Burgers’ equation (18) with viscosity $\nu=0.001$.
Refer to caption: x1.png

Figure: Figure 2: Point-wise solution error for 2D Allen-Cahn equation (21) with $a=15.0$.
Refer to caption: x2.png

Figure: Figure 3: Point-wise solution error for 2D Allen-Cahn equation (21) with $a=20.0$.
Refer to caption: x3.png

### C.2 Ablation Study on Grid Shape

We investigated how the grid shape can influence the performance of our method. To this end, we tested on Burgers’ equation with $\nu=0.02$ and $\nu=0.001$. For the former case, we fixed the number of collocation points to be 4800 and varied the time resolution from 10 to 80, and the spatial resolution is obtained by dividing 4800 by the time resolution and rounding up to an integer. Similarly, for $\nu=0.001$, we fixed the number of the collocation points to 120K, and varied the time resolution from 100 to 600. We show the $L^{2}$ error of using each grid shape in Table 5. It can be seen that the grid shape does influence the error. In particular, when $\nu$ is small, i.e., $\nu=0.001$, the higher the spatial resolution, the smaller the error. The smallest error is achieved when we use the space-time resolution $1200\times 100$. On the other hand, the time resolution seems to have much less effect on the solution accuracy. This is reasonable, because on Burgers’ equation, a smaller viscosity ($\nu$) increases the sharpness of the shock wave (spatial function). Naturally, the higher the spatial resolution, the more accurate the sharpness can be captured. In summary, we believe that in general the grid shape should be viewed as an influence factor in running our method, which needs to be carefully selected. The appropriate choice may also connect to the intrinsic property of the PDE itself.

**(a) Viscosity $\nu=0.02$, with 4800 collocation points.**
| Grid shape | $60\times 80$ | $69\times 70$ | $96\times 50$ | $160\times 30$ | $480\times 10$ |
| --- | --- | --- | --- | --- | --- |
| $L^{2}$ error | 2.27E-03 | 4.10E-04 | 7.54E-05 | 1.72E-04 | 5.98E-04 |

**Table 6: $L^{2}$ error of solving PDEs on irregularly-shaped domains.**
| $L^{2}$ Error | SKS | DAKS | PINN |
| --- | --- | --- | --- |
| Nonlinear Elliptic | 8.40E-04 | 4.86E-05 | 4.20E-02 |
| Allen-Cahn ($a=15$) | 8.30E-02 | 6.06E-01 | 1.00E+00 |

### C.3 Irregularly-Shaped Domains

We tested on solving the nonlinear elliptic PDE (19) and the Allen-Cahn equation (21) with $a=15$. For the nonlinear elliptic PDE, the domain is an inscribed circle within $[0,1]\times[0,1]$. 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 in our paper, with boundary conditions derived from the solution.
For both PDEs, our method (SKS) used a virtual grid on $[0,1]\times[0,1]$ that covers the domain. For DAKS and PINN, we sampled the same number of collocation points from the domain. For a fair comparison, all the methods used the same set of 192 uniformly sampled collocation points on the boundary.
The error of each method is given in Table 6. The point wise error is shown in Fig. 4. As we can see, on irregularly-shaped domains, our method SKS still obtains a reasonably good accuracy for both cases (note that the Allen-Cahn case is much more challenging).

Figure: Figure 4: Point-wise solution error for Nonlinear Elliptic and 2D Allen-Cahn equation (21) with $a=15.0$ on irregularly-shaped domains.
Refer to caption: x4.png

**(a) Nonlinear Elliptic PDE (19).**
| Length-scale | 0.05 | 0.1 | 0.2 | 0.3 |
| --- | --- | --- | --- | --- |
| $L^{2}$ error | 1.19E-02 | 6.80E-06 | 4.62E-05 | 9.14E-04 |

### C.4 Sensitivity to Hyper-Parameters

To examine the sensitivity to the choice of kernel parameters, we run SKS to solve nonlinear Elliptic PDE (19) and Allen-Cahn equation (21) with $a=15$, with a varying set of length-scale parameters. For the nonlinear Elliptic PDE, we employed the grid of size $35\times 35$ while for the Allen-Cahan equation, we used the grid of size $49\times 49$. The results are 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^{2}$ error for solving the nonlinear elliptic PDE decreases from 1.129E-02 to 6.80E-06. Hence, our method is sensitive to the choice of the kernel parameters.

As a comparison, we also examined the sensitivity of PINN to the choice of the architectures. To this end, we fixed the depth at 8 and varied the layer width from 10 to 100, and also fixed the width at 30 while varying the depth from 3 to 30. We then tested PINN on solving the same PDEs. The solution error is reported in Table 8. 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.

**(a) $L^{2}$ error for depth=10 and layer width from 10 to 50.**
| Layer width | 10 | 20 | 30 | 40 | 50 |
| --- | --- | --- | --- | --- | --- |
| Nonlinear Elliptic | 1.75E-01 | 7.68E-02 | 7.01E-02 | 3.00E-02 | 4.06E-02 |
| Allan-cahn ($a=15$) | 6.31E+00 | 5.61E+00 | 7.08E+00 | 7.65E+00 | 1.03E+01 |

**(a) The Burgers’ equation (18) with $\nu=0.001$.**
| Method | 2400 | 4800 | 43200 |
| --- | --- | --- | --- |
| SKS (per-iter) | 4.6E-4 | 9.8E-4 | 6.8E-3 |
| DAKS (per-iter) | 7.43 | 38.5 | N/A |
| PINN (per-iter) | 2.7E-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 |

## Appendix 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 a 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 nor accelerate the convergence. We have also tried the relaxed Gauss-Newton approach as 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.

## Peer Reviews

confidence: 5
contribution: 4
initial_rating: 5
presentation: 4
questions: Apart from the points discussed above, there are the following minor points: 
- Around and after (2), there is some notational confusion between $z^j_m$ (a value) and $z^j$ (a function, later evaluated at $x_m$).
- Eq. (23) in Appendix A: The norm in the lhs should be over $\Omega$, not $\mathcal X$. Moreover $u_M$ should be $u_M^*$.
- Eq. (26), first row: Arguments $x_m$ are missing in the lhs.
- Starting from (7) and in (many of) the following occurrences, the index in the second sum starts from $m=M_+1$, instead of $m=M_{\Omega}+1$.
- Use of sampling inequalities: I don't understand why the one of (Batlle et al., 2023) is needed in eq. (21) and later, and not a sampling inequality on the flat Omega (as e.g. [6]).

[6] H. Wendland and C. Rieger, Approximate Interpolation with Applications to Selecting Smoothing Parameters, Numerische Mathematik (2005)
soundness: 2
strengths: - The topic is of current interest, and the results are clearly presented and discussed.
- The new method is computationally efficient and rather simple, when compared to other related approaches.
- A throughout theoretical analysis is provided, even if one step needs clarification.
summary: The paper introduces and studies a novel kernel-based method for the approximation of the solution of a broad class of non-linear PDEs. The authors discuss computational aspects and provide error estimates in Sobolev spaces of appropriate regularity.
The method is a modification (actually, a significant simplification) of an existing method, but the novelty is relevant as it provides a much more efficient solution method, and proves corresponding modified theoretical guarantees on the error.
weaknesses: - The paper lacks a sufficient discussion of the existing work, especially as it does not acknowledge a large body of literature on kernel-based (symmetric and non symmetric) collocation methods. In particular the last paragraph of Section 5 (Related Work) states that the novelty over existing methods is the use of nodal values instead of coefficients, and in the use of an optimization problem instead of the solution of a linear system. Both these aspects are treated in the literature, sometimes also for nonlinear problem, see e.g. [1-5], even if certainly not for such general problems as (1). I would suggest to expand the discussion of existing work, and also clearly articulate how the approach differs from or improves upon these existing methods.



- The relation between Assumption 4.1 and the existence and uniqueness of a strong solution of (1) is unclear. This relation would be needed to interpret Lemma 4.2 and Proposition 4.3 (i.e., given a PDE (1), which $k$, $t$, $\tau$ should one expect in (14)).
This fact is reflected also in the numerical experiments: It's unclear if these examples fit into the assumptions of the theoretical results (having strong solutions of sufficient regularity), and if so, with which values of $k$, $t$, $\tau$. In more general terms, the numerical experiments should address the expected convergence rates, other than just errors with fixed discretisations. 

- I'm not convinced by the argument leading to (32) in the proof of Lemma 4.2 in Appendix A. There is usually an issue in using a sampling inequality in this way. Namely, one has a domain $\mathcal M$ (in this case $\mathcal M=\mathcal T_i$), which has a corresponding critical value $h_0$ (depending on $\mathcal T_i$ via its diameter and boundary). Then, taking sufficiently many points in $\mathcal M$ one can have $h_i<h_0$ and thus apply a sampling inequality like the one of Proposition A.1 in (Batlle et al., 2023). Here however the only collocation point in $\mathcal T_i$ is $x_i$, and $h_i$ can not be made small without changing $\mathcal T_i$ itself, and thus $h_0$. Some oversampling inside $\mathcal T_i$ is usually needed. I suggest to clarify this point.

[1] K. Böhmer, R. Schaback, A Nonlinear Discretization Theory for Meshfree Collocation Methods applied to Quasilinear Elliptic Equations, ZAMM (2020)

[2] V. Bayona et al.,  RBF-FD formulas and convergence properties, Journal of Computational Physics (2010)

[3] I. Tominec, Residual Viscosity Stabilized RBF-FD Methods for Solving Nonlinear Conservation Laws, J Sci. Comp. (2022)

[4] Ka Chun Cheung et al, H^2-Convergence of Least-Squares Kernel Collocation Methods, SINUM (2018)

[5] N. Flyer et al., A guide to RBF-generated finite differences for nonlinear transport: Shallow water simulations on a sphere, J. Comp. Physics (2012)

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: If the authors can rewrite the manuscript more clearly, I would like change the grade.
soundness: 2
strengths: The topic looks very interesting
summary: The author proposes a very interesting kernel-based numerical solvers for solving nonlinear PDE.
weaknesses: 1. The manuscript has not been well-written, so that the reader cannot find their motivation clearly.
2. The theoretical part has been shown rigorously. 
3. The experiment part: description not clear

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 4
questions: - Why learn $\eta$, instead of $K_{MM}^{-1} \eta$? Both are vectors of the same dimension, so it should not make a diference. Evaluating the RHKS regularization term is also possible by calculating $(K_{MM}^{-1} \eta)^T K_{MM} K_{MM}^{-1} \eta$.
- In the numerical experiments, the collocation points should be the same for SKS and DAKS, as I am convinced the higher errors of DAKS arise from random sampling: with random sampling, there may be "holes" in the domain without collocation points, where the error is bigger.
soundness: 3
strengths: The paper is well presented. The literature review seems adequate and I especially appreciated the introduction to the previous approach by (Chen et al., 2021a). Furthermore, I found particularly clever the idea of placing collocation points on a grid and utilizing a product kernel, so that the kernel decomposes into a Kronecker product.
summary: This paper introduces a kernel learning framework for solving nonlinear partial differential equations (PDEs). Unlike a previous kernel solver (Chen et al., 2021a) that embeds differential operators within the kernel matrix, this new approach directly applies the operator to the kernel function, resulting in a hybrid method between traditional physics-informed neural networks (PINNs) and kernel methods. Additionally, the authors propose placing collocation points on a grid and utilizing a product kernel, so that the Gram matrix can be decomposed as a Kronecker product, significantly reducing computational costs.
weaknesses: I did not like the trade-off of making the kernel just a constant smaller (for instance, 1/3 smaller for the Burger's equation) at the cost of having to do an optimization similar to how regular PINNs are trained. Of course a 27 (3x3x3) times speed-up is not negligible, but with the authors approach, we lose the ability to do everything with fast Linear Algebra operations, and still do not get the advantages of training regular PINNs, where the colocation points can be randomly positioned and using stochastic gradients allows for using fewer collocation points per epoch. The product kernel trick was a nice idea, but I am convinced that it can also be applied in Chen et al., 2021a.

I am also skeptical of the timings reported by the authors in the appendix. In the experiments setup they say SKS is run using the ADAM optimizer for 1e6 epochs, but in the table 9 in the appendix, where they report time per iteration and total, the ratio is not 1e6. I also found it very surprising that the time per iteration is so small, especially compared to those of PINNs: I would think evaluating the derivative of the PINNs should be faster than evaluating this derivatives on the kernel version, even with the product trick.

---

confidence: 4
contribution: 1
initial_rating: 1
presentation: 2
questions: Based on the comments above, here are several questions that maybe useful for the authors if they decide to revise the paper:
1. Why solving the minimization problem in (6) is better than solving (2a) needs some clarifications? Can you provide numerical evidence in terms of accuracy and efficiency? The equation numbers here correspond to the equations labelled in the comments above.
2. Are there specific types of PDEs or solution structures for which the proposed tensorial kernel approach might be particularly well-suited?
3. Can you provide any theoretical insights or empirical evidence demonstrating advantages of this tensorial approach in terms of solution quality or generalization ability?
4. Can you clarify what DAKS is and what its relationship to their current work?
5. Can you provide additional state-of-the-art methods beyond PINN and finite-difference schemes?
soundness: 2
strengths: This work considers a slightly different way to solve PDE via the kernel-ridge-regression. See the assessment below.
summary: This work considers a slightly different way to solve PDE via the kernel-ridge-regression.  To put this work in context, see below.
weaknesses: To put this work in the context of RBF (e.g., [2,3]), let's consider solving a simple PDE problem,
$$
\\mathcal{P}(u) := -u'' + u = f, \\quad (1) %\\label{PDE}
$$
on a one-dimensional periodic domain (ignoring the boundary condition for simplicity). Using the notation in the paper, suppose we let the solution be
$$
u(x) = k(x,\mathcal{M}) \alpha, \quad (2)%\\label{ansatz}
$$
where $$\\mathcal{M} = \\{\mathbf{x}\_1,\ldots,\mathbf{x}\_M\\}$$ and $k$ denotes any kernel (which is positive definite and it corresponds to the RKHS space $\mathcal{H}$), which norm is defined as,
$$
\| u \|^2_{\mathcal{H}} = \alpha \mathbf{K} \alpha,
$$

where $\mathbf{K}= k(\mathcal{M},\mathcal{M})\in \mathbb{R}^{M\times M}$ is the Gram matrix of $k$. 

Inserting the ansatz in (2) into the PDE in (1), one can deduce the following linear system
$$
\mathbf{A}\alpha := (-\mathbf{K}''+c\mathbf{K}) \alpha = \mathbf{f}.
$$
where $\mathbf{K}''$ is a shorthand notation for the Gram matrix corresponds to $k"$ and $\mathbf{f} = (f(\mathbf{x}\_1),\ldots,f(\mathbf{x}\_M))$. At this point, let's just ignore the computational feasibility. Since $\mathbf{A}$ is invertible, one can simply solve this $M\times M$ linear problem. One can write the solution as $u(x) = k(x,\mathcal{M})\mathbf{A}^{-1}\mathbf{f}.$  I am not sure why solving this problem severely worse than the proposed method if the problem is invertible as noted in the last paragraph in Section 5. Many convergence results have been reported in literature (see [2,3] and the references therein).

Now, let's look at the Kernel Ridge Regression for this PDE problem, that is, we solve:
$$
\\min_{u\\in \\mathcal{H}}  \\frac{1}{M}\\sum_{i=1}^M (\\mathcal{P}(u)(\\mathbf{x}\_i) - f\_i)^2 +  \\lambda \\|u\\|^2_{\\mathcal{H}},
$$
Inserting the ansatz in in (2), we rewrite this optimization problem as,
$$
\min_\alpha \frac{1}{M} (\mathbf{A}\alpha - \mathbf{f})^\top (\mathbf{A}\alpha - \mathbf{f}) + \lambda \alpha \mathbf{K}\alpha. (2a)
$$
Taking derivative and set the equation to zero, we arrive at solving the following linear problem,
$$
(\mathbf{A}^2 + M\lambda \mathbf{K}) \alpha = \mathbf{A} f, \quad (3)%\label{linearproblem1}
$$
using the fact that $\mathbf{A}$ is symmetric. If this is invertible, one can simply write the solution as,
$$
u(x) = k(x,\mathcal{M})(\mathbf{A}^2 + M\lambda \mathbf{K})^{-1}\mathbf{A} f.\quad (4) %\label{sol1}
$$
The key idea in this paper is to consider the solution of the following form, 
$$
u(x) = k(x,\mathcal{M}) \mathbf{K}^{-1}\eta.\quad (5)%\label{ansatz2}
$$
In such case, the minimization problem becomes,
$$
\min_\eta \frac{1}{M} (\mathbf{A}\mathbf{K}^{-1}\eta - \mathbf{f})^\top (\mathbf{A}\mathbf{K}^{-1}\eta - \mathbf{f}) + \lambda \eta \mathbf{K}^{-1}\eta,\quad (6)
$$
and following the standard calculus, the solution satisfies,
$$
(\mathbf{K}^{-1}\mathbf{A}^2\mathbf{K}^{-1} + M \lambda \mathbf{K}^{-1}) \eta =\mathbf{K}^{-1} \mathbf{A}\mathbf{f},
$$
which is equivalent to multiplying (3) from the left by $\mathbf{K}^{-1}$ and letting $\alpha = \mathbf{K}^{-1}\eta$. If we solve this problem, we end up with
$$
u(x) = k(x,\mathcal{M})\mathbf{K}^{-1}(\mathbf{K}^{-1}\mathbf{A}^2\mathbf{K}^{-1} + M\lambda \mathbf{K}^{-1})^{-1}\mathbf{K}^{-1} \mathbf{A} f, %\quad (7)\label{soln2}
$$
which I believe is identical to (4) if $\mathbf{K}$ is invertible. In order to write the proposed ansatz in (5), $\mathbf{K}$ is assumed to be invertible. In fact, the author considers speeding up the inversion of $\mathbf{K}$ using an existing method in literature as reported in p.4.

While I fully understand that the authors consider a minimization algorithm to solve this problem, I just cannot see any new idea with this formulation. BTW, if the PDE is nonlinear, I believe RBF will result in a minimization problem as well. 

Beyond this issue, I also do not understand why the choice of tensorial kernel makes any sense except for computational convenience. Basically, the proposed idea in Eq.(9) in the manuscript is to consider the tensor product of kernels that compares scalars (a component of the data). So, the solutions are chosen in the space of functions induced by a kernel that only compares distances in one dimension. The classical approach is to choose a kernel that compares distances in $\mathbb{R}^d$. Mathematically, it is unclear why the proposed choice should always be a better choice aside from the numerical advantage.  

Finally, the numerical simulations are not convincing. I am not sure what DAKS is. It is not surprising the scheme is more accurate than PINN as it is well known that neural-network solutions are not accurate anyway. Lastly, the fact that this scheme beats the finite-difference scheme is also not so surprising based on the classical Numerical Analysis understanding in approximation of derivatives. One can check also the following paper that reported the advantage of RBF over finite-difference [1].

Based on these (comparisons with well-established classical literature), I am not sure this paper merits consideration for publication in ICLR.

References.
[1] B. Fornberg. The pseudospectral method: Comparisons with finite differences for the elastic wave equation.
 Geophysics, 52(4):483--501, 1987.

[2] B. Fornberg and N. Flyer.  A primer on radial basis functions with applications to the geosciences.
SIAM, 2015.

[3] B. Fornberg and N. Flyer. Solving pdes with radial basis functions. Acta Numerica, 24:215--258, 2015.

Meta-review: The paper introduces and studies a novel kernel-based method for the approximation of the solution of a broad class of non-linear PDEs. The authors discuss computational aspects and provide error estimates in Sobolev spaces of appropriate regularity. Unfortunately, two reviewers argue that the proposed methodology does not introduce fundamentally new ideas. The formulations presented appear to replicate well-established methods, particularly in the context of radial basis functions. One of these two reviewers is more positive, acknowledging that the method is a modification (actually, a significant simplification) of an existing method, but the novelty is relevant as it provides a much more efficient solution method, and proves corresponding modified theoretical guarantees on the error. Nevertheless, they point out that the paper lacks a sufficient discussion of the existing work, especially as it does not acknowledge a large body of literature on kernel-based (symmetric and non symmetric) collocation methods. Finally, all reviewers have serious concerns about the experiments.

Decision: Reject