# Nonparametric learning of kernels
in nonlocal operators

Fei Lu and Qingci An
Department of Mathematics
Johns Hopkins University
Baltimore, MD, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS>&Yue Yu
Department of Mathematics
Lehigh University
Bethlehem, PA, USA
<EMAIL_ADDRESS>
###### Abstract

Nonlocal operators with integral kernels have become a popular tool for
designing solution maps between function spaces, due to their efficiency in
representing long-range dependence and the attractive feature of being
resolution-invariant. In this work, we provide a rigorous identifiability
analysis and convergence study for the learning of kernels in nonlocal
operators. It is found that the kernel learning is an ill-posed or even ill-
defined inverse problem, leading to divergent estimators in the presence of
modeling errors or measurement noises. To resolve this issue, we propose a
nonparametric regression algorithm with a novel data adaptive RKHS Tikhonov
regularization method based on the function space of identifiability. The
method yields a noisy-robust convergent estimator of the kernel as the data
resolution refines, on both synthetic and real-world datasets. In particular,
the method successfully learns a homogenized model for the stress wave
propagation in a heterogeneous solid, revealing the unknown governing laws
from real-world data at microscale. Our regularization method outperforms
baseline methods in robustness, generalizability and accuracy.

##  1 Introduction

During the last 20 years there has been a lot of progress in the design of
machine learning models; however, their employment in scientific discovery of
hidden physical laws of complex system is relatively recent. One of the most
classical examples has been learning the material constitutive law directly
from experimental data, where the governing partial differential equation
(PDE) is unknown and the data features noises from measurements. This task can
be formulated as identifying an operator that continuously maps the
displacement field to the loading field. Moreover, such a need for learning
operators between function spaces has also become vital in other applications
ranging from homogenization problems [54, 53, 35, 34], fast PDE solvers [38,
39, 30, 32, 33, 25], to control problems [39, 21], just to name a few. Among
these work, models with nonlocal operators have have received increasing
attention, since they can describe physical phenomena that classical
differential operators fail to capture and provide a powerful model for a
large class of complex engineering and scientific applications [25, 46, 52,
6]. However, despite a long line of work on nonlocal model learning and its
applications, there is limited theoretical characterization of the underlying
inverse problem, even in the linear setting.

In this paper, we study the learning of kernels in linear nonlocal diffusion
operators from data. Suppose that we are given data

$\mathcal{D}=\\{u_{i},f_{i}\\}_{i=1}^{N}=\\{u_{i}(x_{j}),f_{i}(x_{j}):j=1,\ldots,J\\}_{i=1}^{N},$
(1)

where $u_{i},f_{i}$ are real-valued continuous functions on a bounded open
connected set $\Omega\subset\mathbb{R}^{d}$ and $\\{x_{j}\in\Omega\\}$ are
spatial mesh points. The task is to learn the kernel $\phi$ in a nonlocal
diffusion operator $L_{\phi}$ mapping $u_{i}$ to $f_{i}$ of the form:

$L_{\phi}[u](x)=\int_{\Omega}\phi(|y-x|)[u(y)-u(x)]dy=f(x),\forall
x\in\Omega.$ (2)

This operator is simple yet flexible for nonlocal models: it has been employed
in [54, 53] to provide a homogenized material model from microscale
measurements, and can be seen as a simplified linear version of the integral
neural operator in [51].

Our goal is to infer the kernel $\phi$ from data via nonparametric regression,
so as to address the general situations that there is limited information to
derive a parametric form or constraints for the kernel, which can be either
smooth or singular. The regression utilizes the linear dependence of the
operator on the kernel, making it possible to treat the large size functional
data in a scalable fashion.

Three challenges are to be overcome. First, the function space of
identifiability (FSOI) is yet to be specified properly, otherwise the inverse
problem can be ill-defined in the sense that there are multiple kernels
fitting the data. This is fundamentally different from classical nonparametric
regression that learns a function $Y=\phi(X)$ from random samples
$\\{(X_{i},Y_{i})\\}$ from the joint distribution of $(X,Y)$, for which the
FSOI is $L^{2}(\rho)$ with $\rho$ being the distribution of $X$ and the
optimal estimator is the conditional expectation. Second, the kernel estimator
should be resolution independent and converge in a proper function space when
the data resolution refines, so that it can be applied to problems and
simulation tasks with different grids or discretization methods and provides a
guaranteed modeling accuracy. Third, beyond resolution invariance, the
estimator should be robust to imperfect data so as to be applicable in real
applications.

We introduce a nonparametric learning method producing a convergent estimator
for the kernel as the data mesh refines. Our main result has two folds.
Firstly, we specify a data adaptive function space of learning and provide an
identifiability theory, which shows that the function space of identifiability
is the closure of a system intrinsic data adaptive reproducing kernel Hilbert
space (SIDA-RKHS), which is equipped with an exploration measure indicating
the information from data. Beyond this space, the inverse problem may become
undetermined. Furthermore, the identifiability theory shows that the inverse
problem is ill-posed and it becomes more ill-posed when the data resolution
refines, which contradicts the intuition that learning from more datasets and
refined mesh grids should provide an improved estimator. Thus, regularization
is necessary. Second, we introduce a new regularization method that plays a
key role in identifying a convergent estimator in the presence of model errors
and/or measurement noises. It uses the norm of the SIDA-RKHS. In experimental
studies, we compare our proposed SIDA-RKHS regularization method with two
common Tikhonov/ridge regularizers that use $l^{2}$ and $L^{2}$ norms. Results
on both benchmark problems with synthetic data and real-world data show that
only the SIDA-RKHS regularizer can consistently obtain convergent estimators
for all types of kernels, especially when the data is noisy.

We summarize our major contributions below:
1) We establish a rigorous identifiability theory for the nonparametric
learning of kernels in nonlocal operators, and for the first time specifying a
system-intrinsic data-adaptive function space of identifiability (see Lemma
2.3 and Theorem 2.1). The theory also indicates a pitfall of the nonlocal
kernel learning problem: the inverse problem is ill-posed.
2) We introduce a nonparametric regression algorithm equipped with a novel
regularization method based on the SIDA-RKHS (see Section 2.3), which
overcomes the ill-posedness to yield a convergent estimator robust to noise.
3) We validate the theory and the proposed algorithm on a number of benchmark
problems, including various synthetic datasets and a real-world dataset where
the governing law is unknown (see Section 3). Results show that the proposed
algorithm provides a stable and converging estimator, while the common
Tikhonov/ridge regularizers with $l^{2}$ or $L^{2}$-norm fail this task.

###  1.1 Related Work

Nonlocal operators: Nonlocal operators arise in various areas such as nonlocal
and fractional diffusions [45, 12, 11, 2, 1, 5, 8, 50, 48], de-noising and
regularization by nonlocal kernels [24, 15, 19], multi-agent systems with
nonlocal interaction [37, 36, 26] and nonlocal networks [49, 31]. The inverse
problem for nonlocal diffusions has been studied in [22, 28] from a single
solution. To discover hidden nonlocal physical laws from data, a parametric
nonlocal kernel learning approach has been proposed in [52, 53], where the
coefficients of Bernstein polynomials are learnt with physics-based
constraints and a Tikhonov regularization. Beyond the linear nonlocal model
and the nonlocal kernel regression methods, nonlocal operators were further
combined with neural networks, and nonlocal kernel networks were developed for
learning maps between high-dimensional variables in dynamical systems [33, 32]
or function spaces [30, 51]. An attractive feature of these nonlocal
kernel/operator learning methods is the generalizability among approximations
corresponding to different underlying levels of resolution and discretization.
However, as seen in [54, 53, 30, 51], neither the nonlocal kernel learning
methods nor nonlocal kernel networks yield estimator convergence when trained
on finer resolution, and the test error may even increase. This fact indicates
the possible ill-posedness of the learning problem. In this work, we tackle
this issue by introducing a new regularization method based on a data-adaptive
RKHS in a nonparametric learning approach.

Functional data analysis: Functional data analysis (see e.g., [20, 23, 14] and
the references therein) studies the learning an infinite-dimensional operator
from functional data. In contrast, we focus on learning a radial kernel in an
operator, exploiting the low-dimensional structure of the operator, which
enables us to learn the kernel (hence the operator) from limited data.

Regularization methods: Our SIDA-RKHS regularization is a type of
Tikhonov/ridge regularization that adds a penalty term to the loss function.
It differs from previous methods at the penalty term. The commonly used
penalty terms include the Euclidean norm in the classical Tikhonov
regularization [18, 17], the RKHS norm with an ad hoc reproducing kernel
(often the Gaussian kernel) [9, 10], the total variation norm in the Rudin-
Osher-Fatemi method [43], or the $L^{1}$ norm in LASSO [47]. Whereas each of
these penalty terms has their specific applications, none of them take into
account of the FSOI, which is fundamental for learning kernels in operators.
Also, our regularization method is inspired by the kernel flow method that
learns hyper-parameters of the reproducing kernel [41, 16, 7], but our
reproducing kernel is determined by the system and the data. Given the
importance of regularization to overcome ill-posedness and overfitting, we
expect our SIDA-RKHS regularization method to be applicable to a wide range of
linear inverse problems and machine learning methods.

##  2 Learning theory and algorithm

###  2.1 Nonparametric regression with regularization

We construct an estimator by minimizing the loss functional of mean square
error:

$\mathcal{E}(\phi)=\frac{1}{N}\sum_{i=1}^{N}\|L_{\phi}[u_{i}]-f_{i}\|_{L^{2}}^{2}.$
(3)

Here we consider only the $L^{2}$ norm is
$\|f\|_{L^{2}}^{2}=\int|f(x)|^{2}dx$, which has minimal requirements on the
data. Other norms (e.g., the Sobolev norms) can also be used when the data are
smooth.

