# Data adaptive RKHS Tikhonov regularization
for learning kernels in operators

Fei Lu  Quanjun Lang  Qingci An

###### Abstract

We present DARTR: a Data Adaptive RKHS Tikhonov Regularization method for the
linear inverse problem of nonparametric learning of function parameters in
operators. A key ingredient is a system intrinsic data adaptive (SIDA) RKHS,
whose norm restricts the learning to take place in the function space of
identifiability. DARTR utilizes this norm and selects the regularization
parameter by the L-curve method. We illustrate its performance in examples
including integral operators, nonlinear operators and nonlocal operators with
discrete synthetic data. Numerical results show that DARTR leads to an
accurate estimator robust to both numerical error due to discrete data and
noise in data, and the estimator converges at a consistent rate as the data
mesh refines under different levels of noises, outperforming two baseline
regularizers using $l^{2}$ and $L^{2}$ norms.

Key words: ill-posed inverse problem, Tikhonov regularization, RKHS,
identifiability

##  1 Introduction

Regularization plays a crucial role in inverse and machine learning problems
that aim to construct robust generalizable models. The learning of kernel
functions in operators is such a problem: given data
$\\{(u_{k},f_{k})\\}_{k=1}^{N}$ in suitable function spaces, we would like to
learn an optimal kernel function $\phi$ fitting the operator $R_{\phi}(u)=f$
to the data. Such a need for learning operators between function spaces has
become vital in applications ranging from integral operators solving PDEs and
image processing (see e.g., [12, 24, 20, 29]), nonlinear operators in mean-
field equation of interacting particle systems in [26, 22], homogenized
nonlocal operators (see e.g., [39, 38, 25]), just to name a few. Since there
is often limited information to derive a parametric form, the kernel has to be
learnt in a nonparametric fashion. More importantly, the goal is a consistent
estimator that converges as data mesh refines and is robust to noise in data.
Without proper regularization, the estimator often oscillates largely from
data to data due to overfitting. Thus, regularization is crucial for the
discovery of the best kernel.

We present DARTR, a data adaptive RKHS Tikhonov regularization (DARTR) method,
for the linear inverse problem of learning of kernels in operators from data.
That is, the operator $R_{\phi}(u)$, which can be either linear or nonlinear
in $u$, depends linearly on the kernel $\phi$. We learn the kernel by
nonparametric regression that minimizes a loss functional of the mean square
error. With DARTR, our nonparametric regression algorithm produces an
estimator that converges as the data mesh refines and the rate of convergence
is robust to different levels of white noise in data. In numerical examples
including integral operators, nonlinear operators and nonlocal operators with
discrete noisy synthetic data, DARTR consistently leads to accurate
estimators, and the estimator converges at a consistent rate as the data mesh
refines under different levels of noises, outperforming two baseline
regularizers using $l^{2}$ and $L^{2}$ norms.

The major novelty of this method is the construction of a system (the
operator) intrinsic data adaptive (SIDA) RKHS, whose reproducing kernel is
encoded in the loss functional. DARTR takes the norm of this RKHS as the the
penalty norm of regularization, and ensures the learning to take place in the
function space of identifiability. Additionally, we introduce a novel
exploration measure quantifying the exploration of the kernel by the data, and
it allows a unified framework to treat SIDA-RKHS with either discrete or
continuous functions.

###  1.1 Related work

Relation to classical regression. When the data
$\\{(u_{k},f_{k})\\}_{k=1}^{N}$ are scalars instead of functions and the
operator $R_{\phi}(u)=\phi(u)$, we get back to the classical regression
problem (see e.g., [7, 13]). Our data adaptive RKHS reduces to the empirical
$L^{2}(\rho)$ space with $\rho$ being the distribution of data $\\{u\\}_{i}$,
and the DARTR regulation reduces to the classical $L^{2}$ Tikhonov/ridge
regularization (see e.g., [35]).

Relation to functional data analysis. The problem of learning kernel in
operators can be viewed as a functional data analysis (see e.g., [17, 19]),
where the data are samples from distributions on function spaces. Our ur DARTR
method is applicable to this setting. However, this study focuses on the
situation of limited deterministic data (with only a few pairs of data) and on
discovering an intrinsic low-dimensional kernel function.

Tikhonov regularization methods. DARTR differs from other Tikhonov/ridge
regularization methods at the penalty term. The commonly used penalty terms
include the Euclidean norm in the classical Tikhonov regularization ([14, 16,
11]), the RKHS norm with an ad hoc reproducing kernel ([8, 3]), the total
variation norm in the Rudin-Osher-Fatemi method in [32], or the $L^{1}$ norm
in LASSO (see e.g., [34]). Whereas each of these penalty terms has their
specific reasoning and applications, none of them take into account of the
function space of identifiability, which is fundamental for the learning of
kernels in operators.

Data-dependent function spaces. Data-dependent strategies have been explored
in the context of classical nonparametric regression, such as data-dependent
hypothesis space with an $l^{1}$ regularizer in [36, 33] and data-dependent
early stopping rule in [31]. While all strategies achieve data-dependent
regularization, only our DARTR takes into account the function space of
identifiably, which is fundamental for the learning of kernels in operators.

##  2 The inverse problem and the need of regularization

###  2.1 Problem statement: learning function parameters in operators

We consider the linear inverse problem of identifying function parameters in
operators from data. That is, given data

$\mathcal{D}=\\{(u_{k},f_{k})\\}_{k=1}^{N},\quad(u_{k},f_{k})\in\mathbb{X}\times\mathbb{Y}$
(2.1)

where $\mathbb{X}$ and $\mathbb{Y}$ are Hilbert spaces, our goal is to find a
function parameter $\phi$ in an operator $R_{\phi}:\mathbb{X}\to\mathbb{Y}$ so
that $R_{\phi}$ best fits the data pairs $\\{(u_{k},f_{k})\\}_{k=1}^{N}$:

$R_{\phi}[u]=f,$ (2.2)

where operator $R_{\phi}$ can be either linear or nonlinear in $u$ but it
depends linearly on $\phi$. In this study, we focus on such operators in the
form

$R_{\phi}[u](x)=\int_{\Omega}\phi(|y|)g[u](x,y)dy,\,\forall x\in\Omega,$ (2.3)

where $\Omega\subset\mathbb{R}^{d}$ is a bounded connected open set, $\phi$ is
a radial kernel function, $\mathbb{X}=H^{1}_{0}(\Omega)$ and
$\mathbb{Y}=L^{2}(\Omega)$. The functional $g$, which may depend on the
derivatives of $u$, is assumed to known and it specifies the form of the
operator. Examples are as follows: (see more details in Section 5)

  * •

$R_{\phi}$ is an _integral operator_ with $g[u](x,y)=u(x+y)$ and $\phi$ is
called an integral kernel.

  * •

$R_{\phi}$ is a _nonlinear operator_ with $g[u](x,y)=u^{\prime}(x+y)u(x)$ and
$\phi$ is called an interaction kernel in the mean-field equation of
interacting particles.

  * •

$R_{\phi}$ is a _nonlocal operator_ with $g[u](x,y)=u(x+y)-u(x)$ with $\phi$
called a nonlocal kernel.

These inverse problems share three common features: First, the pointwise
values of the function $\phi$ are undetermined from data, because the data
depends on $\phi$ non-locally. Also, the support of $\phi$ is unknown and is
to be learnt from data. Second, the data are discrete and can be noisy. Thus,
the inverse problem has to overcome the numerical error in the approximation
of integrals, as well as the measurement noise. Third, the inverse problem can
be extended to a homogenization problem where the operator aims to fit the
data that are not generated from the equation (2.2). In this case, the inverse
problem has to overcome the model error to identify a best fit.

###  2.2 Nonparametric regression and regularization

Our goal is to infer the kernel function $\phi$ from data in a nonparametric
fashion, so as to address the general situations that there is limited
information to derive a parametric form for the kernel. Thus, we will not
assume any constraint on the function $\phi$. More importantly, we aim for an
estimator that is consistent and resolution independent, i.e., converges in a
proper function space to the true kernel as data mesh refines and is robust to
treat noisy data.

We construct a variational estimator that minimizes loss functional (the mean
square error),

$\widehat{\phi}=\underset{\phi\in\mathcal{H}}{\operatorname{arg}\operatorname{min}}\;\mathcal{E}(\phi),\quad\text{
where
}\mathcal{E}(\phi)=\frac{1}{N}\sum_{k=1}^{N}\|R_{\phi}[u_{k}]-f_{k}\|_{\mathbb{Y}}^{2},$
(2.4)

where the hypothesis space $\mathcal{H}$ is to be selected adaptive to data.
Note that the loss functional $\mathcal{E}(\phi)$ is quadratic in $\phi$ since
the operator $R_{\phi}$ depends linearly on $\phi$. Thus, the minimizer of the
loss functional is a least square estimator. Suppose the hypothesis space is
$\mathcal{H}_{n}=\mathrm{span}\\{\phi_{i}\\}_{i=1}^{n}$ with basis functions
$\\{\phi_{i}\\}$. Then for each
$\phi=\sum_{i=1}^{n}c_{i}\phi_{i}\in\mathcal{H}_{n}$, noticing that
$R_{\phi}=\sum_{i=1}^{n}c_{i}R_{\phi_{i}}$, we can write the loss functional
in (2.4) as

$\mathcal{E}(c)=\mathcal{E}(\phi)=c^{\top}{\overline{A}}_{n}c-2c^{\top}{\overline{b}}_{n}+C_{N}^{f},$
(2.5)

where $C_{N}^{f}=\frac{1}{N}\sum_{k=1}^{N}\int_{\Omega}|f_{k}(x)|^{2}dx$ and
the normal matrix ${\overline{A}}_{n}$ and vector ${\overline{b}}_{n}$ are
given by

${\overline{A}}_{n}(i,j)=\langle\langle{\phi_{i},\phi_{j}}\rangle\rangle,\quad{\overline{b}}_{n}(i)=\frac{1}{N}\sum_{k=1}^{N}\langle{R_{\phi_{i}}[u_{k}],f_{k}}\rangle_{\mathbb{Y}},$
(2.6)

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

$\displaystyle\langle\langle{\phi,\psi}\rangle\rangle=$
$\displaystyle\frac{1}{N}\sum_{k=1}^{N}\langle{R_{\phi}[u_{k}],R_{\psi}[u_{k}]}\rangle_{\mathbb{Y}}.$
(2.7)