Note that the loss functional is quadratic in $\phi$ because the nonlocal
operator is linear in $\phi$. Thus, the minimizer of the loss functional is
the least squares estimator (LSE), which is handy once one selects a set of
basis functions for a hypothesis space. More specifically, suppose the
hypothesis space is $\mathcal{H}_{n}=\mathrm{span}\\{\phi_{i}\\}_{i=1}^{n}$,
for each $\phi=\sum_{i=1}^{n}c_{i}\phi_{i}\in\mathcal{H}_{n}$, we can write
the loss functional in (3) as
$\mathcal{E}(c)=\mathcal{E}(\phi)=c^{\top}{\overline{A}}_{n}c-2c^{\top}{\overline{b}}_{n}+C_{N}^{f}$,
where $C_{N}^{f}=\frac{1}{N}\sum_{k=1}^{N}\int|f_{i}(x)|^{2}dx$, and the
normal matrix ${\overline{A}}$ and vector ${\overline{b}}$ are given by

${\overline{A}}_{n}(i,j)=\langle\langle{\phi_{i},\phi_{j}}\rangle\rangle,\,{\overline{b}}_{n}(i)=\frac{1}{N}\sum_{k=1}^{N}\int
L_{\phi_{i}}[u_{k}](x)f_{k}(x)dx,$ (4)

and the bilinear form $\langle\langle{\cdot,\cdot}\rangle\rangle$ is defined
by

${\langle\langle{\phi,\psi}\rangle\rangle=\frac{1}{N}\sum_{k=1}^{N}\int_{\mathbb{R}^{d}}L_{\phi}[u_{k}](x)L_{\psi}[u_{k}](x)dx.}$
(5)

The least squares estimator is computed directly from the minimizer of the
quadratic function $\mathcal{E}(c)$:

$\widehat{\phi}_{\mathcal{H}_{n}}=\sum_{i}\widehat{c}_{i}\phi_{i},\quad\text{
where }\widehat{c}={\overline{A}}_{n}^{-1}{\overline{b}}_{n},$ (6)

where ${\overline{A}}_{n}^{-1}$ is the inverse (or pseudo-inverse when the
inverse does not exist) of ${\overline{A}}_{n}$.

However, the above least squares regression encounters a big challenge in
obtaining convergent estimators for this ill-posed inverse problem (see
Section 2.2). As a nonparametric method, it is often necessary to select a
relatively large hypothesis space to make the model flexible enough. However,
the large hypothesis space leads to a normal matrix that is often severely
ill-conditioned. As a result, the estimator in (6) oscillates violently when
the data is imperfect due to either measurement noise or model error, and the
estimator does not converge when the data mesh refines.

Regularization methods overcome the ill-posedness by adding a penalty term to
the loss functional:

$\displaystyle\mathcal{E}_{\lambda}(\phi)=\mathcal{E}(\phi)+\lambda\mathcal{R}(\phi),$
(7)

where $\mathcal{R}(\phi)$ is a regularization term, and $\lambda$ is a hyper-
parameter controlling the contribution of the regularization term. Various
penalty terms have been proposed, however, none of them take into account of
the function space of identifiability, which is at the foundation of learning
(see Section 2.2). Based on it, we will introduce a data-adaptive RKHS
regularization method (in Section 2.3). Thus, it is different from classical
regularization using an ad hoc RKHS [10, 3].

###  2.2 Function space of identifiability

The identifiability theory characterizes the function space of learning. There
are two key elements in our identifiability theory: 1) an exploration measure,
which is a probability measure that quantifies the exploration of the kernel’s
variable by the data, and 2) the function space of identifiability, in which
the loss functional has a unique minimizer. They are described as follows.

The exploration measure. As the first key element, we introduce first a novel
measure on $\mathbb{R}_{+}$ that quantifies the exploration of the independent
variable of the kernel by the data. We assume the radial kernel’s support to
be in an interval $[0,R_{0}]$. A given dataset may only explore part of this
interval. More specifically, the discrete data set in (1) explores only the
pairwise distances $|x_{j}-x_{k}|$ in
$\mathcal{R}_{N}^{J}=\\{r_{ijk}=|x_{j}-x_{k}|\leq
R_{0}:u_{i}(x_{j})-u_{i}(x_{k})\neq 0\text{ for some }i,j,k\\}$, the set of
all the pairwise distances $|x_{j}-x_{k}|$ with repetition. We define an
empirical measure and its continuous limit

$\displaystyle\rho_{N}^{J}(dr)$
$\displaystyle=\frac{1}{|\mathcal{R}_{N}^{J}|}\sum_{i=1}^{N}\sum_{j,k=1}^{J}\delta_{|x_{j}-x_{k}|}(r)w_{i}(x_{j},x_{k}),$
(8) $\displaystyle\rho_{N}(dr)$
$\displaystyle=\frac{1}{ZN}\sum_{i=1}^{N}\int_{\Omega}\int_{\Omega}\delta_{|x-y|}(r)w_{i}(x,y)dxdy$

for $r\in[0,R_{0}]$, where $|\mathcal{R}_{N}^{J}|$ is the cardinality of the
set $\mathcal{R}_{N}^{J}$, $\delta_{s}(r)$ is the Dirac distribution with
point mass at $s$, and $Z$ is the normalizing constant. Here the weight
function is $w_{i}(x,y)=|u_{i}(x)-u_{i}(y)|$.

The exploration measure plays an important role in the learning of the kernel.
It reflects the strength of exploration to $|x-y|$ by the data
$|u_{i}(x)-u_{i}(y)|$ in the loss function, and it will act as a re-weighting
factor through the SIDA-RKHS regularization to be introduced in Section 2.3.
Thus, we will use it to quantify the accuracy of the kernel’s estimator in
$L^{2}(\rho_{N})$ (or $L^{2}(\rho_{N}^{J})$ for discrete data).

Main result: function space of identifiability. We define the function space
of identifiability (FSOI) as the largest linear space in which the loss
functional has a unique minimizer. In other words, the variational inverse
problem of finding a unique minimizer of the loss functional is well-defined
in this space. In the following, we write only the continuous function space
$L^{2}(\rho_{N})$, but all the arguments apply to the discrete function space
$L^{2}(\rho_{N}^{J})$ in an obvious manner (see Remark 2.4).

######  Theorem 2.1 (Function space of identifiability).

Consider the problem of learning the kernel $\phi$ by minimizing the loss
functional $\mathcal{E}$ in (3) with $\\{u_{i},f_{i}\\}_{i=1}^{N}$ being
continuous in a bounded domain $\Omega$. Then, the function space of
identifiability (FSOI), the largest subspace of $L^{2}(\rho_{N})$ in which
$\mathcal{E}$ has a unique minimizer, is the eigen-space of nonzero
eigenvalues of ${\mathcal{L}_{\overline{G}}}$, an integral operator defined by

${\mathcal{L}_{\overline{G}}}\phi(r)=\int_{0}^{\infty}\phi(s)\overline{G}(r,s)\rho_{N}(ds).$
(9)

Here the integral kernel ${\overline{G}}$ comes from data:

${\overline{G}}(r,s)=[\rho_{N}^{\prime}(r)\rho_{N}^{\prime}(s)]^{-1}G(r,s),$
(10)

where $\rho_{N}^{\prime}$ is the density of $\rho_{N}$ and $G$ is

$\displaystyle G(r,s)=\frac{1}{N}\sum_{i=1}^{N}$
$\displaystyle\int_{|\eta|=1}\int_{|\xi|=1}\left[\int[u_{i}(x+r\xi)-u_{i}(x)]\right.\left.[u_{i}(x+s\eta)-u_{i}(x)]dx\right]d\xi
d\eta,$ (11)

for $r,s\in\mathrm{supp}(\rho_{N})$, and $G(r,s)=0$ otherwise. Furthermore,
the minimizer of $\mathcal{E}$ is

$\widehat{\phi}={\mathcal{L}_{\overline{G}}}^{-1}P\phi_{N}^{f},$

where $P$ is the projection to the FSOI. Here $\phi_{N}^{f}\in
L^{2}(\rho_{N})$ is the Riesz representation of the bounded linear functional
defined by
$\langle{\phi_{N}^{f},\psi}\rangle_{L^{2}(\rho_{N})}=\frac{1}{N}\sum_{i=1}^{N}\int
2L_{\psi}[u_{i}](x)f_{i}(x)dx,\,\forall\psi\in L^{2}(\rho_{N}).$

When the data is continuous and noiseless, the true kernel is the unique
minimizer, i.e., it is identifiable by the loss functional, if it is in the
FSOI. Also, note that the FSOI is data-dependent. When the data is discrete or
noisy, the unique minimizer is an optimal estimator in the FSOI, and it
converges to the true kernel as the data mesh refines (see Remark 2.4) and the
noise to signal ratio reduces.

The proof of Theorem 2.1 is based on the uniqueness of zero of the Fréchet
derivative of the loss functional, which becomes clear from the following
lemma. Their proofs are deferred to Appendix A.

######  Lemma 2.2 (The Fréchet derivative of the loss functional).

The Fréchet derivative of the loss functional $\mathcal{E}$ in
$L^{2}(\rho_{N})$, with ${\mathcal{L}_{\overline{G}}}$ defined in (9) and
$\phi_{N}^{f}$ defined in Theorem 2.1, is

$\nabla\mathcal{E}(\phi)=2({\mathcal{L}_{\overline{G}}}\phi-\phi_{N}^{f}).$

System-intrinsic data-adaptive RKHS. Theorem 2.1 highlights two fundamental
challenges: the inverse problem is well-defined only in the FSOI, and it is
ill-posed in the FSOI because it involves the inverse of a compact operator
${\mathcal{L}_{\overline{G}}}$ (as shown in the next lemma). Fortunately, the
integral kernel ${\overline{G}}$ defines a reproducing kernel Hilbert space
(RKHS), which provides a regularization norm to ensure the learning to take
place in the FSOI and to overcome the ill-posedness. This RKHS is system
intrinsic as it depends on the structure of the system of nonlocal operators,
and it is data-adaptive, utilizing both the exploration measure and the data
$\\{u_{i}\\}$. Thus, we call it SIDA-RKHS.

######  Lemma 2.3 (Characterization of the SIDA-RKHS).

Suppose that the data $\\{u_{i}\\}$ are continuous in $\Omega$. Then, the
following statements hold true.

  * (a)

The integral kernel ${\overline{G}}$ defined in (10) is positive semi-
definite.

  * (b)

The integral operator ${\mathcal{L}_{\overline{G}}}:L^{2}(\rho_{N})\to
L^{2}(\rho_{N})$ defined in (9) is compact and positive semi-definite, and we
have, for any $\phi,\psi\in L^{2}(\rho_{N})$,

$\langle\langle{\phi,\psi}\rangle\rangle=\langle{{\mathcal{L}_{\overline{G}}}\phi,\psi}\rangle_{L^{2}(\rho_{N})}.$
(12)

  * (c)

The RKHS $H_{G}$ with ${\overline{G}}$ as reproducing kernel satisfies
$H_{G}={\mathcal{L}_{\overline{G}}}^{1/2}(L^{2}(\rho_{N}))$, and its inner
product satisfies
$\langle{\phi,\psi}\rangle_{H_{G}}=\langle{{\mathcal{L}_{\overline{G}}}^{-1/2}\phi,{\mathcal{L}_{\overline{G}}}^{-1/2}\psi}\rangle_{L^{2}(\rho_{N})}$
for any $\phi,\psi\in H_{G}$.

  * (d)

The eigenvalues of ${\mathcal{L}_{\overline{G}}}$ converges to zero, and its
eigen-functions $\\{\psi_{k}\\}_{k}$ form a complete orthonormal basis of
$L^{2}(\rho_{N})$. For any $\phi=\sum_{k}c_{k}\psi_{k}$, we have

$\displaystyle\langle\langle{\phi,\phi}\rangle\rangle=\sum_{k}\lambda_{k}c_{k}^{2},\quad\|\phi\|^{2}_{L^{2}(\rho_{N})}=\sum_{k}c_{k}^{2},\quad\|\phi\|^{2}_{H_{G}}=\sum_{k}\lambda_{k}^{-1}c_{k}^{2},$
(13)

where the last equation is restricted to $\phi\in H_{G}$.

Algorithm 1 Nonparametric learning of the nonlocal kernel with SIDA-RKHS
regularization

Input: The data
$\\{u_{i},f_{i}\\}_{i=1}^{N}=\\{u_{i}(x_{j}),f_{i}(x_{j})\\}_{i,j=1}^{N,J}$ to
construct the nonlocal model $L_{\phi}[u]=f$.

Output: Estimator $\widehat{\phi}$

1\. Estimate the exploration measure $\rho_{N}^{J}$ as in (8), and denote $R$
the upper bound of its support.

2\. Get regression data (see Appendix B).

3\. Select a class of hypothesis spaces
$\mathcal{H}_{n}=\mathrm{span}\\{\phi_{i}\\}_{i=1}^{n}$ with $n$ in a proper
range.

4\. For $n$ in the range

4a) Compute $({\overline{A}}_{n},{\overline{b}}_{n},B_{n})$ for
$\mathcal{H}_{n}=\mathrm{span}\\{\phi_{i}\\}_{i=1}^{n}$ with
$B_{n}=(\langle{\phi_{i},\phi_{j}}\rangle_{L^{2}(\rho_{N}^{J})})_{1\leq
i,j\leq n}$;

4b) If the basis matrix $B_{n}$ is singular, stop and remove $n$ from the
range;

4c) Solve the generalized eigenvalue problem ${\overline{A}}_{n}V=B_{n}\Lambda
V$, where $\Lambda$ is the diagonal matrix of eigenvalues and
$V^{\top}B_{n}V=I_{n}$;

4d) Compute the RKHS-norm matrix $B_{rkhs}=(V\Lambda V^{\top})^{-1}$;

4e) Use the L-curve method find an optimal estimator
$\widehat{\phi}_{\lambda_{n}^{*}}$.

5\. Select the optimal dimension $n^{*}$ (and degree if using B-spline basis)
that has the minimal loss value (along with other cross-validation criteria if
available). Return the estimator
$\widehat{\phi}=\sum_{i=1}^{n^{*}}c^{i}_{n^{*}}\phi_{i}$.

######  Remark 2.4 (Discrete data).

When the space $L^{2}(\rho_{N}^{J})$ is a discrete vector space due to
discrete data, we learn the kernel on finitely many points
$\\{r_{k}\\}_{k=1}^{n}$ explored by the data. In this case, the integral
kernel $G$ in (11) becomes a positive semi-definite matrix in
$\mathbb{R}^{n\times n}$, so is ${\overline{G}}$ in (10). Now the operator
${\mathcal{L}_{\overline{G}}}$ is defined by the matrix ${\overline{G}}$ on
the weighted vector space $\mathbb{R}^{n}$ and its eigenvalues is the
generalized eigenvalue of $({\overline{G}},B_{n})$ with $B_{n}$ being the
diagonal matrix of $\rho_{N}^{J}$. As a result, the SIDA-RKHS $H_{G}$ is the
vector space spanned by the eigenvectors with nonzero eigenvalues.
Furthermore, its norm in (13) can be computed directly from the eigen-
decomposition. This norm is better suited for regularization even when the
SIDA-RKHS has the same dimension as $L^{2}(\rho_{N}^{J})$ (or dense in it). As
data mesh refines, these vector spaces converges to the corresponding function
spaces when the data is smooth.

###  2.3 Algorithm: LSE with SIDA-RKHS regularization

Based on the function space of identifiability, we introduce next a
nonparametric learning algorithm with SIDA-RKHS regularization. The algorithm
consists of three steps. First, we utilize the data to estimate the
exploration measure and the support of the kernel. Based on them, we set a
class of hypothesis spaces, with their dimensions, i.e., the number of basis
functions, in a proper range moving from under-fitting to over-fitting. For
the hypothesis space $\mathcal{H}_{n}=\mathrm{span}\\{\phi_{i}\\}_{i=1}^{n}$,
we compute the basis matrix
$B_{n}=(\langle{\phi_{i},\phi_{j}}\rangle_{L^{2}(\rho_{N}^{J})})_{1\leq
i,j\leq n}\in\mathbb{R}^{n\times n}$. Second, we assemble the regression
matrices from data for each of these hypothesis spaces. We approximate the
integrals by Riemann sum or other numerical integrator. Finally, we identify
an estimator with SIDA-RKHS regularization for each of these hypothesis spaces
by the L-curve method [17] and select the one with the best fitting. We
summarize the method in Algorithm 1, with its full details provided in Section
B.

The core innovations are the exploration measure and the regularization using
the SIDA-RKHS norm. Importantly, they bring little extra computational cost.
The exploration measure is available directly from data. The SIDA-RKHS norm is
computed directly from the triplet
$({\overline{A}}_{n},{\overline{b}}_{n},B_{n})$ using the generalized
eigenvalue problem as detailed in the algorithm.

Our SIDA-RKHS regularization uses the RKHS norm
$\mathcal{R}(\phi)=c^{\top}B_{rkhs}c$, where $B_{rkhs}$ is defined in (4d) in
Algorithm 1. It differs from the commonly-used Tikhonov/ridge regularization
using either the $l^{2}$-norm that sets $\mathcal{R}(\phi)=\sum_{i}c_{i}^{2}$
or the $L^{2}(\rho_{N}^{J})$-norm that sets
$\mathcal{R}(\phi)=c^{\top}B_{n}c$. We note that the three norms become the
same when $B_{n}=I_{n}$ and all the eigenvalue of $A_{n}$ are $1$.

##  3 Tests on synthetic and real-world data

We test our nonparametric learning method on both synthetic data and real-
world data in 1D examples. On each dataset, we compare our SIDA-RKHS
regularizer with two baseline regularizers using the $l^{2}$ and $L^{2}$ norm
(denoted as l2 and L2, respectively). All these regularizers use the same
L-curve method to select the hyper-parameter $\lambda$ as described in
Appendix C.2. In the case of synthetic data, we systematically examine the
method with three types of kernels in the regimes of noiseless and noisy data.
Since the ground-truth kernel is known, we study the convergence of estimators
to the true kernel as the data mesh refines. We also apply our method to a
real-world dataset for stress wave propagation in a heterogeneous bar, with
the goal of constructing a homogenized model from microscale data. Since there
is no ground-truth, we examine the performance of estimators by studying their
physical stability and capability of reproducing the wave motion on a cross-
validation dataset. All datasets and codes used will be publicly released on
GitHub.

Settings for the learning algorithm. In implementation of Algorithm 1, we use
B-spline basis functions consisting of piece-wise polynomials with degree 2 so
that the estimated kernel is twice differentiable (see Section C.1 for a brief
introduction of B-splines). The knots of B-splines are evenly spaced on
interval $[0,R]$, with one additional knots at $0$ to make the first basis
nonzero at $x=0$. We select the dimension with minimal loss from a sequence of
dimensions in the range $\lfloor{\frac{R}{\Delta x}}\rfloor\times[0.2,1]$ as
long as the basis matrix $B_{n}$ is well-conditioned.

###  3.1 Examples with synthetic data

Numerical settings. We consider three kernels: a sine kernel, a Gaussian
kernel, and a fractional Laplacian kernel (specified below). They act on the
same set of functions $\\{u_{i}\\}_{i=1,2}$ with
$u_{1}=\sin(x)\mathbf{1}_{[-\pi,\pi]}(x)$ and
$u_{2}(x)=\cos(x)\mathbf{1}_{[-\pi,\pi]}(x)$. In the ground-truth model, the
integral $L_{\phi}[u_{i}]$ is computed by the adaptive Gauss-Kronrod
quadrature method, which is much more accurate than the Riemann sum integrator
that we will use in the learning stage. To create discrete datasets with
different resolutions, for each $\Delta x\in 0.0125\times\\{1,2,4,8,16\\}$, we
take values
$\\{u_{i},f_{i}\\}_{i=1}^{N}=\\{u_{i}(x_{j}),f_{i}(x_{j}):x_{j}\in{[-40,40]},j=1,\ldots,J\\}_{i=1}^{N}$,
where $x_{j}$ is a point on the uniform grid with mesh size $\Delta x$.

For each kernel, we consider both noiseless and noisy data with different
noise levels, with a noise-to-signal-ratio ($nsr$) taking values
$\\{0,0.5,1,2\\}$. Here the noise is added to each spatial mesh point,
independent and identically distributed centered Gaussian with standard
deviation $\sigma$, and the noise-to-signal-ratio is the ratio between
$\sigma$ and the average $L^{2}$ norm of $f_{i}$.