The least square estimator is minimizes the quadratic loss function
$\mathcal{E}(c)$:

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

where ${\overline{A}}_{n}^{-1}$ is the inverse of ${\overline{A}}_{n}$ or
Moore–Penrose pseudo-inverse when ${\overline{A}}_{n}$ is singular.

A major challenge is to find an optimal estimator capable of avoiding either
under-fitting or over-fitting, being robust to imperfect data and model error
and in particular, converging in synthetic tests when the data mesh refines.
Unfortunately, this is an ill-posed inverse problem (see Section 3.2) and the
normal matrix ${\overline{A}}_{n}$ is often highly ill-conditioned or
singular. As a result, the estimator in (2.8) oscillates largely and fails to
converge when the data mesh refines.

Various regularization methods have been introduced to prevent over-fitting in
such ill-posed inverse problems. The idea is to add a penalty term to the loss
functional:

$\mathcal{E}_{\lambda}(\phi)=\mathcal{E}(\phi)+\lambda\mathcal{R}(\phi)$ (2.9)

where $\mathcal{R}(\phi)$ is a regularization term and $\lambda$ is a
parameter which controls the importance of the regularization. Various penalty
terms have been proposed, including, for example, the Euclidean norm
$\mathcal{R}(\phi)=\|c\|^{2}$ for $\phi=\sum_{i=1}^{n}c_{i}\phi_{i}$ in the
classical Tikhonov regularization (see e.g., [35, 15]), the RKHS norm
$\mathcal{R}(\phi)=\|\phi\|_{H}^{2}$ with $H$ being a reproducing kernel
Hilbert space with an artificial reproducing kernel (see e.g., [8, 3]) , the
total variation norm $\mathcal{R}(\phi)=\|\phi^{\prime}\|_{L^{1}}$ in
Rudin–Osher–Fatemi method or the $L^{1}$ norm
$\mathcal{R}(\phi)=\|\phi\|_{L^{1}}$ in LASSO (see e.g., [34]).

Whereas each of these penalty terms has their specific reasoning and
applications, none of them take into account of the function space of
identifiability (by the loss functional, see Section 3.2), only in which the
inverse problem is well-defined. Our DARTR method (see Section 4.1) will
utilize a norm that restricts the learning in the function space of
identifiability, thus providing the most suitable regularization.

##  3 Identifiability theory and regularization

The foundation of learning is the function space of identifiability. We show
that the function space of identifiability derived by the loss functional in
(2.4) is a system intrinsic data adaptive (SIDA) RKHS. This space is the image
of the square root of the Fréchet derivative of the loss functional, which is
a compact operator. Thus the inverse problem is ill-posed since it requires
the inversion of a compact operator.

The main theme our the identifiability theory is to find the function space on
which the quadratic loss functional has a unique minimizer. In other words, we
seek the function space in which the Fréchet derivative of the loss functional
is invertible. Using the bilinear form
$\langle\langle{\cdot,\cdot}\rangle\rangle$ in (2.7), we can rewrite the loss
functional in (2.4) as

$\displaystyle\mathcal{E}(\phi)=\langle\langle{\phi,\phi}\rangle\rangle-2\frac{1}{N}\sum_{k=1}^{N}\langle{R_{\phi}[u_{k}],f_{k}}\rangle_{\mathbb{Y}}+C_{f},$
(3.1)

where $C_{N}^{f}=\frac{1}{N}\sum_{k=1}^{N}\int|f_{k}(x)|^{2}dx$. However,
there is no function space for $\phi$ yet. To start with, we introduce two key
elements: a data-adaptive exploration measure that leads to a default function
space of learning and an integral operator which plays a crucial role in our
DARTR. Throughout this section, we assume continuous data to simplify the
notation. All the integrals will be numerically approximated from discrete
data in the next section.

###### Assumption 3.1

The data $\mathcal{D}=\\{u_{k},f_{k}\\}_{k=1}^{N}$ in (2.1) are continuous
with compact support.

###  3.1 An integral operator and the SIDA-RKHS

#### The exploration measure.

We introduce first a probability measure that quantifies the exploration of
the variable of $\phi$ by the data. Given data in (2.1), we define an
empirical measure

$\displaystyle\rho(dr)=\frac{1}{ZN}\sum_{k=1}^{N}\int_{\Omega}\int_{\Omega}\delta(|y|-r)\left|g[u_{k}](x,y)\right|dxdy,$
(3.2)

where
$Z=\int_{0}^{\infty}\frac{1}{N}\sum_{k=1}^{N}\int_{\Omega}\int_{\Omega}\delta(|y|-r)\left|g[u_{k}](x,y)\right|dxdydr$
is the normalizing constant. By definition, this measure reflects the weight
being put by the loss function on $|y|$ through the data
$\\{g[u_{k}](x,y)\\}_{k=1}^{N}$.

The exploration measure plays an important role in the learning of the
function $\phi$. Its support is the region inside of which the learning
process ought to work and outside of which we have limit information from the
data to learn the function $\phi$. Thus, it defines a default function space
of learning: $L^{2}(\rho)$.

#### An integral operator.

The loss functional’s Fréchet derivative in $L^{2}(\rho)$ comes directly from
the bilinear form $\langle\langle{\cdot,\cdot}\rangle\rangle$ in (2.7). To see
this, we rewrite the bilinear form as

$\displaystyle\langle\langle{\phi,\psi}\rangle\rangle=$
$\displaystyle\frac{1}{N}\sum_{k=1}^{N}\int\left[\int\int\phi(|z|)\psi(|y|)g[u_{k}](x,z)g[u_{k}](x,y)dydz\right]dx$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\int_{0}^{\infty}\phi(r)\psi(s)G(r,s)drds=\int_{0}^{\infty}\int_{0}^{\infty}\phi(r)\psi(s){\overline{G}}(r,s)\rho(dr)\rho(ds),$
(3.3)

where the second-to-last equation follows from a change of order of
integration and a change of variables to polar coordinates with the integral
kernel $G$ given by

$G(r,s)=\frac{1}{N}\sum_{k=1}^{N}\int_{|\eta|=1}\int_{|\xi|=1}\left[\int
g[u_{k}](x,r\xi)g[u_{k}](x,s\eta)dx\right]d\xi d\eta,$ (3.4)

for $r,s\in\mathrm{supp}(\rho)$ and $G(r,s)=0$ otherwise. The last equality in
(3.3) is a re-weighting by $\rho$ with

${\overline{G}}(r,s)=\frac{G(r,s)}{\rho(r)\rho(s)},$ (3.5)

where, by an abuse of notation, we also use $\rho(r)$ to denote the density of
the probability measure $\rho$ defined in (3.2).

The next lemma shows that ${\overline{G}}$ defines a positive semi-definite
integral operator. Its proof, as well as proofs to later lemmas and theorems,
are presented in Appendix LABEL:sec:proofs.

######  Lemma 3.2 (The integral operator)

Under Assumption 3.1, the integral kernel ${\overline{G}}$ is positive semi-
definite and the integral operator
${\mathcal{L}_{\overline{G}}}:L^{2}(\rho)\to L^{2}(\rho)$

${\mathcal{L}_{\overline{G}}}\phi(r)=\int_{0}^{\infty}\phi(s)\overline{G}(r,s)\rho(s)ds$
(3.6)

is compact and positive semi-definite. Further more, for any $\phi,\psi\in
L^{2}(\rho)$,

$\langle\langle{\phi,\psi}\rangle\rangle=\langle{{\mathcal{L}_{\overline{G}}}\phi,\psi}\rangle_{L^{2}(\rho)};$
(3.7)

Proof of Lemma 3.2. Recall that a bi-variate function ${\overline{G}}$ is
positive semi-definite if 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_{i=1}^{m}\sum_{j=1}^{m}c_{i}c_{j}{\overline{G}}(r_{i},r_{j})\geq 0$.
(see e.g. [5, 8, 23]). Using (3.4) and (3.5), we have

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

Thus ${\overline{G}}$ is positive semi-definite. The operator
${\mathcal{L}_{\overline{G}}}$ is compact because ${\overline{G}}\in
L^{2}(\rho\times\rho)$, which follows from the fact that each $u_{k}$ is
bounded and the definition of $\rho$ (see also in [21]). Also, since
${\overline{G}}$ is positive semi-definite, so is
${\mathcal{L}_{\overline{G}}}$. The equation (3.7) follows from (3.3).

The next lemma provides an operator characterization of the RKHS with
${\overline{G}}$ as the reproducing kernel [2]. This RKHS is system(the
operator $R_{\phi}$) intrinsic data adaptive (SIDA), and we refer it as SIDA-
RKHS. It is the data adaptive RKHS in our DARTR.

######  Lemma 3.3 (The SIDA-RKHS)

Assume Assumption 3.1. Then the following statements hold.

  * (a)

The RKHS $H_{G}$ with ${\overline{G}}$ as the reproducing kernel satisfies
$H_{G}={\mathcal{L}_{\overline{G}}}^{1/2}(L^{2}(\rho))$ 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)}$
for any $\phi,\psi\in H_{G}$.

  * (b)

The eigen-functions of ${\mathcal{L}_{\overline{G}}}$, denoted by
$\\{\psi_{i},\psi^{0}_{j}\\}_{i,j}$ with $\\{\psi_{i}\\}$ corresponding to
positive eigenvalues $\\{\lambda_{i}\\}$ in decreasing order and
$\\{\psi^{0}_{j}\\}$ corresponding to zero eigenvalues (if any), form an
orthonormal basis of $L^{2}(\rho)$ and $\lambda_{i}$ converges to $0$.
Furthermore, for any $\phi=\sum_{i}c_{i}\psi_{i}$, we have

$\langle\langle{\phi,\phi}\rangle\rangle=\sum_{i}\lambda_{i}c_{i}^{2},\quad\|\phi\|^{2}_{L^{2}(\rho)}=\sum_{i}c_{i}^{2},\quad\|\phi\|^{2}_{H_{G}}=\sum_{i}\lambda_{i}^{-1}c_{i}^{2},$
(3.8)

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

  * (c)

For any $\phi\in L^{2}(\rho)$ and $\psi\in H_{G}$, we have

$\langle{\phi,\psi}\rangle_{L^{2}(\rho)}=\langle{{\mathcal{L}_{\overline{G}}}\phi,\psi}\rangle_{H_{G}},\quad\langle\langle{\phi,\psi}\rangle\rangle=\langle{{\mathcal{L}_{\overline{G}}}^{2}\phi,\psi}\rangle_{H_{G}}.$
(3.9)