The three ground-truth kernels are specified as follows.
$\bullet$ _Sine kernel._ The sine kernel is
$\phi_{true}(r)=\sin(6r)\mathbf{1}_{[0,10]}(r)$. This sine kernel represent a
smooth oscillating kernel in the same class as the data $u_{i}$. The estimated
support is in $[0,R]$ with $R=11.02$.
$\bullet$ _Gaussian kernel._ The Gaussian kernel $\phi_{true}$ is the Gaussian
density centered at 5 with standard deviation 1. This kernel represents a
smooth kernel. It has $R=11.58$.
$\bullet$ _Fractional Laplacian kernel._ It is a truncated version of the
fractional Laplacian kernel that has been widely studied in fractional and
nonlocal diffusions (see e.g., [5, 2, 12, 48]). We set
$\phi_{true}(r)=c_{d,s}r^{-(d+2s)}\mathbf{1}_{[0.1,6]}(x)+10^{d+2s}\mathbf{1}_{[0,0.1]}(x)$
with exponent $s=0.5$ and $d=1$, where
$c_{d,s}=4^{s}\pi^{-d/2}\Gamma(d/2+s)\Gamma(-s)$. It is almost singular with
multiscale values and its values near the singularity are crucial to the
operator. It has $R=6.51$.

Figure 1: (a): Typical estimators from noisy data with noise-to-signal-ratio
$\mathrm{nsr}=1$ and $\Delta x=0.025$. (b): the relative $L^{2}(\rho_{N}^{J})$
errors of these estimators. Bold numbers highlight the best method. The SIDA-
RKHS regularizer consistently obtains accurate estimators in all three cases.

Performance of the regularizers. We present the typical estimators and the
convergent rate of the estimator as data mesh refines. Figure 1 shows typical
estimators for the three examples from noisy data with a noise-to-signal ratio
nsr=1 and $\Delta x=0.025$. The hypothesis space’s dimension is selected by
minimal loss value. All three regularizers are able to estimate the Sine
kernel accurately and the Fractional Laplacian kernel reasonably. The SIDA-
RKHS regularizer significantly outperforms the regularizers with $l^{2}$ or
$L^{2}$-norm in the example of the Gaussian kernel.

Figure 2: (a) The means and standard deviations of the convergence rates as
mesh refines in 100 independent simulations. The SIDA-RKHS regularizer obtains
consistent rates for noisy data. (b) Error decay as the number of data pairs
increases when $nsr=1$ and $\Delta x=0.0125$.

The SIDA-RKHS regularizer’s superior performance is further validated by the
rates of convergence when $\Delta x$ decreases, from 100 independent
simulations with noises with $nsr\in\\{0,0.5,1,2\\}$, as shown in Figure 2.
The SIDA-RKHS regularizer has rates generally higher than those of the other
two regularizers when the data gets more noisy. We note that the SIDA-RKHS
regularizer has deceivingly lower rates for noiseless data, even though it
actually has more accurate estimators (see Figure 4 in appendix). Thus, the
goal is to seek an accurate estimator with a consistent rate. Here the rates
for the smooth kernels are higher than the rate for the singular kernel,
because the order of numerical error in the Riemann sum integrator are higher
(see [13]).

Increasing the number of data pairs. Sine the operator is linear, only
linearly independent data brings new information for the learning. Figure 2(b)
shows that as $N$ of data
$\\{u_{i}(x)\\}_{i=1}^{N}=\\{\sin(ix),\cos(ix)\\}_{i=1}^{N/2}$ increases, the
estimators become more accurate but without a convergence rate. Note that the
data pairs do not provide independent (random) samples of $\rho_{N}$, which
also varies with data. Thus, our learning problem is fundamentally different
from regression for random samples and we do not expect a convergence rate
$N^{-1/2}$. An interesting future direction is to design experiments to
collect informative data to enlarge the FSOI and accelerate the convergence.

In summary, the SIDA-RKHS regularizer consistently obtains accurate convergent
estimators when data mesh refines for either noiseless or noisy data. On the
contrary, the regularizers with $l^{2}$ norm or $L^{2}$ norm, are not robust
to noise and may fail to converge, due to their negligence of the FSOI.

###  3.2 Homogenization of wave propagation in meta-material

Figure 3: Real-world application: wave propagation in a heterogeneous bar with
ordered microstructure of period $L=0.4$, and the estimated support of the
kernel has a bound $R=1.65$.

We seek a nonlocal homogenized model for the stress wave propagation in a one-
dimensional heterogeneous bar with a periodic microstructure. For this
problem, the goal is to obtain an effective surrogate model from high-fidelity
(HF) datasets generated by solving classical wave equation, acting at a much
larger scale than the size of the microstructure. Differing from previous
examples, this problem has no ground-truth kernel. Therefore, we evaluate the
estimator by measuring its effectiveness of reproducing HF data in
applications that are subject to different loading conditions with a much
longer time from the problems used as training data.

For both training and validation purposes we use the HF dataset generated by
the direct numerical solver (DNS) introduced in [44], which provides exact
solutions of velocities including the appropriate jump conditions for the
discontinuities in stress that occur at waves. Although the DNS has high
accuracy on wave velocity, it is not suitable for long-term prediction because
it requires the modeling of wave propagation through thousands of
microstructural interfaces, which makes the computational cost prohibitive. To
accelerate the computation, we approximate the HF model by a nonlocal model:

$\partial_{tt}{u}(x,t)-L_{\phi}[{u}](x,t)=g(x,t),\text{ for
}(x,t)\in\Omega\times[0,T],$ (14)

where $L_{\phi}$ is a nonlocal operator in the form of (2) with a kernel
$\phi$ being supported in $[0,R]$.

Experiment settings. We consider four types of data: three for training and
one for validation of our algorithm. Three types of training datasets are
employed: In Type 1 dataset, the bar is subject to an oscillating source
$g(x,t)$; In Type 2 dataset, a boundary velocity loading
$\partial_{t}{u}(-50,t)=\cos(jt)$ is applied; In Type 3 dataset, all settings
are the same as in Type 2, except that the $\cos(jt)$ type loading is replaced
by $\sin(jt)$. In all training datasets we consider a relatively small domain
$\Omega=[-50,50]$ and short time $t\in[0,2]$. Two spatial resolutions, $\Delta
x=0.05$ and $\Delta x=0.025$ are considered, which we denote as the “coarse”
and “fine” datasets, respectively.

With these three types of training datasets, we design three experiment
settings to validate our method:
$\bullet$ _Coarse dataset_ 1: we train the estimator using “coarse” dataset of
Types 1 and 2.
$\bullet$ _Coarse dataset_ 2: we train the estimator using “coarse” dataset of
Types 1 and 3. By comparing the learnt estimator from this setting with the
result from setting 1, we mean to investigate the sensitivity of the inverse
problem with respect to the choice of datasets.
$\bullet$ _Fine dataset:_ we train the estimator using “fine” dataset of Types
1 and 2. By comparing the learnt estimator from this setting with the result
from setting 1, we aim to check the convergence of the estimator with
increasing data resolution. Note that the problem might becomes more ill-posed
when decreasing $\Delta x$. Therefore, proper regularization is expected to
become more important.

Additionally, we create a validation dataset, denoted as Type 4 dataset, very
different from the training dataset. It considers a much longer bar
($\Omega=[-133.3,133.3]$), under a different loading condition from the
training dataset, and with a $50$ times longer simulation time
($t\in[0,100]$). Therefore, the cross-validation error checks the
generalizability of the estimators.

Results assessment. We present the learnt estimators in Figure 3. Since there
is no ground-truth kernel, we assess the performance of each estimator based
on three criteria. Firstly, we report in Figure 3(b) the prediction $L^{2}$
error of displacement on the cross validation dataset at $T=100$. Secondly, we
report in Figure 3(d) the resultant estimators the group velocity curves from
our model and compare them with the curves computed with DNS. These curves
directly depicts how much our surrogate model reproduces the dispersion
properties in the heterogeneous material. At last, the learnt model should
provide a physically stable material model. To check this, we also report the
dispersion curve in 3(e). Its positivity indicates that the learnt nonlocal
model is physically stable.

Performance of the estimators. Comparing the three estimators in Figure 3(c),
one can see that only the SIDA-RKHS regularizer obtains consistent estimators
in all three experiment settings. The oscillatory estimators of regularizers
with $l^{2}$ or $L^{2}$-norm verify the ill-posedness, and highlight the
importance of using proper regularizers in nonlocal operator learning methods.
The dispersion curves in Figure 3(e) stress the importance of regularizer from
another aspect of view: our SIDA-RKHS regularizer provides physically stable
material models in all settings, while the regularizers with $l^{2}$ or
$L^{2}$-norm may result in highly oscillatory and non-physical models.

We further examine the regularized estimator in terms of its capability in
reproducing DNS simulations through the prediction error of $u$ on the cross
validation dataset. When $\Delta x=0.025$, it takes about $48$ hours for the
DNS simulation to generate one sample, while the homogenized nonlocal model
only requires less than $20$ minutes. From 3(b), we can see that when $\Delta
x=0.05$, all three regularizers are robust and able to reproduce the DNS
simulation with a reasonable accuracy ($\sim 20\%$). When we increase the data
resolution to $\Delta x=0.025$, the estimated nonlocal model from $l^{2}$
regularizer becomes unstable, which again verifies our analysis: when the data
mesh refines, the kernel learning problem becomes more ill-posed and a good
regularizer becomes a necessity. Meanwhile, both the $L^{2}$ and SIDA-RKHS
regularizers lead to a more accurate estimator, indicating a trend of
convergence. On both datasets, the SIDA-RKHS regularizer obtains the most
accurate estimators.

###  3.3 Limitations and future directions

_Non-radial high-dimensional kernels._ When the kernel is radial, our
algorithm is readily applicable to higher dimensions (see Appendix B). When
the kernel is non-radial high-dimensional, however, the regression will face
the well-known curse-of-dimensionality, but our identifiability theory remains
valid. Thus, a future direction is to utilize methods such as kernel-
regression or neural networks and further develop the SIDA-RKHS
regularization.
_Convergence analysis._ We have obtained convergent regularized estimators,
but a convergence analysis is left as future work. The main difficulty to
overcome is the complex combination of three factors: operator spectrum decay,
the errors from numerical integration and noise, and regularization.

##  4 Conclusion

We have characterized the identifiability pitfall in the learning of kernels
in nonlocal operators, and proposed a new regularization method to fix this
issue and achieve estimator convergence. In particular, we have established a
rigorous identifiability theory for the nonparametric learning of kernels in
nonlocal diffusion operators, specifying the function space of
identifiability. Based on the theory, we have introduced a nonparametric
regression algorithm equipped with a data adaptive RKHS regularization method.
Tests on synthetic and real-world datasets show that the our algorithm
consistently obtains accurate and convergent estimator, outperforming common
benchmark regularizers. Our method addresses the critical estimator diverging
phenomena observed in previous nonlocal operator learning methods, and the
proposed framework provides a promising new direction towards overcoming the
ill-posedness to achieve convergence in operator learning.

##  Appendix A Proofs

###### Proof of Lemma 2.3.

Part (a) follows directly from the definition of ${\overline{G}}$. Recall that
a bivariate function ${\overline{G}}$ is positive semi-definite iff for any
$(c_{1},\ldots,c_{m})\in\mathbb{R}^{m}$ and any
$\\{r_{j}\\}_{j=1}^{m}\subset\mathbb{R}^{d}$, the sum
$\sum_{k=1}^{m}\sum_{j=1}^{m}c_{k}c_{j}{\overline{G}}(r_{k},r_{j})\geq 0$ (see
e.g. [4, 10, 29]). Then, noting that from (11) and (10) we have

$\displaystyle\sum_{k=1}^{m}\sum_{j=1}^{m}c_{k}c_{j}{\overline{G}}(r_{k},r_{j})$
$\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\int_{|\eta|=1}\int_{|\xi|=1}\left[\int\sum_{k=1}^{m}\sum_{j=1}^{m}c_{k}c_{j}\frac{[u_{i}(x+r_{k}\xi)-u_{i}(x)][u_{i}(x+r_{j}\eta)-u_{i}(x)]}{\rho_{N}^{\prime}(r_{j})\rho_{N}^{\prime}(r_{k})}dx\right]d\xi
d\eta$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\int_{|\eta|=1}\int_{|\xi|=1}\left[\int\left|\sum_{k=1}^{m}c_{k}\frac{[u_{i}(x+r_{k}\xi)-u_{i}(x)]}{\rho_{N}^{\prime}(r_{k})}dx\right|^{2}\right]d\xi
d\eta\geq 0.$

Thus, ${\overline{G}}$ is positive semi-definite.

For Part (b), the operator ${\mathcal{L}_{\overline{G}}}$ is compact because
${\overline{G}}\in L^{2}(\rho_{N}\times\rho_{N})$, which follows from the fact
that each $u_{i}$ is bounded (thus, ${\overline{G}}$ is also bounded). Also,
since ${\overline{G}}$ is positive semi-definite, so is
${\mathcal{L}_{\overline{G}}}$. The equation (12) follows from (16).

Part (c) is a standard operator characterization of the RKHS $H_{G}$ (see
e.g., [10]).

For Part (d), the eigenfunctions are orthonormal and the eigenvalues decay to
zero because the operator ${\mathcal{L}_{\overline{G}}}$ is positive semi-
definite and compact, as shown in Part (b). The first equation in (13) follows
from (12), and the second equation follows from the orthogonality of the
eigenfunctions. At last, if $\phi\in H_{G}$, by the characterization of
$H_{G}$’s inner product in Part (c), we have the third equation in (13). ∎

###### Proof of Lemma 2.2.

Recall that with the bilinear form
$\langle\langle{\cdot,\cdot}\rangle\rangle$, defined in (5), we can rewrite
the loss functional as

$\displaystyle\mathcal{E}(\phi)$
$\displaystyle=\langle\langle{\phi,\phi}\rangle\rangle-\frac{1}{N}\sum_{i=1}^{N}\int
2L_{\phi}[u_{i}](x)f_{i}(x)dx+C_{f},$ (15)

where $C_{N}^{f}=\frac{1}{N}\sum_{k=1}^{N}\int|f_{i}(x)|^{2}dx$. Then, the
derivative $\nabla\mathcal{E}(\phi)$ follows from (12) and a rewriting of the
bilinear form:

$\displaystyle\langle\langle{\phi_{1},\phi_{2}}\rangle\rangle=\frac{1}{N}\sum_{i=1}^{N}\int\left[\int\int\phi_{1}(|z|)[u_{i}(x+z)-u_{i}(x)]\phi_{2}(|y|)[u_{i}(x+y)-u_{i}(x)dydz]\right]dx$
$\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\int\int\phi_{1}(|z|)\phi_{2}(|y|)\left[\int[u_{i}(x+z)-u_{i}(x)][u_{i}(x+y)-u_{i}(x)]dx\right]dydz$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\int_{0}^{\infty}\phi_{1}(r)\phi_{2}(s)G(r,s)drds=\int_{0}^{\infty}\int_{0}^{\infty}\phi_{1}(r)\phi_{2}(s){\overline{G}}(r,s)\rho_{N}(dr)\rho_{N}(ds),$
(16)

with $G$ and ${\overline{G}}$ given in (11) and (10), where the last equality
is a re-weighting by $\rho_{N}$. ∎

###### Proof of Theorem 2.1 .

By Lemma 2.2, the Fréchet derivative of the loss functional is
$\nabla\mathcal{E}(\phi)=2({\mathcal{L}_{\overline{G}}}\phi-\phi_{N}^{f})$.
Thus, the loss functional has a unique minimizer only in the function space
where $\nabla\mathcal{E}(\phi)$ has a unique zero, that is, the operator
${\mathcal{L}_{\overline{G}}}$ has an inversion. The largest such a space is
the eigenspace expanded by all eigenfunctions with non-zero eigenvalues of
${\mathcal{L}_{\overline{G}}}$. Furthermore, projecting $\phi_{N}^{f}$ to this
sapce, we have the the minimizer
$\widehat{\phi}={\mathcal{L}_{\overline{G}}}^{-1}P\phi_{N}^{f}$ as given in
the theorem. ∎

##  Appendix B Algorithm: nonparametric regression with SIDA-RKHS
regularization

In this section we provide detailed description of the algorithm proposed in
Section 2.3.

Our algorithm consists of three steps. First, we utilize the data to estimate
the exploration measure and the support of the kernel. Based on them, we set a
class of hypothesis spaces, with their dimensions i.e., the number of basis
functions, in a proper range moving from under-fitting to over-fitting.
Second, we assemble the regression matrices vectors from data for each of the
hypothesis spaces. Finally, we identify the estimators with SIDA-RKHS
regularization for these hypothesis spaces and select the one with the best
fitting.

To start, we assume that the discrete data
$\\{u_{i}(x_{j}),f_{i}(x_{j})\\}_{i=1}^{N}$ comes with equidistant mesh points
$\\{x_{j}=j\Delta x\\}_{j=0}^{J}$. For simplicity, we consider only the 1D
case, and the extension to multi-dimensional cases is straightforward. We note
that the current problem setting assumes data on mesh points, thus the data
size increases exponentially as the dimension increases, which is the well-
known curse-of-dimensionality. To overcome this curse-of-dimensionality, one
can consider other settings with mesh-free representation of data by random
samples and a loss functional based on expectations (see, e.g. [27]), and this
is beyond the scope of the current study.

#### Step 1: Set a class of hypothesis spaces.

We set a class of data-adaptive hypothesis spaces with their dimensions set to
range from under-fitting to over-fitting. The key is the exploration measure
and the support of the kernel estimated data. The exploration measure
$\rho_{N}^{J}$ is computed from data as in (8), which uses only the
information from $u_{i}$. To estimate the support of the kernel, we extract
the additional information from $\\{f_{i}\\}$ as follows. We set the data-
adaptive support of the kernel to be $[0,R]$ with $R$ defined by

$R=1.1\min\\{R_{\rho},\max\\{|L^{f}_{i}-L^{u}_{i}|,|R^{f}_{i}-R^{u}_{i}|\\}_{i=1}^{N}\\},$
(17)

where $(L^{u}_{i},R^{u}_{i})$ and $(L^{f}_{i},R^{f}_{i})$ are the lower and
upper bounds of the supports $\mathrm{supp}(u_{i})$ and $\mathrm{supp}(f_{i})$
respectively, and $R_{\rho}$ is the maximum of the support of $\rho_{N}^{J}$.
That is, the support of the kernel lies inside the support of the exploration
measure, and it is the maximal interaction range indicated by the difference
between supports of $u_{i}$ and $f_{i}$, which extracts the additional
information in the data $\\{f_{i}\\}$. Here the multiplicative factor 1.1 is
an artificial factor to enlarge the range, so that the supports of the basis
functions will fully cover the explored region. To avoid unbounded support in
the data-based estimation in (17), in numerical experiments we set a threshold
to be $10^{-8}$ when estimating supports of $u_{i}$, $f_{i}$ and
$\rho_{N}^{J}$. This truncation narrows the interaction range.

The estimated support of the kernel is the region explored by data. Outside of
the region, the data provides little information about the kernel. Thus, we
focus on learning the kernel in this region and set the local basis functions
to be supported in it. Furthermore, we constrain the exploration measure to be
supported in $[0,R]$. For simplicity of notation, we still denote it by
$\rho_{N}^{J}$ or $\rho_{N}$.

With the exploration measure and the support of the kernel, we select a class
of basis functions $\\{\phi_{i}\\}_{i=1}^{n}$ and a range of $n$ for the
hypothesis space $\mathcal{H}_{n}=\mathrm{span}\\{\phi_{i}\\}_{i=1}^{n}$. The
basis function can be either global basis functions such as Bernstein
polynomials as those used in [52, 54] and trigonometric functions, or local
basis functions such B-spline polynomials (see Appendix C.1 for a brief
introduction). We focus on local basis functions because they are more
flexible to adaptive to local structure of the kernel. To set the range for
$n$, we note that the mesh points of the kernel’s independent variable
explored by data are $\\{k\Delta x:k=1,\ldots,\lfloor{\frac{R}{\Delta
x}}\rfloor\\}$. Meanwhile, the basis function should be linearly independent
in $L^{2}(\rho_{N}^{J})$ so that the basis matrix