Proof of Lemma 3.3. Part (a) is a standard operator characterization of the
RKHS $H_{G}$ (see e.g., [8, Section 4.4]).

For Part (b), since the operator ${\mathcal{L}_{\overline{G}}}$ is symmetric
positive semi-definite and compact as shown in Lemma 3.2, the eigenfunctions
are orthonormal and the eigenvalues decay to zero. The first equation in (3.8)
follows from (3.7) and the second equation follows from the orthonormality of
the eigenfunctions. At last, if $\phi\in H_{G}$, by the characterization of
the inner product of $H_{G}$ in Part (a), we have the third equation in (3.8).

The first equality in Part (c) follows from Part (a) and that
${\mathcal{L}_{\overline{G}}}^{-1/2}$ is self-adjoint, which implies that
$\langle{{\mathcal{L}_{\overline{G}}}\phi,\psi}\rangle_{H_{G}}=\langle{{\mathcal{L}_{\overline{G}}}^{1/2}\phi,{\mathcal{L}_{\overline{G}}}^{-1/2}\psi}\rangle_{L^{2}(\rho)}=\langle{\phi,\psi}\rangle_{L^{2}(\rho)}$.
The second equality in (3.9) follows from the first equality and (3.7).

###### Remark 3.4

The space $L^{2}(\rho)$ can be a discrete vector space with the function
$\phi$ defined only on finitely many points $\\{r_{i}\\}_{i=1}^{n}$ that are
explored by the data. In this setting, the integral kernel $G$ in (3.4)
becomes a positive semi-definite matrix in $\mathbb{R}^{n}$, so does
${\overline{G}}$ in (3.5). Now the integral operator
${\mathcal{L}_{\overline{G}}}$ is defined by the matrix ${\overline{G}}$ on
the weighted vector space $\mathbb{R}^{n}$ and its eigenvalues are the
generalized eigenvalues of $(G,B)$ with
$B=\mathrm{Diag}(\rho(r_{1}),\cdots,\rho(r_{n}))$. As a result, the SIDA-RKHS
$H_{G}$ is the vector space spanned by the eigenvectors with nonzero
eigenvalues. Furthermore, the norms in (3.8) can be computed directly from the
eigen-decomposition. These discrete values can be viewed piecewise constant
approximations to the functions, and the numerical algorithm in Section 4.1
applies. When $n\to\infty$, they converge to the corresponding functions under
suitable regularity conditions. Thus, the measure $\rho$ allows for a unified
framework to treat the SIDA-RKHS with either discrete or continuous functions.

###  3.2 Function space of identifiability and regularizations

We show that the function space of identifiability, i.e., on which the loss
functional has a unique minimizer (see Definition 3.5), is the subspace of
$L^{2}(\rho)$ spanned by the eigenfunctions of ${\mathcal{L}_{\overline{G}}}$
with positive eigenvalues. When zero is an eigenvalue of
${\mathcal{L}_{\overline{G}}}$, this function space is a proper subspace of
$L^{2}(\rho)$ and the loss functional has multiple minimizers in
$L^{2}(\rho)$. Thus, the inverse problem is well-defined only on this function
space. Furthermore, we show that the regularization by the SIDA-RKHS norm
enforces the regularized minimizer to be in it.

###### Definition 3.5

The function space of identifiability is the largest subspace of $L^{2}(\rho)$
in which true function $\phi_{true}$ is the unique minimizer of the loss
functional $\mathcal{E}$ with continuous noiseless data.

The next theorem characterizes the function space of identifiability.
Furthermore, it shows that this inverse problem is ill-posed since the
estimator requires the inverse of a compact operator.

######  Theorem 3.6 (Function space of identifiabilty)

Suppose that Assumption 3.1 holds. Let $\phi_{N}^{f}\in L^{2}(\rho)$ be the
Riesz representation of the bounded linear functional:

$\langle{\phi_{N}^{f},\psi}\rangle_{L^{2}(\rho)}=\frac{1}{N}\sum_{k=1}^{N}\int
2R_{\psi}[u_{k}](x)f_{k}(x)dx,\,\forall\psi\in L^{2}(\rho).$ (3.10)

Then the following statements hold.

  * (a)

The Fréchet derivative of $\mathcal{E}(\phi)$ in $L^{2}(\rho)$ is
$\nabla\mathcal{E}(\phi)=2({\mathcal{L}_{\overline{G}}}\phi-\phi_{N}^{f})$.

  * (b)

The function space of identifiability is
$H=\overline{\mathrm{span}\\{\psi_{i}\\}}$ with closure in $L^{2}(\rho)$,
where $\\{\psi_{i}\\}$ are eigenfunctions of ${\mathcal{L}_{\overline{G}}}$
with positive eigenvalues. Furthermore, the minimizer of $\mathcal{E}(\phi)$
in $H$ is $\widehat{\phi}={\mathcal{L}_{\overline{G}}}^{-1}\phi_{N}^{f}$ if
$\phi_{N}^{f}\in{\mathcal{L}_{\overline{G}}}(L^{2}(\rho))$. In particular, if
the data is perfect and generated from a true function $\phi_{true}$, we have
$\phi_{N}^{f}={\mathcal{L}_{\overline{G}}}\phi_{true}$ and
$\widehat{\phi}={\mathcal{L}_{\overline{G}}}^{-1}\phi_{N}^{f}=\phi_{true}$.

  * (c)

The Fréchet derivative of $\mathcal{E}$ in $H_{G}$ is
$\nabla^{H_{G}}\mathcal{E}(\phi)=2({\mathcal{L}_{\overline{G}}}^{2}\phi-{\mathcal{L}_{\overline{G}}}\phi_{N}^{f})$.
Its zero leads to another estimator
$\widehat{\phi}={\mathcal{L}_{\overline{G}}}^{-2}{\mathcal{L}_{\overline{G}}}\phi_{N}^{f}$
if $\phi_{N}^{f}\in{\mathcal{L}_{\overline{G}}}(L^{2}(\rho))$.

Proof of Theorem 3.6. From (3.7), we can write the loss functional in (3.1) as

$\mathcal{E}(\phi)=\langle{{\mathcal{L}_{\overline{G}}}\phi,\phi}\rangle_{L^{2}(\rho)}-2\langle{\phi_{N}^{f},\phi}\rangle_{L^{2}(\rho)}+C_{N}^{f}.$

Then we can compute the Fréchet derivative directly from definition, and Part
(a) follows.

For Part (b), first note that for any
$\phi_{N}^{f}\in{\mathcal{L}_{\overline{G}}}(L^{2}(\rho))$, the estimator
$\widehat{\phi}={\mathcal{L}_{\overline{G}}}^{-1}\phi_{N}^{f}$ is the unique
zero of the loss functional’s Fréchet derivative in $H$, hence it is the
unique minimizer of $\mathcal{E}(\phi)$ in $H$. In particular, when the
perfect data is generated from $\phi_{true}$, i.e.
$R_{\phi_{true}}[u_{k}]=f_{k}$, by (3.7) and the definition of the bilinear
form $\langle\langle{\cdot,\cdot}\rangle\rangle$ in (2.7), we have

$\langle{\phi_{N}^{f},\psi}\rangle_{L^{2}(\rho)}=\langle{{\mathcal{L}_{\overline{G}}}\phi_{true},\psi}\rangle_{L^{2}(\rho)}$

for any $\psi\in L^{2}(\rho)$. Thus,
$\phi_{N}^{f}={\mathcal{L}_{\overline{G}}}\phi_{true}$ and
$\widehat{\phi}={\mathcal{L}_{\overline{G}}}^{-1}\phi_{N}^{f}=\phi_{true}$.
That is, $\phi_{true}\in H$ is the unique minimizer of the loss functional
$\mathcal{E}$ for perfect data. Meanwhile, note that $H$ is the orthogonal
complement of the null space of ${\mathcal{L}_{\overline{G}}}$, and
$\mathcal{E}(\phi_{true}+\phi^{0})=\mathcal{E}(\phi_{true})$ for any
$\phi^{0}$ such that ${\mathcal{L}_{\overline{G}}}\phi^{0}=0$. Thus, $H$ is
the largest such function space, and we conclude that $H$ is the function
space of identifiability.

To prove Part (c), we further re-write the loss functional as

$\displaystyle\mathcal{E}(\phi)=$
$\displaystyle\langle{{\mathcal{L}_{\overline{G}}}\phi,{\mathcal{L}_{\overline{G}}}\phi}\rangle_{H_{G}}-2\langle{{\mathcal{L}_{\overline{G}}}^{1/2}\phi_{N}^{f},{\mathcal{L}_{\overline{G}}}^{1/2}\phi}\rangle_{H_{G}}+C_{N}^{f},$

which follows from (3.9) and the definition of
$\langle{\cdot,\cdot}\rangle_{H_{G}}$. Thus, by definition, the Fréchet
derivative of $\mathcal{E}(\phi)$ in the direction of $\psi\in H_{G}$ is

$\displaystyle\langle{\nabla^{H_{G}}\mathcal{E}(\phi),\psi}\rangle_{H_{G}}$
$\displaystyle=\lim_{\epsilon\to
0}\frac{1}{\epsilon}[\mathcal{E}(\phi+\epsilon\psi)-\mathcal{E}(\phi)]$
$\displaystyle=2\langle{{\mathcal{L}_{\overline{G}}}\phi,{\mathcal{L}_{\overline{G}}}\psi}\rangle_{H_{G}}-2\langle{{\mathcal{L}_{\overline{G}}}^{1/2}\phi_{N}^{f},{\mathcal{L}_{\overline{G}}}^{1/2}\psi}\rangle_{H_{G}}$
$\displaystyle=2\langle{{\mathcal{L}_{\overline{G}}}^{2}\phi-{\mathcal{L}_{\overline{G}}}\phi_{N}^{f},\psi}\rangle_{H_{G}},$

which gives the Fréchet derivative $\nabla^{H_{G}}\mathcal{E}(\phi)$.

######  Remark 3.7 (Regularization with the $L^{2}$ and the SIDA-RKHS norms)

In practice, due to the discrete and/or noisy data, we often have
$\phi_{N}^{f}={\mathcal{L}_{\overline{G}}}\phi_{true}+\phi_{1}^{\delta}+\phi_{2}^{\delta}$,
where the perturbation from the true function is decomposed to
$\phi_{1}^{\delta}\in{\mathcal{L}_{\overline{G}}}(L^{2}(\rho))$ and
$\phi^{\delta}_{2}\in{\mathcal{L}_{\overline{G}}}(L^{2}(\rho))^{\perp}$.
Clearly, when $\phi^{\delta}_{2}\neq 0$, the estimator
$\widehat{\phi}={\mathcal{L}_{\overline{G}}}^{-1}\phi_{N}^{f}$ does not exist
and regularization is necessary. Next, we compare the $L^{2}$ and the SIDA-
RKHS regularizers, i.e., consider the regularized loss functional with
$\mathcal{R}(\phi)$ being $\lambda\|\phi\|_{L^{2}}^{2}$ and
$\lambda\|\phi\|_{H_{G}}^{2}$. Then, their minimizers are

$\widehat{\phi}_{\lambda}^{L^{2}}=({\mathcal{L}_{\overline{G}}}+\lambda
I)^{-1}\phi_{N}^{f},\quad\widehat{\phi}_{\lambda}^{H_{G}}=({\mathcal{L}_{\overline{G}}}^{2}+\lambda
I)^{-1}{\mathcal{L}_{\overline{G}}}\phi_{N}^{f}.$

Plugging in
$\phi_{N}^{f}={\mathcal{L}_{\overline{G}}}\phi_{true}+\phi_{1}^{\delta}+\phi_{2}^{\delta}$,
we have

$\displaystyle\widehat{\phi}_{\lambda}^{L^{2}}$
$\displaystyle=\phi_{true}+({\mathcal{L}_{\overline{G}}}+\lambda
I)^{-1}(\phi_{1}^{\delta}-\lambda\phi_{true}+\phi_{2}^{\delta}),$
$\displaystyle\widehat{\phi}_{\lambda}^{H_{G}}$
$\displaystyle=\phi_{true}+({\mathcal{L}_{\overline{G}}}^{2}+\lambda
I)^{-1}({\mathcal{L}_{\overline{G}}}\phi_{1}^{\delta}-\lambda\phi_{true}).$

A regularizer then selects the optimal $\lambda$ to balance the errors,

$\displaystyle\|\widehat{\phi}_{\lambda}^{L^{2}}-\phi_{true}\|^{2}_{L^{2}(\rho)}$
$\displaystyle\leq\|({\mathcal{L}_{\overline{G}}}+\lambda
I)^{-1}(\phi_{1}^{\delta}+\phi_{2}^{\delta})\|^{2}+\|({\mathcal{L}_{\overline{G}}}+\lambda
I)^{-1}\lambda\phi_{true}\|^{2},$
$\displaystyle\|\widehat{\phi}_{\lambda}^{H_{G}}-\phi_{true}\|^{2}_{L^{2}(\rho)}$
$\displaystyle\leq\|({\mathcal{L}_{\overline{G}}}^{2}+\lambda
I)^{-1}{\mathcal{L}_{\overline{G}}}\phi_{1}^{\delta}\|^{2}+\|({\mathcal{L}_{\overline{G}}}^{2}+\lambda
I)^{-1}\lambda\phi_{true}\|^{2},$

where in each of them, the first term on the right hand side requires a large
$\lambda$, whereas the second term requires a small $\lambda$. In practice,
the errors $\phi_{i}^{\delta}$ are much smaller than $\phi_{true}$, and the
optimal $\lambda$ should be small so that the second term is negligible. In
this case, the bias in $\widehat{\phi}_{\lambda}^{L^{2}}$ is about
${\mathcal{L}_{\overline{G}}}^{-1}(\phi_{1}^{\delta})+\lambda^{-1}\phi_{2}^{\delta}$,
whereas the bias in $\widehat{\phi}_{\lambda}^{H_{G}}$ is about
${\mathcal{L}_{\overline{G}}}^{-1}(\phi_{1}^{\delta})$. Thus, the SIDA-RKHS
regularized estimator $\widehat{\phi}_{\lambda}^{H_{G}}$ is more accurate than
the $L^{2}$ regularized estimator. To avoid amplifying the error
$\phi_{2}^{\delta}$, a projection is necessary for the $L^{2}$ regularizer,
and we will compare the projected $L^{2}$ regularizer with the SIDA-RKHS
regularizer in Section 4.2.

##  4 Learning algorithm

###  4.1 Algorithm: nonparametric regression with DARTR

Based on the identifiability theory in Section 3.2, we introduce next a
nonparametric learning algorithm with Data Adaptive RKHS Tikhonov
Regularization (DARTR). We briefly sketch the algorithm in the following four
steps, whose details are presented in Appendix A.1.

  1. 1.

Estimate the exploration measure $\rho$. We utilize the data to estimate the
support of the true kernel and the exploration measure $\rho$. The support of
the true kernel lies in $[0,R_{0}]$ with $R_{0}$ being the diameter of the
domain $\Omega$, and it is further confined from a comparison between the
supports of $f_{k}$ and $g[u_{k}]$ (see Appendix A.1 for more details). Then,
we constrain the discrete approximation of $\rho$ defined (3.2) on the support
of $\phi$. In this process, we also assemble the regression data that will be
repeatedly used.

  2. 2.

Assemble the regression matrices and vectors. We select a class of hypothesis
spaces $\mathcal{H}_{n}=\mathrm{span}\\{\phi_{i}\\}_{i=1}^{n}$ with basis
functions $\\{\phi_{i}\\}$ and with dimension $n$ in a proper range. Then, we
compute the regression normal matrices and vectors, as well as the basis
matrix,

$\displaystyle{\overline{A}}_{n}(i,j)$
$\displaystyle=\langle\langle{\phi_{i},\phi_{j}}\rangle\rangle,\quad{\overline{b}}_{n}(i)=\langle{\phi_{N}^{f},\phi_{i}}\rangle_{L^{2}(\rho)},\quad
B_{n}(i,j)=\langle{\phi_{i},\phi_{j}}\rangle_{L^{2}(\rho)}.$ (4.1)

from data for each of these hypothesis spaces.

  3. 3.

For each triplet $({\overline{A}}_{n},{\overline{b}}_{n},B_{n})$, find the
best regularized estimator $\widehat{c}_{\lambda_{n}}$ by DARTR in Algorithm
1, as well as corresponding loss value
$\mathcal{E}(\widehat{c}_{\lambda_{n}})$.

  4. 4.

From the estimators $\\{\widehat{c}_{\lambda_{n}}\\}_{n}$, we select the one
with the smallest loss value $\mathcal{E}(\widehat{c}_{\lambda_{n}})$.

1:The regression triplet $({\overline{A}},{\overline{b}},B)$ consisting of
normal matrix ${\overline{A}}$, vector ${\overline{b}}$ and basis matrix $B$
as in (4.1).

2:SIDA-RKHS regularized estimator $\widehat{c}_{\lambda_{0}}$ and loss value
$\mathcal{E}(\widehat{c}_{\lambda_{0}})$.

3:Solve the generalized eigenvalue problem ${\overline{A}}V=BV\Lambda$, where
$\Lambda$ is the diagonal matrix of eigenvalues and the matrix $V$ has columns
being eigenvectors orthonormal in the sense that $V^{\top}BV=I$.

4:Compute the RKHS-norm matrix $B_{rkhs}=(V\Lambda V^{\top})^{-1}$, using
pseudo inverse when $\Lambda$ is singular. We refer to Remark A.1 on a
computational technique to avoid the inverse matrix.

5:Use the L-curve method to find an optimal estimator
$\widehat{c}_{\lambda_{0}}$: select $\lambda_{0}$ maximizing the curvature of
the $\lambda$-curve
$(\log\mathcal{E}(\widehat{c}_{\lambda}),\log(\widehat{c}_{\lambda}^{\top}B_{rkhs}\widehat{c}_{\lambda}))$,
where the least squares estimator
$\widehat{c}_{\lambda}=({\overline{A}}+\lambda B_{rkhs})^{-1}{\overline{b}}$
minimizes the regularized loss function
$\mathcal{E}_{\lambda}(c)=\mathcal{E}(c)+\lambda c^{\top}B_{rkhs}c\,\,\text{
with
}\,\mathcal{E}(c)=c^{\top}{\overline{A}}c-2c^{\top}{\overline{b}}+{\overline{b}}^{\top}{\overline{A}}^{-1}{\overline{b}},$
where the matrix inversion is a pseudo-inverse when it is singular.

Algorithm 1 Data Adaptive RKHS Regularization (DARTR).

In comparison to the classical nonparametric regression using
$({\overline{A}}_{n},{\overline{b}}_{n})$, we need only an additional basis
matrix $B_{n}$. The novelty of our algorithm is the data adaptive components,
such as the exploration measure $\rho$, the basis matrix $B_{n}$ in
$L^{2}(\rho)$ and the norm of the SIDA-RKHS for regularization. The
computation of the SIDA-RKHS norm is based on the generalized eigenvalues
problem with the pair $({\overline{A}}_{n},B_{n})$, whose eigenvalues
approximate the eigenvalues of ${\mathcal{L}_{\overline{G}}}$ in (3.6) and
$\widehat{\psi}_{k}=V_{jk}\phi_{j}$ approximates the eigenfunctions of
${\mathcal{L}_{\overline{G}}}$ (see Theorem 4.1). The additional computational
cost is only the generalized eigenvalue problem which can be solved
efficiently.

###### Theorem 4.1

Let $\mathcal{H}_{n}=\mathrm{span}\\{\phi_{i}\\}_{i=1}^{n}\subset L^{2}(\rho)$
and let $({\overline{A}}_{n},B_{n})$ be the normal and basis matrix in (4.1).
Assume that $\mathcal{H}_{n}$ is large enough so that
${\mathcal{L}_{\overline{G}}}(L^{2}(\rho))\subset\mathcal{H}_{n}$ (which is
true, for example when $\rho$ is a discrete-measure on a discrete set
$\mathcal{R}$ and $\\{\phi_{n}\\}$ are piecewise constant functions with
$n=|\mathcal{R}|$). Then, the operator ${\mathcal{L}_{\overline{G}}}$ in (3.6)
has eigenvalues $(\lambda_{1},\ldots,\lambda_{n})$ solved by the generalize
eigenvalue problem

${\overline{A}}_{n}V=B_{n}\Lambda V,\quad
s.t.,V^{\top}B_{n}V=I_{n},\quad\Lambda=\mathrm{Diag}(\lambda_{1},\ldots,\lambda_{n}),$
(4.2)