$B_{n}=(\langle{\phi_{i},\phi_{j}}\rangle_{L^{2}(\rho_{N}^{J})})_{1\leq
i,j\leq n}\in\mathbb{R}^{n\times n}$ (18)

is non-singular. Thus, we set the range of $n$ to be in
$\lfloor{\frac{R}{\Delta x}}\rfloor\times[0.2,1]$ such that $B_{n}$ is non-
singular while covering a wide range of dimensions. For example, when we use
piecewise constant basis, we can set $n=\lfloor{\frac{R}{\Delta x}}\rfloor$,
and we get $B_{n}=\mathrm{Diag}(\rho_{N}^{J})$. Thus, we estimate the kernel
as a vector of its values on the mesh points, with $L^{2}(\rho_{N}^{J})$ being
a vector space with a discrete-measure $\rho_{N}^{J}$.

#### Step 2: Assemble regression matrices and vectors.

We assemble the regression matrix ${\overline{A}}_{n}$ and vector
${\overline{b}}_{n}$, as defined in (4), for each hypothesis spaces
$\mathcal{H}_{n}=\mathrm{span}\\{\phi_{i}\\}_{i=1}^{n}$. Together with the
basis matrix $B_{n}$ in (18), the triplet
$({\overline{A}}_{n},{\overline{b}}_{n},B_{n})$ is all we need for regression
with SIDA-RKHS regularization in the next step.

To avoid repeated reading of data, we extract the regression data that can be
used for all hypothesis spaces by utilizing the regression structure, which
requires reading the data only once. Note that to compute
${\overline{A}}_{n}(i,j)=\langle\langle{\phi_{i},\phi_{j}}\rangle\rangle$ for
any pair of basis functions, with the bilinear form defined in (16), we only
need $G$ defined in (11). We note that when $d=1$, the integral
$\int_{|\eta|=1}g(\eta)d\eta=g(\eta)+g(-\eta)$, therefore, we have

$\displaystyle
G(r,s)=\frac{1}{N}\sum_{i=1}^{N}\int[u_{i}(x+r)+u_{i}(x-r)-2u_{i}(x)][u_{i}(x+s)+u_{i}(x-s)-2u_{i}(x)]dx$
(19)

for $r,s\in\mathrm{supp}(\rho_{N})$. Similarly, for a basis function
$\phi_{i}$, to compute ${\overline{b}}(i)$ in (4), which can be re-written as
${\overline{b}}_{n}(i)=\frac{1}{N}\sum_{k=1}^{N}\int
L_{\phi_{i}}[u_{k}](x)f_{k}(x)dx=\int_{0}^{R}\phi_{i}(r)g_{N}^{f}(r)dr$, we
only need the function $g_{N}^{f}$ defined by

$g_{N}^{f}(r)=\frac{1}{N}\sum_{i=1}^{N}\int_{\Omega}\int_{|\xi|=1}[u_{i}(x+r\xi)-u_{i}(x)]f_{i}(x)d\xi\,dx.$
(20)

Let $r_{k}=k\Delta x$ for $k=1,\ldots,\lfloor{\frac{R}{\Delta x}}\rfloor$,
which are all the mesh points the data explore. Then, all the regression data
we need in the original data (1) are

$\left\\{G(r_{k},r_{l}),g_{N}^{f}(r_{k}),\rho_{N}^{J}(r_{k}),\text{ with
}k,l=1,\ldots,\lfloor{\frac{R}{\Delta x}}\rfloor\right\\},$ (21)

where $G$, $g_{N}^{f}$ and $\rho_{N}^{J}$ are defined respectively in (11),
(20) and (8).

With these regression data, the triplet
$({\overline{A}}_{n},{\overline{b}}_{n},B_{n})$ can be efficiently evaluated
for any basis functions using a numerical integrator to approximate the
corresponding integrals. For example, with Riemann sum approximation, we
compute the normal matrix ${\overline{A}}_{n}$ and vector ${\overline{b}}_{n}$
and the basis matrix $B_{n}$ as

$\displaystyle{\overline{A}}_{n}(i,j)$
$\displaystyle=\langle\langle{\phi_{i},\phi_{j}}\rangle\rangle\approx\sum_{k,l}\phi_{i}(r_{k})\phi_{j}(r_{l})G(r_{k},r_{l}))\Delta
x^{2},$ (22) $\displaystyle{\overline{b}}_{n}(i)$
$\displaystyle\approx\sum_{k}\phi_{i}(r_{k})g_{N}^{f}(r_{k}))\Delta x,$
$\displaystyle B_{n}(i,j)$
$\displaystyle\approx\sum_{k}\phi_{i}(r_{k})\phi_{j}(r_{k})\rho_{N}^{J}(r_{k})\Delta
x.$

#### Step 3: Regress with SIDA-RKHS regularization.

Our SIDA-RKHS regularization method uses the norm of the SIDA-RKHS so as to
ensure the learning to take space in the function space of identifiability as
discussed in Section 2.2. That is, our estimator is the minimizer of the
regularized loss in (7) with the regularization norm
$\mathcal{R}(\phi)=\|\phi\|^{2}_{H_{G}}$ defined in (13).

_Computation of the RKHS norm_. We can effectively approximate the RKHS norm
$\|\phi\|^{2}_{H_{G}}$ using the triplet
$({\overline{A}}_{n},{\overline{b}}_{n},B_{n})$. It proceeds in two steps.
First, we solve the generalized eigenvalue problem
${\overline{A}}_{n}V=B_{n}V\Lambda$, where $\Lambda$ is a diagonal matrix of
the generalized eigenvalues and the matrix $V$ has columns being eigenvectors
orthonormal in the sense that $V^{\top}B_{n}V=I_{n}$. Here these eigenvalues
approximate the eigenvalue of ${\mathcal{L}_{\overline{G}}}$ in (9), and
$\widehat{\psi}_{k}=V_{jk}\phi_{j}$ approximates the eigenfunctions of
${\mathcal{L}_{\overline{G}}}$. Then, we compute the square RKHS norm of
$\phi=\sum_{i}c_{i}\phi_{i}$ as

$\|\phi\|^{2}_{H_{G}}=c^{\top}B_{rkhs}c,\,\text{ with }B_{rkhs}=(V\Lambda
V^{\top})^{-1},$ (23)

where the inverse is taken as pseudo-inverse, particularly when $\Lambda$ has
zero eigenvalues.

With the RKHS-norm ready, we write the regularized loss for each function
$\phi=\sum_{i}c_{i}\phi_{i}$ as
$\mathcal{E}_{\lambda}(\phi)=c^{\top}({\overline{A}}_{n}+\lambda
B_{rkhs})c-2c^{\top}{\overline{b}}_{n}+C_{N}^{f}.$ The regularized estimator
is

$\displaystyle\widehat{\phi_{\lambda}}=\sum_{i=1}^{n}c^{i}_{\lambda}\phi_{i},\
c_{\lambda}=({\overline{A}}_{n}+\lambda B_{rkhs})^{-1}{\overline{b}}_{n}.$
(24)

We will select the hyper-parameter that balances the loss $\mathcal{E}$ and
the regularization term by the widely-used L-curve method [17]. It identifies
the optimal hyper-parameter as the maximizer of the curvature of the curve
(see Section C.2).

##  Appendix C B-spline basis functions and the L-curve method

###  C.1 B-spline basis functions

B-spline is a class of piecewise polynomials, and is capable of representing
the local information of the target function. Here we review briefly the
recurrence definition and properties of the balanced B-splines, for more
details we refer to the Chapter 2 of [42] and [40].

Given a non-decreasing sequence of real numbers
$\\{r_{0},r_{1},\ldots,r_{m}\\}$ (called knots), the B-spline basis functions
of degree $p$, denoted by $\\{N_{i,p}\\}_{i=0}^{m-p}$, is defined recursively
as