and the corresponding eigenfunctions of ${\mathcal{L}_{\overline{G}}}$ are
$\\{\psi_{k}=\sum_{j=1}^{n}V_{jk}\phi_{j}\\}$.

Proof of Theorem 4.1. Let $\psi_{k}=\sum_{j=1}^{n}V_{jk}\phi_{j}$ with
$V^{\top}B_{n}V=I_{n}$. Then, $\psi_{k}$ is an eigenfunction of
${\mathcal{L}_{\overline{G}}}$ with eigenvalue $\lambda_{k}$ if and only if
for each $i$,

$\displaystyle\langle{\phi_{i},\lambda_{k}\psi_{k}}\rangle_{L^{2}(\rho)}=\langle{\phi_{i},{\mathcal{L}_{\overline{G}}}\psi_{k}}\rangle_{L^{2}(\rho)}=\sum_{j=1}^{n}\langle{\phi_{i},{\mathcal{L}_{\overline{G}}}\phi_{j}}\rangle_{L^{2}(\rho)}V_{jk}=\sum_{j=1}^{n}{\overline{A}}_{n}(i,j)V_{jk},\vspace{-2mm}$

where the last equality follows from the definition of ${\overline{A}}_{n}$ in
(4.1). Meanwhile, by the definition of $B_{n}$ we have
$\langle{\phi_{i},\lambda_{k}\psi_{k}}\rangle_{L^{2}(\rho)}=\sum_{j=1}^{n}B_{n}(i,j)\lambda_{k}V_{jk}$
for each $i$. Then, Equation (4.2) follows.

###  4.2 Comparison with projected $l^{2}$ and $L^{2}$ regularizers

Our DARTR method differs from other regularizers in its use of the SIDA-RKHS
norm, which restricts the function to be in the function space of
identifiability. In the following, we compare it with the $l^{2}$ and $L^{2}$
regularizers that apply regularization terms
$\mathcal{R}(\phi)=\|\phi\|_{l^{2}}^{2}=\sum_{i}c_{i}^{2}$ or
$\mathcal{R}(\phi)=\|\phi\|_{L^{2}}^{2}=c^{\top}B_{n}c$. In fact, a direct
application of these two regularization terms would lead to problematic
regularizers with largely biased estimators when ${\overline{A}}_{n}$ is
singular, i.e., when the function space of identifiability is a proper
subspace of $L^{2}(\rho)$, because the inverse problem is ill-defined on
$L^{2}(\rho)$. Thus, in practice, one makes a projection to the function space
of identifiability (i.e., the image of ${\overline{A}}_{n}$ in computation)
before adding these regularization terms, and we call them projected $l^{2}$
and $L^{2}$ regularizers.

Table 1: The SIDA-RKHS regularizer v.s. the projected $l^{2}$,$L^{2}$
regularizers∗.

| $l^{2}$ | $L^{2}$ | SIDA-RKHS
---|---|---|---
$\mathcal{R}(\phi)$ | $\|c\|^{2}=c^{\top}c$ | $\|c\|_{B_{n}}^{2}=c^{\top}B_{n}c$ | $\|c\|_{H_{G}}^{2}=c^{\top}B_{rkhs}c$
$c_{\lambda}$ | $c_{\lambda}=\sum_{i=1}^{k}\frac{1}{\sigma_{i}+\lambda}u_{i}^{\top}{\overline{b}}$ | $c_{\lambda}=\sum_{i=1}^{k}\frac{1}{\lambda_{i}+\lambda}v_{i}^{\top}{\overline{b}}$ | $c_{\lambda}=\sum_{i=1}^{k}\frac{1}{\lambda_{i}+\lambda\lambda_{i}^{-1}}v_{i}^{\top}{\overline{b}}$
SVD |  ${\overline{A}}_{n}=\sum_{i=1}^{n}\sigma_{i}u_{i}u_{i}^{\top}$, $u_{i}^{\top}u_{j}=\delta_{ij}$ |  ${\overline{A}}_{n}=\sum_{i=1}^{n}\lambda_{i}v_{i}v_{i}^{\top}$, $v_{i}^{\top}B_{n}v_{j}=\delta_{ij}$
|  $U^{\top}{\overline{A}}_{n}U=\Sigma$, $U^{\top}U=I$ |  $V^{\top}{\overline{A}}_{n}V=\Lambda$, $V^{\top}B_{n}V=I$

*All regularizers estimate $\phi=\sum_{i=1}^{n}c_{i}\phi_{i}$ from ${\overline{A}}_{n}c={\overline{b}}_{n}$ with basis matrix $B_{n}$ (see (4.1)). The projected $l^{2}$ and $L^{2}$ regularizers use only the non-zero eigenvalues $\\{\sigma_{i}\\}_{i=1}^{k}$ and $\\{\lambda_{i}\\}_{i=1}^{k}$ and their eigenvectors.

Table 1 compares our SIDA-RKHS regularizer with the projected $l^{2}$ and
$L^{2}$ regularizers. We note that there are the following connections:

  * •

The $L^{2}$ regularizer is a basis-adaptive generalization of the $l^{2}$
regularizer. When $B_{n}=I$ (i.e., the basis $\\{\phi_{i}\\}$ are orthonormal
in $L^{2}(\rho)$), the two are the same. When $B_{n}$ is not the identity
matrix (i.e., the basis $\\{\phi_{i}\\}$ are not orthonormal in
$L^{2}(\rho)$), which happens often, the $L^{2}$ regularizer takes it into
account through the generalized eigenvalue problem.

  * •

The SIDA-RKHS regularizer is an improvement over the $L^{2}$ regularizer. When
all the generalized eigenvalues are $\lambda_{i}\equiv 1$ (e.g.,
${\mathcal{L}_{\overline{G}}}$ is an identity operator or when
${\overline{A}}_{n}=B_{n}$ as in classical regression), the two are the same.
Otherwise, the SIDA-RKHS regularizer not only takes into account of the
function space of identifiability but also a balance between $\lambda_{i}$ and
$\lambda_{i}^{-1}$.

  * •

The SIDA-RKHS regularizer restricts the learning to be in the function space
of identifiability by definition, while the other two regularizers, if not
projected, miss this fundamental issue.

##  5 Numerical results

We test our learning method on three types of operators: linear integral
operators, nonlocal operators and nonlinear operators. For each type, we
systematically examine the method in the regimes of noiseless and noisy data,
with kernels in and out of the SIDA-RKHSs. Since the ground-truth kernel is
known, we study the convergence of estimators to the true kernel as the data
mesh refines. Thus, the regularization has to overcome both numerical error
and noise in the imperfect data. The codes implementing the proposed algorithm
are shared on GitHub111URL: https://github.com/LearnDynamics/DARTR_kernel.

###  5.1 Settings and main results

The settings of the numerical tests for all three types of operators are as
below.

Comparison with baseline methods. On each dataset, we compare our SIDA-RKHS
regularizer with two baseline regularizers using the projected $l^{2}$ and
$L^{2}$ norm (denoted by l2 and L2 in the figures below, respectively) defined
in Section 4.2. All three regularizers use the same L-curve method to select
the hyper-parameter $\lambda$ as described in Appendix A.1. They differ only
at the regularization norm.

Settings of synthetic data. We test two kernels for each type of operators:
$\bullet$ _Truncated sine kernel._ The truncated sine kernel is
$\phi_{true}=\sin(2x)\mathbf{1}_{[0,3]}(x)$. It represents a kernel with
discontinuity. Due to the nonlocal dependency of the operator on the kernel,
this discontinuity can cause a global bias to the estimator.
$\bullet$ _Gaussian kernel._ The kernel $\phi_{true}$ is the Gaussian density
centered at 3 with standard deviation 0.75. It represents a smooth kernel
whose interaction concentrated in the middle of its support.

The kernels act on the same set of function $\\{u_{k}\\}_{k=1,2}$ with
$u_{1}=\sin(x)\mathbf{1}_{[-\pi,\pi]}(x)$ and
$u_{2}(x)=\sin(2x)\mathbf{1}_{[-\pi,\pi]}(x)$. When generating the data for
learning, the integral $R_{\phi}[u_{k}]=f_{k}$ is computed by the adaptive
Gauss-Kronrod quadrature method. This integrator 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 of
$\\{u_{k},f_{k}\\}_{k=1}^{N}=\\{u_{k}(x_{j}),f_{k}(x_{j}):x_{j}\in{[-40,40]},j=1,\ldots,J\\}_{k=1}^{N}$,
where $x_{j}$ is a point on the uniform grid with mesh size $\Delta x$. For
the nonlinear operator, to avoid the inverse problem being ill-defined, we
introduce add an additional pair of data $(u_{3},f_{3})$ with
$u_{3}(x)=x\mathbf{1}_{[-\pi,\pi]}(x)$ (see Section 5.3 for more details). In
short, the discrete data $\\{u_{k}\\}_{k=1,2}$ are continuous functions and
the discrete data $u_{3}$ is a piece-wise smooth function.

For each kernel, we consider both noiseless and noisy data with different
noise levels by taking values of noise-to-signal-ratio ($\mathrm{nsr}$) in
$\\{0,0.5,1,1.5,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_{k}$.

Settings for the learning algorithm. When estimating the kernels from the
discrete data, we estimate the values of the kernel on the points
$\mathcal{S}=\\{r_{j}\\}_{j=1}^{J}$ with $r_{j}=j\Delta x$, the support of the
empirical exploration measure $\rho$. When the data mesh refines, the size of
this set increases. In terms of the algorithm in Section 4.1, such a discrete
estimation uses a hypothesis space with B-spline basis functions consisting of
piece-wise constants with knots being the points in $\mathcal{S}$. Thus, this
hypothesis space has the largest dimension for the basis matrix $B_{n}$ in
(4.1) being non-singular, and there is no need to select an optimal dimension.
In this setting, the regularizer is the only source of regularization and
there is no regularization from basis functions. Hence, this setting
highlights the role of the regularizers.

Performance assessment. We assess the performance of the regularizers by their
ability to consistently identify the true kernels in the presence of numerical
error (in the Riemann sum approximation of the integrals due to discrete data)
and noise (due to noisy data). We present typical estimators, the
$L^{2}(\rho)$ errors of the estimators as data mesh refines, as well as the
statistics (mean and standard deviation) of the rates of convergence that are
computed from 20 independent simulations.

#### Summary of main results

Our main finding are as follows.

  * •

The SIDA-RKHS regularizer’s estimators are the most accurate as long as the
optimal regularization parameter $\lambda$ is properly selected. However, the
selection of $\lambda$ depends on multiple factors, ranging from the form of
the operator, the numerical approximation, the noise and the treatment of the
singular or ill-conditioned normal matrix. As a result, it occasionally
happens that the $l^{2}$ or $L^{2}$ regularizer performs slightly better.
Thus, in addition to accuracy of the estimator, it is important to also
compare the consistency of convergence rates for different levels of noise.

  * •

The SIDA-RKHS regularizer robustly leads to estimators converging at a
consistent rate for all levels of noises for each operator, while the other
two regularizers cannot.

  * •

The rate of convergence of the SIDA-RKHS regularizer’s estimator from noisy
data depends on both the continuity of the kernel and the continuity of the
discrete data: when the two matches, the rate is higher and closer to 1, as
shown in Table 2.

Table 2: Rate of convergence of the SIDA-RKHS regularizer’s estimators from
noisy data.

| Linear Integral Operator | Nonlinear Operator | Nonlocal Operator
---|---|---|---
Kernel | Data continuity(C) | Data continuity (D) | Data continuity (C)
Truncated Sine (D) | 0.29 | 0.94 | 0.29
Gaussian (C) | 0.62 | 0.66 | 1.01

* Here “C” stands for continuous, and “D” stands for discontinuous. When the continuity of the kernel and data matches, the rates are closer to $1$ than when the two dis-matches. The rates are the average of the mean rates for $\mathrm{nsr}\in\\{0.1,0.5,1,2\\}$ in the right columns of Figure 1-3. We do not report the rate for the $l^{2}$ and $L^{2}$ regularizers because they do not have a consistent rate.

###  5.2 Linear integral operators

We consider first the integral operator with kernel $\phi$:

$R_{\phi}[u](x)=\int_{\Omega}\phi(|y-x|)u(y)dy=f(x).$

After a change of variables in the integral, it is the operator $R_{\phi}$ in
(2.3) with $g[u](x,y)=u(x+y)$. Such kernels in operators arise in a wide range
of applications, such as the Green’s function in PDEs (see e.g., [10, 12]) and
convolution kernels in image processing in [29], to name just a few.

For this operator, the exploration measure $\rho$ (defined in (3.2)) is a
uniform measure, since each data $g[u_{k}]$ interacts with the kernel
uniformly. Furthermore, since each $g[u_{k}]$ is continuous, the reproducing
kernel ${\overline{G}}$ in (3.4) is continuous on the support of $\rho$, thus
the SIDA-RKHS consists of continuous functions. As a result, we expect the
algorithm to learn the smooth Gaussian kernel better than the discontinuous
truncated sine kernel.

Figure 1: Linear integral operators with the sine kernel (top row) and
Gaussian kernel (bottom row). Left Column: typical estimators by the three
regularizers, in comparison of the true kernel, superimposed with the
exploration measure $\rho$ (in cyan color), when $\Delta x=0.05$ and noise-to-
signal-ratio $\mathrm{nsr}=1$. Middle 3-Columns: convergence of estimators as
the data mesh-size $\Delta x$ refines, along with values of the loss function.
Right Column: the mean and standard deviation of the convergence rates in 20
independent simulations, with five levels of noise (with
$\mathrm{nsr}\in\\{0,0.1,0.5,1,2\\}$). Only the SIDA-RKHS regularizer’s
estimator consistently converges for all levels of noise, and its estimators
are mostly more accurate than those of the other two regularizers.

Figure 1 shows the results. The left column shows the typical estimators by
the three regularizers, in comparison of the true kernel, when $\Delta x=0.05$
and noise-to-signal-ratio $\mathrm{nsr}=1$. The exploration measure $\rho$ (in
light cyan color) is uniform for each kernel, and its support, estimated from
the difference between the supports of $g[u_{k}]$ and $f_{k}$, is slightly
larger than the support of the true kernel. All three regularizers lead to
accurate estimators. The RKHS regularizer’s estimators are the closest to the
true kernel and this is further verified in the middle 3-column panel with
$\Delta x=0.05$ add $\mathrm{nsr}=1$: for the sine kernel, all three
estimators’ $L^{2}(\rho)$ errors are about $10^{-1}$; but for the Gaussian
kernel, the RKHS’s estimator has an error close to $10^{-2.5}$ while the other
two regularizers’ error are about $10^{-2}$.

The middle 3-column panel shows the convergence of the estimator’s
$L^{2}(\rho)$ error as the data mesh refines when $\mathrm{nsr}=0.1$ and
$\mathrm{nsr}=1$, superimposed with the corresponding values of the loss
function. When $\mathrm{nsr}=1$, all three regularizers’ estimators converge
for both kernels, at rates that are close to the rates of the loss function,
and their errors are comparable. However, when $\mathrm{nsr}=0.1$, the RKHS
regularizer continues to yield converging estimators, whereas the other two
regularizers have flat error lines even though the corresponding loss values
keep decaying. In particular, those flat error lines are above those errors
for $\mathrm{nsr}=1$ with $\Delta x\leq 0.025$, i.e., when the numerical error
is small. Thus, these results demonstrates the importance to take into account
the function space of learning via SIDA-RKHS, particularly when the noise
level is relatively low.

The right column shows the mean and standard deviations of the rates of
convergence in 20 independent simulations. The RKHS regularizer has consistent
rates of convergence for all levels of noises. The rates are closer to 1 for
the smooth Gaussian kernel (which matches the continuity of data) than the
rates for the discontinuous truncated sine kernel when the data are noisy. The
rates are close to 1 when the data are noiseless. On the other hand, the
$l^{2}$ and $L^{2}$ regularizers fails to have consistent rates when the noise
level reduces. In particular, for the sine kernel, they present deceivingly
higher rates than the RKHS regularizer when $\mathrm{nsr}\in\\{0.5,1,2\\}$,
and the middle 3-column panel reveals the facts: they often have much larger
errors than the RKHS when $\Delta x=0.2$, thus leading to deceiving better
rates even when their errors remains large as $\Delta x$ decreases.

In short, the RKHS regularizer leads to estimators that converge consistently,
at lower rates for the discontinuous sine kernel (which is discontinuous,
different from the dta) and at higher rates for the smooth Gaussian kernel
(which match the continuity of the data), while the $l^{2}$ and $L^{2}$
regularizers cannot. Furthermore, RKHS regularizer’s estimators are often more
accurate than those of the other two regularizers.

###  5.3 Nonlinear operators

Next we consider the nonlinear operator $R_{\phi}$ with
$g[u](x,y)=\partial_{x}[u(x+y)u(x)]$:

$R_{\phi}[u](x)=\int_{\Omega}\phi(|y|)\partial_{x}[u(x+y)u(x)]dy=[u*\phi(|\cdot|)u]^{\prime}(x).$

Such nonlinear operators arise in the mean-field equations of interaction
particles (see e.g.,[18, 28, 26, 22]), and the function $\phi$ is called an
interaction kernel. More precisely, the mean-field equations are of the form
$\partial_{t}u=\nu\Delta u+\mathrm{div}(u*K_{\phi}u)$ on $\mathbb{R}^{d}$,
where $K_{\phi}(y)=\phi(|y|)\frac{y}{|y|}$. Here we consider only $d=1$ and
neglect the ratio $\frac{y}{|y|}$ to obtain the above operator.

We add an additional pair of data $(u_{3},f_{3})$ with
$u_{3}(x)=x\mathbf{1}_{[-\pi,\pi]}(x)$, so as to avoid the issue that the
value of $[u*\phi(|\cdot|)u](x)$ is under-determined from the data
$f(x)=[u*\phi(|\cdot|)u]^{\prime}(x)$ due to the differential. Here we set the
derivative of $u_{3}$ to be $u_{3}^{\prime}(x)=\mathbf{1}_{[-\pi,\pi]}(x)$.
These derivatives are approximated by finite difference when learning the
kernel from discrete data. Note that the $u_{3}$ and its derivative have jump
discontinuities. As a result, the reproducing kernel ${\overline{G}}$ in (3.4)
also has discontinuity, and the SIDA-RKHS contains discontinuous functions.

Figure 2: Nonlinear operators, in the same setting as in Figure 1. The SIDA-
RKHS regularizer’s estimators are significantly more accurate than those of
the $l^{2}$ and $L^{2}$ regularizers in the left column. The middle 3-column
panel shows that the SIDA-RKHS regularizer leads to consistently converging
estimators as the data mesh refines, for both levels of noise, while the other
two regularizers have slower and less consistent error decay and their error
lines flatten when the noise level is $\mathrm{nsr}=1$. The right column shows
that only the SIDA-RKHS regularizer has consistent rates for all levels of
noise, and the other two regularizers’ rates drops significantly when the
noise level increases.

Figure 2 shows the results. The left column shows that the exploration measure
$\rho$ is non-uniform due to the nonlinear function $g[u_{k}]$, and its
density is a decreasing function, suggesting that the data explores the short
range interactions more than the long range interaction. The RKHS
regularizer’s estimators significantly outperforms the other two regularizers,
and they are near smooth and are close to the true kernels. The $l^{2}$ and
$L^{2}$ regularizers have largely oscillating estimators, suggesting an
overfitting. Note that the RKHS estimators also have oscillating parts, but
they are only in the region where the exploration measure has little weight,
due to limited data exploration. The superior performance of RKHS regularizer
is further verified in the middle 3-column panel with $\Delta x=0.05$ add
$\mathrm{nsr}=1$: its errors are much smaller than those of the other two
regularizers.

The middle 3-column panel shows that the RKHS regularizer’s error consistently
decreases as the data mesh refines. In contrast, the other two regularizers
have slower and less consistent error decay, in particular, their error lines
flatten as the noise level increases.

The right column shows that the RKHS regularizer has consistent rates of
convergence for all levels of noises, with all rates close to 1 for the sine
kernel, and slightly above 0.5 for the Gaussian kernel. In comparison, the
other two regularizers’ rates decreases as the noise level increases, dropping
to close zero when the noise level is $\mathrm{nsr}=2$.

In short, the RKHS regularizer’s estimators are more accurate than those of
the $l^{2}$ and $L^{2}$ regularizers. More importantly, the RKHS regularizer
consistently leads to convergent estimators, maintaining similar rates for all
levels of noises, at rates close to $1$ for the truncated sine kernel (which
is discontinuous, matching the discontinuity of data) and at rates slightly
above $0.5$ for the Gaussian kernel (which is smooth, different from the
data). The $l^{2}$ and $L^{2}$ regularizers have convergent estimators, but
the rate of convergence drops when the noise level increases.

###  5.4 Nonlocal operators

At last, we consider nonlocal operators $R_{\phi}$ with
$g[u](x,y)=u(x+y)-u(x)$:

$R_{\phi}[u](x)=\int_{\Omega}\phi(|y|)[u(x+y)-u(x)]dy.$

Such nonlocal operators arise in various areas such as nonlocal and fractional
diffusions (see e.g.,[9, 1, 6]) and such nonlocal operators have been used to
construct homogenized models for peridynamic in [38, 37].

Figure 3: Nonlocal operators, in the same setting as in Figure 1. The left
column show that all regularizers lead to accurate estimators. The middle
3-column panel shows that the SIDA-RKHS regularizer leads to converging
estimators as the data mesh refines for two levels of noise, though at a slow
rate for the sine kernel. The $l^{2}$ and $L^{2}$ regularizers have less
consistent error decay for different noise levels and different kernels.
Overall, the SIDA-RKHS estimators have the smallest error mostly. The right
column shows that only the SIDA-RKHS regularizer has consistent rates for all
levels of noise.

Figure 3 shows the results. The left column shows typical estimators. The
exploration measure $\rho$ shrinks to zero near the origin due to the
difference $g[u]=u(y)-u(x)$ and the continuity of $u$. All three regularizers
lead to accurate estimators, and the RKHS estimator is the most accurate.

In the middle 3-column panel, we observe again that the RKHS regularizer leads
to estimators remain converging as data mesh refines for both noise levels,
even though the errors decay slower than the loss function. On the other hand,
the $l^{2}$ and $L^{2}$ regularizers have inconsistent error decay: the errors
decreasing monotonically when $\mathrm{nsr}=1$, but the error lines oscillate
when $\mathrm{nsr}=0.1$ for the sine kernel, and for the Gaussian kernel, they
present deceiving rates larger than the decay of the loss function due their
large errors when $\Delta x$ is large.

The right column further confirms the consistency of the RKHS regularizer’s
rates and the inconsistency of the $l^{2}$ and $L^{2}$-regularizers’ rates.
When the data is noisy, the rates of the RKHS regularizer are about 0.29 for
the truncated sine kernel (which has a jump discontinuity) and about 1 for the
Gaussian kernel (which is continuous). Meanwhile, the rates for the $l^{2}$
and $L^{2}$-regularizers are about 0.65 for the sine kernel, and about 0.8 for
the Gaussian kernel. We note again that they can have deceivingly better rates
than the RKHS regularizer’s while their errors are larger. Moreover, when the
data is noiseless, RKHS regularizer has rates close to 1 for both kernels,
while the other two regularizers rates are not consistent.

##  6 Discussion and future work

We have proposed a data adaptive RKHS Tikhonov regularization (DARTR) method
for the nonparametric learning of kernel functions in operators. The DARTR
method regularizes the least squares regression by the norm of a system
intrinsic and data adaptive (SIDA) RKHS, which constraints the learning to the
function space of identifiability.

Our numerical tests on synthetic data suggests that DARTR has the following
advantages: (1) it is naturally adaptive to both data and the operator; (2) it
is robust to numerical error due to discrete data and white noise in data,
leading to estimators converging at a consistent rate for different levels of
noises; (3) it is computationally as efficient as classical nonparametric
least squares regression methods, requiring in addition only an exploration
measure and a basis matrix that come with negligible computational cost.

This study presents a preliminary introduction of the DARTR method. There are
several directions for further development and analysis of DARTR in general
settings and applications:

  1. 1.

Convergence analysis. We are in short of a convergence analysis of the
regularized estimators due to the numerical errors in the normal matrix.

  2. 2.

Multivariate kernel functions. When the kernel is a multivariate function,
sparse-grid representation or sparse basis functions (sparse polynomials)
become necessary. A related issue is to select the optimal dimension of the
hypothesis space.

  3. 3.

Applications to Bayesian inverse problems. In a Bayesian perspective, the
Tikhonov regularization can be interpreted as a Gaussian prior with a
covariance matrix corresponding to the penalty term. In this perspective, our
SIDA-RKHS norm coincides with the Zellner’s g-prior ([40, 4]) that uses
${\overline{A}}_{n}^{-1}$ as prior covariance, because we have
$B_{rkhs}={\overline{A}}_{n}^{-1}$ when the basis functions are orthonormal in
$L^{2}(\rho)$.

  4. 4.

The DARTR method is applicable to general linear inverse problems with a
quadratic loss functional. It is particularly useful when the data depends on
the unknown function non-locally.

##  Appendix A Algorithm details

###  A.1 Detailed nonparametric learning algorithm

We consider only discrete data $\\{u_{k}(x_{j}),f_{k}(x_{j})\\}_{k=1}^{N}$ in
1-dimensional and at equidistant mesh points $\\{x_{j}=j\Delta
x\\}_{j=0}^{J}$. The extension to multi-dimensional cases is straightforward.

#### Step 1: Estimate the exploration measure and assemble regression data.

We first estimate the exploration measure and extract the regression data that
can be used for all hypothesis spaces by utilizing the regression structure
and reading the data only once. This step can reduce the computational cost in
orders of magnitude when the data is large with thousands of pairs
$(u_{k},f_{k})$ with fine mesh.

Let $R_{0}$ be the diameter of the set $\Omega$. The discrete data set
$\\{u_{k}(x_{j}),f_{k}(x_{j})\\}_{k=1}^{N}$ explores only the variable $r$ of
$\phi$ in the set $\mathcal{R}_{N}^{J}=\\{r_{ijk}=|y_{i}|\leq
R_{0}:g[u_{k}](x_{i},y_{j})\neq 0\text{ for some }i,j,k\\}$, the set of all
values explored by data with repetition. A discrete approximation of the
exploration measure $\rho$ in (3.2) is

$\displaystyle\rho_{N}^{J}(dr)$
$\displaystyle=\frac{1}{|\mathcal{R}_{N}^{J}|}\sum_{k=1}^{N}\sum_{i,j=1}^{J}\delta_{|y_{i}|}(r)|g[u_{k}](x_{j},y_{i})|.$
(A.1)

This measure $\rho_{N}^{J}$ uses only the information from $u_{k}$ and it does
not reflect the information about the kernel in $f_{k}$. To estimate the
support of the kernel, we extract the additional information from
$\\{f_{k}\\}$ 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}\\},$
(A.2)