$\displaystyle N_{i,0}(r)=\left\\{\begin{array}[]{lr}1,&r_{i}\leq
r<r_{i+1},\\\ 0,&otherwise,\end{array}\right.$ (25) $\displaystyle
N_{i,p}(r)=\frac{r-r_{i}}{r_{i+p}-r_{i}}N_{i,p-1}(r)+\frac{r_{i+p+1}-r}{r_{i+p+1}-r_{i+1}}N_{i+1,p-1}(r).$

The B-spline basis has the following properties:

  * •

Each function $N_{i,p}$ is a nonnegative local polynomial of degree $p$,
supported on $[r_{i},r_{i+p+1}]$;

  * •

At a knot with multiplicity $k$, it is $p-k$ times continuously
differentiable. Hence, the smoothness increases with the degree but decreases
when the knot multiplicity increases;

  * •

The basis satisfies partition unity: for each $r\in[r_{i},r_{i+1}]$,
$\sum_{j}N_{j,p}(r)=\sum_{j=i-p}^{i}N_{j,p}(r)=1$.

We set the knots to be a uniform partition of the support of
$\overline{\rho}$, $[R_{min},R_{max}]$,

$R_{min}=r_{0}\leq r_{1}\leq\cdots\leq r_{m}=R_{min}.$

We set the basis functions of the hypothesis $\mathcal{H}$, whose dimension is
$n=m-p$, to be

$\phi_{i}(r)=N_{i,p}(r),\ i=1,\dots,m-p.$

Thus, the basis functions $\\{\phi_{i}\\}$ are piecewise degree-$p$
polynomials with knots adaptive to $\overline{\rho}$.

###  C.2 Hyper-parameter selection by the L-curve method

We select the parameter $\lambda$ by the L-curve method [17, 27]. Let $l$ be a
parametrized curve in $\mathbb{R}^{2}$:

$l(\lambda)=(x(\lambda),y(\lambda)):=(\text{log}(\mathcal{E}(\widehat{\phi_{\lambda}}),\text{log}(\mathcal{R}(\widehat{\phi_{\lambda}})),$

where
$\mathcal{E}(\widehat{\phi_{\lambda}})=c_{\lambda}^{\top}{\overline{A}}_{n}c_{\lambda}-2c_{\lambda}^{\top}{\overline{b}}_{n}-C_{N}^{f}$,
and $\mathcal{R}(\phi)$ is the regularization term, for example,
$\mathcal{R}(\widehat{\phi_{\lambda}})=\|\widehat{\phi_{\lambda}}\|_{H_{\overline{G}}}^{2}=c_{\lambda}^{\top}B_{rkhs}c_{\lambda}$.
The optimal parameter is the maximizer of the curvature of $l$. In practice,
we restrict $\lambda$ in the spectral range $[\lambda_{min},\lambda_{max}]$ of
the operator ${\mathcal{L}_{\overline{G}}}$,

$\displaystyle\lambda_{0}=\underset{\lambda_{\text{min}}\leq\lambda\leq\lambda_{\text{max}}}{\operatorname{arg}\operatorname{max}}\;\kappa(l(\lambda))=\underset{\lambda_{\text{min}}\leq\lambda\leq\lambda_{\text{max}}}{\operatorname{arg}\operatorname{max}}\;\frac{x^{\prime}y^{\prime\prime}-x^{\prime}y^{\prime\prime}}{(x^{\prime}\,{}^{2}+y^{\prime}\,{}^{2})^{3/2}},$
(26)

where $\lambda_{min}$ and $\lambda_{max}$ are computed from the smallest and
the largest generalized eigenvalues of $({\overline{A}}_{n},B_{n})$. This
optimal parameter $\lambda_{0}$ balances the loss $\mathcal{E}$ and the
regularization (see [17] for more details). In practice, instead of computing
the second order derivatives, we compute the curvature by the reciprocal of
the radius of the interior circle of three consecutive points111Are Mjaavatten
(2022). Curvature of a 1D curve in a 2D or 3D space
(https://www.mathworks.com/matlabcentral/fileexchange/69452-curvature-
of-a-1d-curve-in-a-2d-or-3d-space), MATLAB Central File Exchange..

##  Appendix D Additional numerical results for synthetic data examples

This section provides additional numerical results for the examples with
synthetic data.

Figure 4 shows that the SIDA-RKHS regularizer leads to converging estimators
in all three examples for both noisy and noiseless data, whereas the
$l^{2}$-norm and the $L^{2}$-norm regularizers’ estimators have slow
convergent rates or even no convergence when the data is noisy.

Figure 4: Convergence of function estimators as the data mesh-size $\Delta x$
refines, along with values of the loss function. The SIDA-RKHS regularizer
consistently converges for both noiseless and noisy data, with better rates
(slope) than the other two regularizers for noisy data. Note that for the
fractional kernel, it has a lower rate though being more accurate.

We note that the performance of these regularizers depends on the optimal
regularization strength $\lambda_{0}$, which is selected by the L-curve method
introduced in Section C.2. In our tests, all regularizers can successfully
select the optimal $\lambda_{0}$ for most of the time, and the SIDA-RKHS
regularizer has the most well-shaped L-curve, which leads to the most robust
regularization (see Figure 5 for typical L-curve plots).

Figure 5: Typical L-curve plots for the selection of the optimal
regularization parameter $\lambda_{0}$ for the Gaussian kernel with $\Delta
x=0.05$ and $nsr=1$. From left to right: the $l^{2}$, $L^{2}$ and SIDA-RKHS
regularizers. All regularizers successfully select the optimal $\lambda_{0}$,
and the SIDA-RKHS regularizer has the most well-shaped L-curve.

##  Appendix E Detailed Real-world Dataset Experiment Settings

In this section we provide further experiment details for the real-world
dataset studied in 3.2.

For both training and validation purposes we generate data using high-fidelity
(HF) simulations for the propagation of stress waves within the microstructure
of the heterogeneous, linear elastic bar. In the following, we use $\hat{u}$
to denote the HF solution, to distinguish the HF dataset from the homogenized
solution of (14). The HF-model is a classical wave equation: the displacement
$\hat{u}(x,t)$ satisfies, for $(x,t)\in\Omega\times[0,T]$ with
$\Omega\subset\mathbb{R}$,

$\displaystyle\partial_{tt}\hat{u}(x,t)-L_{HF}[\hat{u}](x,t)=g(x,t),$ (27)

with a force loading term $g(x,t)$, proper boundary conditions and initial
conditions $\hat{u}(x,0)=0$, $\partial_{t}\hat{u}(x,0)=0$. Considering the
heterogeneous bar of two materials depicted in Figure 3, (27) describes the
stress wave propagating with speed $c_{1}=\sqrt{E_{1}/\rho}$ in material 1 and
speed $c_{2}=\sqrt{E_{2}/\rho}$ in material 2. We solve the HF-model (27) by
the direct numerical solver (DNS) introduced in [44]. The DNS employs the
characteristic line method, which provides exact solutions of velocities. For
each grid point $x_{j}\in\Omega$ at time step $t^{n}=n\Delta t$, where $\Delta
t$ is the time step size, with the calculated exact velocity
$\hat{v}(x_{j},t^{n})$ and the estimated displacement from the last time step
$\hat{u}(x_{j},t^{n-1})$ we update the HF displacement by

$\hat{u}(x_{j},t^{n})=\hat{u}(x_{j},t^{n-1})+\Delta t\hat{v}(x_{j},t^{n}).$

With the above procedure, we then consider various boundary velocity loading
$\partial_{t}\hat{u}_{i}(x,t)$, $x\in\partial\Omega$, and force loading
$g_{i}(x,t)$ scenarios, and solve for the corresponding HF displacement field
$\hat{u}_{i}(x,t)$. Resultant data pairs
$\\{\hat{u}_{i},g_{i}\\}_{i=1}^{N}=\\{\hat{u}_{i}(x_{j},t^{n}),g_{i}(x_{j},t^{n}):j=1,\ldots,J\\}_{i=1,n=0}^{N,T/\Delta
t}$ are employed as the training and validation datasets. Discretization
parameters for the DNS solver are set to $\Delta t=0.01$ and $\max{\Delta
x}=0.01$.

The homogenization problem is then to learn the kernel of the nonlocal
operator $L_{\phi}$ that approximates the operator $L_{HF}$ from data
$\\{\hat{u},f\\}$ generated by $L_{HF}[\hat{u}]=f$, where
$f=\partial_{tt}\hat{u}-g$. Discretizing the time derivative in (14) with the
central difference scheme, we obtain

$\dfrac{1}{\Delta
t^{2}}(\hat{u}^{n+1}(x)-2\hat{u}^{n}(x)+\hat{u}^{n-1}(x))-g(x,t^{n}):=f^{n}(x),$

where $\hat{u}^{n}(\cdot):=\hat{u}(\cdot,t^{n})$ denotes the solution at time
$t^{n}$. Given
$\mathcal{D}=\\{\hat{u}_{i}^{n}(x),f_{i}^{n}(x)\\}_{i=1,n=1}^{N,T/\Delta t}$,
our goal is to learn the kernel $\phi$. The loss functional is

$\displaystyle\mathcal{E}(\phi)$ $\displaystyle=\frac{\Delta
t}{NT}\sum_{k=1}^{N}\sum_{n=1}^{T/\Delta
t}\|L_{\phi}[\hat{u}^{n}_{k}]-f^{n}_{k}\|^{2}_{L^{2}(\Omega)}.$ (28)

###  E.1 Settings on real-world data

In the learning problem, we consider four types of data and use the first
three for training and the last one for validation of our algorithm. For all
data we set $L=0.2$, $\Delta t=0.02$, and the symmetric domain
$\Omega=[-b,b]$. Two spatial resolutions, $\Delta=0.05$ and $\Delta=0.025$ are
considered, which we denote as the “coarse” and “fine” datasets, respectively.

  1. Type 1

Oscillating source (20 samples in total). $b=50$, $T=2$,
$g(x,t)=\exp^{-(\frac{2x}{5jL})^{2}}\exp^{-(\frac{t-0.8}{0.8})^{2}}\cos^{2}(\frac{2\pi
x}{jL})$, where $j=1,2,\cdots,20$.

  2. Type 2

Plane wave with $\cos$ loading (11 samples in total). $b=50$, $T=2$,
$g(x,t)=0$ and $\partial_{t}{u}(-50,t)=\cos(jt)$, where the loading frequency
$j=0.35,0.70,\cdots,3.85$.

  3. Type 3

Plane wave with $\sin$ loading (11 samples in total). $b=50$, $T=2$,
$g(x,t)=0$ and $\partial_{t}{u}(-50,t)=\sin(jt)$, where the loading frequency
$j=0.35,0.70,\cdots,3.85$.

  4. Type 4

Wave packet (3 samples in total). $b=133.3$, $T=100$, $g(x,t)=0$ and
$\partial_{t}{u}(-b,t)=\sin(jt)\exp\left(-(t/5-3)^{2}\right)$, for
$j=1,\,2,\,3$.

Notice that the validation dataset (Type 4 dataset) is under a different
loading condition from the training dataset, and with a much longer simulation
time.

## Acknowledgments and Disclosure of Funding

YY are supported by the National Science Foundation under award DMS 1753031,
and the AFOSR grant FA9550-22-1-0197. YY would also like to like to thank Dr.
Stewart Silling for his help on the DNS codes and for valuable discussions. FL
is grateful for supports from NSF-1913243 and FA9550-20-1-0288. FL and QA
would like to thank Quanjun Lang for helpful discussions on regularization.

## References

  * [1] Fuensanta Andreu-Vaillo, José Mazón, Julio Rossi, and J. Julián Toledo-Melero.  Nonlocal Diffusion Problems, volume 165 of Mathematical Surveys and Monographs.  American Mathematical Society, Providence, Rhode Island, 2010.
  * [2] David Applebaum.  Lévy processes and stochastic calculus.  Cambridge university press, 2009.
  * [3] Frank Bauer, Sergei Pereverzev, and Lorenzo Rosasco.  On regularization algorithms in learning theory.  Journal of complexity, 23(1):52–72, 2007.
  * [4] Christian Berg, Jens Peter Reus Christensen, and Paul Ressel.  Harmonic analysis on semigroups: theory of positive definite and related functions, volume 100.  New York: Springer, 1984.
  * [5] Claudia Bucur and Enrico Valdinoci.  Nonlocal Diffusion and Applications, volume 20 of Lecture Notes of the Unione Matematica Italiana.  Springer International Publishing, Cham, 2016.
  * [6] Olena Burkovska, Christian Glusa, and Marta D’Elia.  An optimization-based approach to parameter learning for fractional type nonlocal models.  Computers & Mathematics with Applications, 2021.
  * [7] Yifan Chen, Houman Owhadi, and Andrew M. Stuart.  Consistency of Empirical Bayes And Kernel Flow For Hierarchical Parameter Estimation.  ArXiv200511375 Cs Math Stat, 2021.
  * [8] Zhen-Qing Chen and Xicheng Zhang.  Heat Kernels for Non-symmetric Non-local Operators.  In Giampiero Palatucci and Tuomo Kuusi, editors, Recent Developments in Nonlocal Theory, pages 24–51. De Gruyter Open, 2017\.
  * [9] Felipe Cucker and Steve Smale.  On the mathematical foundations of learning.  Bulletin of the American mathematical society, 39(1):1–49, 2002\.
  * [10] Felipe Cucker and Steve Smale.  Emergent behavior in flocks.  IEEE Transactions on Automatic Control, 52(5):852–862, 2007.
  * [11] Marta D’Elia, Qiang Du, Christian Glusa, Max Gunzburger, Xiaochuan Tian, and Zhi Zhou.  Numerical methods for nonlocal and fractional models.  Acta Numerica, 29:1–124, 2020.
  * [12] Qiang Du, Max Gunzburger, R. B. Lehoucq, and Kun Zhou.  Analysis and Approximation of Nonlocal Diffusion Problems with Volume Constraints.  SIAM Rev., 54(4):667–696, 2012.
  * [13] Yiming Fan, Xiaochuan Tian, Xiu Yang, Xingjie Li, Clayton Webster, and Yue Yu.  An asymptotically compatible probabilistic collocation method for randomly heterogeneous nonlocal problems.  arXiv preprint arXiv:2107.01386, 2021.
  * [14] Frédéric Ferraty and Philippe Vieu.  Nonparametric functional data analysis: theory and practice, volume 76.  Springer, 2006.
  * [15] Guy Gilboa and Stanley Osher.  Nonlocal operators with applications to image processing.  Multiscale Modeling & Simulation, 7(3):1005–1028, 2009.
  * [16] Boumediene Hamzi and Houman Owhadi.  Learning dynamical systems from data: A simple cross-validation perspective, part I: Parametric kernel flows.  Physica D: Nonlinear Phenomena, 421:132817, 2021.
  * [17] P. C. Hansen.  The L-Curve and its Use in the Numerical Treatment of Inverse Problems.  WIT Press, 2000.
  * [18] Per Christian Hansen.  REGULARIZATION TOOLS: A Matlab package for analysis and solution of discrete ill-posed problems.  Numer Algor, 6(1):1–35, 1994.
  * [19] Gernot Holler and Karl Kunisch.  Learning nonlocal regularization operators.  ArXiv200109092 Math, 2020.
  * [20] Tailen Hsing and Randall Eubank.  Theoretical foundations of functional data analysis, with an introduction to linear operators, volume 997.  John Wiley & Sons, 2015.
  * [21] Rakhoon Hwang, Jae Yong Lee, Jin Young Shin, and Hyung Ju Hwang.  Solving pde-constrained control problems using operator learning.  arXiv preprint arXiv:2111.04941, 2021.
  * [22] Bangti Jin and William Rundell.  A tutorial on inverse problems for anomalous diffusion processes.  Inverse problems, 31(3):035003, 2015.
  * [23] Hachem Kadri, Emmanuel Duflos, Philippe Preux, Stéphane Canu, Alain Rakotomamonjy, and Julien Audiffren.  Operator-valued kernels for learning from functional response data.  Journal of Machine Learning Research, 17(20):1–54, 2016.
  * [24] Stefan Kindermann, Stanley Osher, and Peter W Jones.  Deblurring and denoising of images by nonlocal functionals.  Multiscale Modeling & Simulation, 4(4):1091–1115, 2005.
  * [25] Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Neural operator: Learning maps between function spaces.  arXiv preprint arXiv:2108.08481, 2021.
  * [26] Quanjun Lang and Fei Lu.  Identifiability of interaction kernels in mean-field equations of interacting particles.  arXiv preprint arXiv:2106.05565, 2021.
  * [27] Quanjun Lang and Fei Lu.  Learning interaction kernels in mean-field equations of first-order systems of interacting particles.  SIAM Journal on Scientific Computing, 44(1):A260–A285, 2022.
  * [28] Yang Li, Yubin Lu, Shengyuan Xu, and Jinqiao Duan.  Extracting stochastic dynamical systems with $alpha$-stable lévy noise from data.  arXiv preprint arXiv:2109.14881, 2021.
  * [29] Zhongyang Li, Fei Lu, Mauro Maggioni, Sui Tang, and Cheng Zhang.  On the identifiability of interaction functions in systems of interacting particles.  Stochastic Processes and their Applications, 132:135–163, 2021\.
  * [30] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Neural operator: Graph kernel network for partial differential equations.  arXiv preprint arXiv:2003.03485, 2020.
  * [31] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Neural Operator: Graph Kernel Network for Partial Differential Equations.  ArXiv200303485 Cs Math Stat, 2020.
  * [32] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Andrew Stuart, Kaushik Bhattacharya, and Anima Anandkumar.  Multipole graph neural operator for parametric partial differential equations.  Advances in Neural Information Processing Systems, 33, 2020.
  * [33] Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, Anima Anandkumar, et al.  Fourier neural operator for parametric partial differential equations.  In International Conference on Learning Representations, 2020.
  * [34] Chensen Lin, Zhen Li, Lu Lu, Shengze Cai, Martin Maxey, and George Em Karniadakis.  Operator learning for predicting multiscale bubble growth dynamics.  The Journal of Chemical Physics, 154(10):104118, 2021.
  * [35] Chensen Lin, Martin Maxey, Zhen Li, and George Em Karniadakis.  A seamless multiscale operator neural network for inferring bubble dynamics.  Journal of Fluid Mechanics, 929, 2021.
  * [36] Fei Lu, Mauro Maggioni, and Sui Tang.  Learning interaction kernels in stochastic systems of interacting particles from multiple trajectories.  Foundations of Computational Mathematics, pages 1–55, 2021.
  * [37] Fei Lu, Ming Zhong, Sui Tang, and Mauro Maggioni.  Nonparametric inference of interaction laws in systems of agents from trajectory data.  Proc. Natl. Acad. Sci. USA, 116(29):14424–14433, 2019.
  * [38] Lu Lu, Pengzhan Jin, and George Em Karniadakis.  Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators.  arXiv preprint arXiv:1910.03193, 2019.
  * [39] Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis.  Learning nonlinear operators via deeponet based on the universal approximation theorem of operators.  Nature Machine Intelligence, 3(3):218–229, 2021.
  * [40] Tom Lyche, Carla Manni, and Hendrik Speleers.  Foundations of Spline Theory: B-Splines, Spline Approximation, and Hierarchical Refinement, volume 2219, pages 1–76.  Springer International Publishing, Cham, 2018.
  * [41] Houman Owhadi and Gene Ryan Yoo.  Kernel Flows: From learning kernels from data into the abyss.  Journal of Computational Physics, 389:22–47, 2019.
  * [42] Les Piegl and Wayne Tiller.  The NURBS Book.  Monographs in Visual Communication. Springer Berlin Heidelberg, Berlin, Heidelberg, 1997.
  * [43] Leonid I Rudin, Stanley Osher, and Emad Fatemi.  Nonlinear total variation based noise removal algorithms.  Physica D: nonlinear phenomena, 60(1-4):259–268, 1992.
  * [44] S. A. Silling.  Propagation of a stress pulse in a heterogeneous elastic bar.  Sandia Report SAND2020-8197, Sandia National Laboratories, 2020\.
  * [45] S. A. Silling, M. Epton, O. Weckner, J. Xu, and E. Askari.  Peridynamic States and Constitutive Modeling.  J Elasticity, 88(2):151–184, 2007.
  * [46] Yunzhe Tao, Qi Sun, Qiang Du, and Wei Liu.  Nonlocal neural networks, nonlocal diffusion and nonlocal modeling.  In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pages 494–504, 2018.
  * [47] Robert Tibshirani.  Regression Shrinkage and Selection Via the Lasso.  Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267–288, 1996.
  * [48] Hong Wang and Treena S Basu.  A fast finite difference method for two-dimensional space-fractional diffusion equations.  SIAM Journal on Scientific Computing, 34(5):A2444–A2458, 2012.
  * [49] Xiaolong Wang, Ross Girshick, Abhinav Gupta, and Kaiming He.  Non-local Neural Networks.  In 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 7794–7803, Salt Lake City, UT, USA, 2018. IEEE.
  * [50] Jie Xiong, Jiayu Zheng, and Xiaowen Zhou.  Unique strong solutions of lévy processes driven stochastic differential equations with discontinuous coefficients.  Stochastics, 91(4):592–612, 2019.
  * [51] Huaiqian You, Yue Yu, Marta D’Elia, Tian Gao, and Stewart Silling.  Nonlocal kernel network (nkn): a stable and resolution-independent deep neural network.  arXiv preprint arXiv:2201.02217, 2022.
  * [52] Huaiqian You, Yue Yu, Stewart Silling, and Marta D’Elia.  Data-driven learning of nonlocal models: From high-fidelity simulations to constitutive laws.  ArXiv201204157 Cs Math, 2020.
  * [53] Huaiqian You, Yue Yu, Stewart Silling, and Marta D’Elia.  A data-driven peridynamic continuum model for upscaling molecular dynamics.  Computer Methods in Applied Mechanics and Engineering, 389:114400, 2022.
  * [54] Huaiqian You, Yue Yu, Nathaniel Trask, Mamikon Gulian, and Marta D’Elia.  Data-driven learning of nonlocal physics from high-fidelity synthetic data.  Computer Methods in Applied Mechanics and Engineering, 374:113553, 2021.