where $(L^{u}_{i},R^{u}_{i})$ and $(L^{f}_{i},R^{f}_{i})$ are the lower and
upper bounds of the supports $g[u_{k}](x,y)$ and $\mathrm{supp}(f_{k})$
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_{k}$ and $f_{k}$, which extracts the additional
information in the data $\\{f_{k}\\}$. 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.

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. Accordingly, we constrain the exploration measure to be
supported in $[0,R]$, and for simplicity of notation, we still denote it by
$\rho_{N}^{J}$.

Assemble regression data.  Next, we assemble the regression data that will be
used repeatedly, thus saving the computational cost by orders of magnitude,
particularly when the data size is large with thousands of pairs
$(u_{k},f_{k})$. In order to compute the normal matrix
${\overline{A}}(i,j)=\langle\langle{\phi_{i},\phi_{j}}\rangle\rangle$ for any
pair of basis functions, with the bilinear form defined in (3.3), we only need
the integral kernel $G$. In particular, when $d=1$, the integral
$\int_{|\eta|=1}h(\eta)d\eta=h(\eta)+h(-\eta)$, therefore, we have

$\displaystyle
G(r,s)=\frac{1}{N}\sum_{k=1}^{N}\int_{\Omega}\left(g[u_{k}](x,r)+g[u_{k}](x,-r)\right)\left(g[u_{k}](x,s)+g[u_{k}](x,-s)\right)dx$
(A.3)

for $r,s\in\mathrm{supp}(\rho)$. Similarly, for a basis function $\phi_{i}$,
to compute ${\overline{b}}(i)$ in (2.6), which can be re-written as

$\displaystyle{\overline{b}}_{n}(i)=\frac{1}{N}\sum_{k=1}^{N}\int
R_{\phi_{i}}[u_{k}](x)f_{k}(x)dx=\int_{0}^{R}\phi_{i}(r)g_{N}^{f}(r)dr,$ (A.4)

we only need the function $g_{N}^{f}$ defined by

$g_{N}^{f}(r)=\frac{1}{N}\sum_{k=1}^{N}\int_{\Omega}\left(g[u_{k}](x,r)+g[u_{k}](x,-r)\right)f_{k}(x)dx.$
(A.5)

Let $r_{k}=k\Delta x$ for $k=1,\ldots,\lfloor{\frac{R}{\Delta x}}\rfloor$,
which are the mesh points of $\phi$ explored by the data. Then, all the
regression data we need in the original data (2.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\\},$ (A.6)

where $G$, $g_{N}^{f}$ and $\rho_{N}^{J}$ are defined respectively in (A.3),
(A.5) and (A.1).

#### Step 2: Select a class of hypothesis spaces and assemble regression
matrices and vectors.

We set a class of data-adaptive hypothesis spaces
$\mathcal{H}_{n}=\mathrm{span}\\{\phi_{i}\\}_{i=1}^{n}$ with their dimensions
set to range from under-fitting to over-fitting. The basis functions can be
either global basis functions such as polynomials and trigonometric functions,
or local basis functions such B-spline polynomials (see e.g., Chapter 2 of
[30] and [27]). 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}$ (A.7)

is nonsingular. 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
nonsingular 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$ with $\phi_{i}(x)=\delta(x_{i}-x)$, 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}$.

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},$ (A.8) $\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_{i}(r_{k})\phi_{j}(r_{k})\rho_{N}^{J}(r_{k})\Delta
x.$

The triplet $({\overline{A}}_{n},{\overline{b}}_{n},B_{n})$ is all we need for
regression with regularization in the next step.

#### Step 3: Regression with DARTR.

Our DARTR method uses the norm of the SIDA-RKHS, which is the function space
of identifiability as discussed in Section 3.2. That is, our estimator is the
minimizer of the regularized loss in (2.9) with the regularization norm
$\mathcal{R}(\phi)=\|\phi\|^{2}_{H_{G}}$ defined in (3.8).

Computation of the RKHS norm In practice, 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 three 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 (3.6), 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},$ (A.9)

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}.$
(A.10)

Then, we select the hyper-parameter $\lambda$ by the L-curve method (see
Section A.2).

######  Remark A.1 (Least squares to avoid matrix inverse)

The matrix inverses can cause numerical issues when the normal matrix
${\overline{A}}$ is ill-conditioned or singular. Fortunately, the matrix
inversions in $B_{rkhs}$ and in solving $({\overline{A}}_{n}+\lambda
B_{rkhs})c_{\lambda}={\overline{b}}_{n}$ can be avoided by using minimum norm
least squares solution. Note that this linear equation is equivalent to
$(B_{rkhs}^{-T/2}{\overline{A}}_{n}B_{rkhs}^{-1/2}+\lambda
I)\widetilde{c}_{\lambda}=B_{rkhs}^{-T/2}{\overline{b}}_{n}$ with
$\widetilde{c}_{\lambda}=B_{rkhs}^{-1/2}c_{\lambda}$, where $B_{rkhs}^{-T/2}$
is the transpose of the square root matrix $B_{rkhs}^{-1/2}$. Meanwhile, the
square root $B_{rkhs}^{-1/2}=(V\Lambda V^{\top})^{1/2}$ comes directly from
(A.9). Thus, these treatments avoid the matrix inversions and lead to more
robust estimators.

We summarize the method in Algorithm 2.

Algorithm 2 Nonparametric learning of the nonlocal kernel with spare-aware
regularization

1:The data
$\\{u_{k},f_{k}\\}_{k=1}^{N}=\\{u_{k}(x_{j}),f_{k}(x_{j})\\}_{k,j=1}^{N,J}$
with $x_{j}=j\Delta x$ to construct the nonlocal model $R_{\phi}[u]=f$.

2:Estimator $\widehat{\phi}$

3:Estimate the exploration measure $\rho_{N}^{J}$ from data as in (A.1), and
estimate the support of the kernel from data as in (A.2). Let $R$ be the upper
bound of the support.

4:Get regression data $(G,g_{N}^{f})$ in (A.6).

5:Select a class of hypothesis spaces
$\mathcal{H}_{n}=\mathrm{span}\\{\phi_{i}\\}_{i=1}^{n}$ by selecting a type of
basis functions, e.g., polynomials or B-splines, $n$ in the range
$\lfloor{\frac{R}{\Delta x}}\rfloor\times[0.2,1]$.

6:For each $n$, compute $({\overline{A}}_{n},{\overline{b}}_{n},B_{n})$ as in
(A.8) for $\mathcal{H}_{n}=\mathrm{span}\\{\phi_{i}\\}_{i=1}^{n}$, using
$(G,g_{N}^{f},\rho_{N}^{J})$ obtained above. If the basis matrix $B_{n}$ is
singular, remove $n$ from the range. For the
$({\overline{A}}_{n},{\overline{b}}_{n},B_{n})$, find the best regularized
estimator $\widehat{c}_{\lambda_{n}}$ by DARTR in Algorithm 1, as well as
corresponding loss value $\mathcal{E}(\widehat{c}_{\lambda_{n}})$.

7: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}$.

###  A.2 Hyper-parameter by the L-curve method

We select the parameter $\lambda$ by the L-curve method [16, 22]. 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\smash{\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}},}$
(A.11)

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 [16] for more details).

#### Acknowledgments

The work of F.L. is partially funded by NSF Award DMS-1913243. FL would like
to Professor Yue Yu and Professor Mauro Maggioni for inspiring discussions.

## References

  * [1] D. Applebaum.  Lévy processes and stochastic calculus.  Cambridge university press, 2009.
  * [2] N. Aronszajn.  Theory of reproducing kernels.  Transactions of the American mathematical society, 68(3):337–404, 1950.
  * [3] F. Bauer, S. Pereverzev, and L. Rosasco.  On regularization algorithms in learning theory.  Journal of complexity, 23(1):52–72, 2007.
  * [4] M. J. Bayarri, J. O. Berger, A. Forte, and G. García-Donato.  Criteria for bayesian model choice with application to variable selection.  The Annals of Statistics, 40(3), Jun 2012.
  * [5] C. Berg, J. P. R. Christensen, and P. Ressel.  Harmonic analysis on semigroups: theory of positive definite and related functions, volume 100.  New York: Springer, 1984.
  * [6] C. Bucur and E. Valdinoci.  Nonlocal Diffusion and Applications, volume 20 of Lecture Notes of the Unione Matematica Italiana.  Springer International Publishing, Cham, 2016.
  * [7] F. Cucker and S. Smale.  On the mathematical foundations of learning.  Bulletin of the American mathematical society, 39(1):1–49, 2002\.
  * [8] F. Cucker and S. Smale.  Emergent behavior in flocks.  IEEE Transactions on Automatic Control, 52(5):852–862, 2007.
  * [9] Q. Du, M. Gunzburger, R. B. Lehoucq, and K. Zhou.  Analysis and Approximation of Nonlocal Diffusion Problems with Volume Constraints.  SIAM Rev., 54(4):667–696, 2012.
  * [10] L. C. Evans.  Partial differential equations, volume 19.  American Mathematical Soc., 2010.
  * [11] S. Gazzola, P. C. Hansen, and J. G. Nagy.  Ir tools: a matlab package of iterative regularization methods and large-scale test problems.  Numerical Algorithms, 81(3):773–811, 2019.
  * [12] C. R. Gin, D. E. Shea, S. L. Brunton, and J. N. Kutz.  Deepgreen: deep learning of green’s functions for nonlinear boundary value problems.  Scientific reports, 11(1):1–14, 2021.
  * [13] L. Györfi, M. Kohler, A. Krzyzak, and H. Walk.  A distribution-free theory of nonparametric regression.  Springer Science & Business Media, 2006.
  * [14] P. C. Hansen.  REGULARIZATION TOOLS: A Matlab package for analysis and solution of discrete ill-posed problems.  Numer Algor, 6(1):1–35, 1994.
  * [15] P. C. Hansen.  Rank-deficient and discrete ill-posed problems: numerical aspects of linear inversion.  SIAM, 1998.
  * [16] P. C. Hansen.  The L-curve and its use in the numerical treatment of inverse problems.  pages 119–142, 2000.
  * [17] T. Hsing and R. Eubank.  Theoretical foundations of functional data analysis, with an introduction to linear operators, volume 997.  John Wiley & Sons, 2015.
  * [18] P.-E. Jabin and Z. Wang.  Mean field limit for stochastic particle systems.  In Active Particles, Volume 1, pages 379–402. Springer, 2017.
  * [19] H. Kadri, E. Duflos, P. Preux, S. Canu, A. Rakotomamonjy, and J. Audiffren.  Operator-valued Kernels for Learning from Functional Response Data.  page 54.
  * [20] N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, and A. Anandkumar.  Neural operator: Learning maps between function spaces.  arXiv preprint arXiv:2108.08481, 2021.
  * [21] Q. Lang and F. Lu.  Identifiability of interaction kernels in mean-field equations of interacting particles.  arXiv preprint arXiv:2106.05565, 2021.
  * [22] Q. Lang and F. 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.
  * [23] Z. Li, F. Lu, M. Maggioni, S. Tang, and C. Zhang.  On the identifiability of interaction functions in systems of interacting particles.  Stochastic Processes and their Applications, 132:135–163, 2021\.
  * [24] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar.  Neural operator: Graph kernel network for partial differential equations.  arXiv preprint arXiv:2003.03485, 2020.
  * [25] C. Lin, Z. Li, L. Lu, S. Cai, M. Maxey, and G. E. Karniadakis.  Operator learning for predicting multiscale bubble growth dynamics.  The Journal of Chemical Physics, 154(10):104118, 2021.
  * [26] F. Lu, M. Maggioni, and S. Tang.  Learning interaction kernels in stochastic systems of interacting particles from multiple trajectories.  Foundations of Computational Mathematics, pages 1–55, 2021.
  * [27] T. Lyche, C. Manni, and H. Speleers.  Foundations of Spline Theory: B-Splines, Spline Approximation, and Hierarchical Refinement, volume 2219, pages 1–76.  Springer International Publishing, Cham, 2018.
  * [28] S. Motsch and E. Tadmor.  Heterophilious Dynamics Enhances Consensus.  SIAM Rev, 56(4):577 – 621, 2014.
  * [29] H. Owhadi and G. R. Yoo.  Kernel flows: From learning kernels from data into the abyss.  Journal of Computational Physics, 389:22–47, 2019.
  * [30] L. Piegl and W. Tiller.  The NURBS Book.  Monographs in Visual Communication. Springer Berlin Heidelberg, Berlin, Heidelberg, 1997.
  * [31] G. Raskutti, M. J. Wainwright, and B. Yu.  Early stopping and non-parametric regression: an optimal data-dependent stopping rule.  The Journal of Machine Learning Research, 15(1):335–366, 2014.
  * [32] L. I. Rudin, S. Osher, and E. Fatemi.  Nonlinear total variation based noise removal algorithms.  Physica D: nonlinear phenomena, 60(1-4):259–268, 1992.
  * [33] L. Shi, Y.-L. Feng, and D.-X. Zhou.  Concentration estimates for learning with $l^{1}$-regularizer and data dependent hypothesis spaces.  Applied and Computational Harmonic Analysis, 31(2):286–302, 2011\.
  * [34] R. Tibshirani.  Regression Shrinkage and Selection Via the Lasso.  Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267–288, 1996.
  * [35] A. N. Tihonov.  Solution of incorrectly formulated problems and the regularization method.  Soviet Math., 4:1035–1038, 1963.
  * [36] H. Wang.  Analysis of statistical learning algorithms in data dependent function spaces.  PhD thesis, City University of Hong Kong, 2009.
  * [37] H. You, Y. Yu, S. Silling, and M. D’Elia.  Data-driven learning of nonlocal models: From high-fidelity simulations to constitutive laws.  ArXiv201204157 Cs Math, 2020.
  * [38] H. You, Y. Yu, S. Silling, and M. D’Elia.  A data-driven peridynamic continuum model for upscaling molecular dynamics.  Computer Methods in Applied Mechanics and Engineering, 389:114400, 2022.
  * [39] H. You, Y. Yu, N. Trask, M. Gulian, and M. D’Elia.  Data-driven learning of nonlocal physics from high-fidelity synthetic data.  Computer Methods in Applied Mechanics and Engineering, 374:113553, 2021.
  * [40] A. Zellner and A. Siow.  Posterior odds ratios for selected regression hypotheses.  Trabajos de Estadistica Y de Investigacion Operativa, 31:585–603, 1980.