###### Abstract

Surrogate modeling based on Gaussian processes (GPs) has received increasing
attention in the analysis of complex problems in science and engineering.
Despite extensive studies on GP modeling, the developments for functional
inputs are scarce. Motivated by an inverse scattering problem in which
functional inputs representing the support and material properties of the
scatterer are involved in the partial differential equations, a new class of
kernel functions for functional inputs is introduced for GPs. Based on the
proposed GP models, the asymptotic convergence properties of the resulting
mean squared prediction errors are derived and the finite sample performance
is demonstrated by numerical examples. In the application to inverse
scattering, a surrogate model is constructed with functional inputs, which is
crucial to recover the reflective index of an inhomogeneous isotropic
scattering region of interest for a given far-field pattern.

Functional-Input Gaussian Processes with Applications to Inverse Scattering
Problems
Chih-Li Sung∗,a,†, Wenjia Wang***These authors contributed equally to the
manuscript.,b, Fioralba Cakonic,†, Isaac Harrisd,†, Ying Hungc†††The authors
gratefully acknowledge funding from NSF DMS-2113407, DMS-2107891, and
CCF-1934924.
aMichigan State University
bHong Kong University of Science and Technology
cRutgers, the State University of New Jersey
dPurdue University

Keywords: Computer Experiments; Emulation; Surrogate Model; Uncertainty
Quantification

##  1 Introduction

Computer experiments, the studies of real systems using mathematical models
such as partial differential equations, have received increasing attention in
science and engineering for the analysis of complex problems. Typically,
computer experiments require a great deal of time and computing. Therefore,
based on a finite sample of computer experiments, it is crucial to build a
surrogate for the actual mathematical models and use the surrogate for
prediction, inference, and optimization. The Gaussian process (GP) model, also
called kriging, is a widely used surrogate model due to its flexibility,
interpolating property, and the capability of uncertainty quantification
through the predictive distribution. More discussions on computer experiments
and surrogate modeling using GP models can be found in Santner et al., (2018)
and Gramacy, (2020).

This paper is motivated by an inverse scattering problem in computer
experiments, where the computer experiments involve functional inputs and
therefore the analysis and inference rely on a surrogate model that can take
into account functional inputs. Figure 1 illustrates the idea of inverse
scattering. Let the functional input $g$ represent the material properties of
an inhomogeneous isotropic scattering region of interest shown in the middle
of Figure 1. For a given functional input, the far-field pattern, $u^{s}$, is
obtained by solving partial differential equations (Cakoni et al.,, 2016)
which is computationally intensive. Given a new far-field pattern, the goal of
inverse scattering is to recover the functional input using a surrogate model.
Therefore, a crucial step to address this problem is to develop a surrogate
model applicable to functional inputs. Problems with functional inputs are
frequently found in engineering applications of non-destructive testing where
measurements on the surface or exterior of an object is used to infer the
interior structure. Similar problems also come up in electrical impedance
tomography where one wishes to recover the functional input representing the
electric conductivity, from the measured current to voltage mapping; see,
e.g., Mueller and Siltanen, (2020), for the electrical impedance tomography
model and Cakoni et al., (2016) for the model studied here for inverse
scattering.

Figure 1: Illustration of the inverse scattering problem.

Despite extensive studies on surrogate modeling using GPs, the developments
for functional inputs are scarce. Analysis of functional input is often
addressed by the techniques developed in functional data analysis literature
including functional linear models (Malfait and Ramsay,, 2003; Hall and
Horowitz,, 2007), generalized functional linear models (James,, 2002; Wang et
al.,, 2010; Chen et al.,, 2011), and varying-coefficient models (Fan and
Zhang,, 1999). There are recent developments on GP models. For example, Nguyen
and Peraire, (2015) propose a functional-input GP with bilinear covariance
operators and apply it to linear partial differential equations. Morris,
(2012) develops a kriging model with a covariance function specifically for
time-series data. Tan, (2019) introduces a functional-input GP model for
finite element models. Chen et al., (2021) propose a spectral-distance
correlation function and apply it to 3D printing. To the best of our
knowledge, most of the existing methods are restrictive to specific
applications, while a systematical study is still lacking.

To address this problem, a general framework of GP modeling with a new class
of kernel functions is introduced for the construction of surrogate models
with functional inputs. The remainder of the paper is organized as follows. In
Section 2, a functional-input GP model is introduced. A new class of kernel
functions, including a linear and a nonlinear kernel, and their theoretical
properties are discussed in Section 3. Numerical analysis is conducted in
Section 4 to examine the prediction accuracy of the proposed models. In
Section 5, the proposed framework is applied to construct a surrogate model
for an inverse scattering problem. Concluding remarks are given in Section 6.
Detailed theoretical proofs are provided in Supplementary Materials.

##  2 Functional-Input Gaussian Process

Suppose that $V$ is a functional space consisting of functions defined on a
compact and convex region $\Omega\subseteq\mathbb{R}^{d}$, and all functions
$g\in V$ are continuous. A functional-input GP, $f:V\rightarrow\mathbb{R}$, is
denoted by

$\displaystyle f(g)\sim\mathcal{FIGP}(\mu,K(g,g^{\prime})),$ (1)

where $\mu$ is an unknown mean and $K(g,g^{\prime})$ is a semi-positive kernel
function for $g,g^{\prime}\in V$. A new class of kernel $K(g,g^{\prime})$ for
functional inputs is discussed in Section 3.

Given a properly defined kernel function, the estimation and prediction
procedures are similar to the conventional GP. Assume that there are $n$
realizations from the functional-input GP, where $g_{1},\ldots,g_{n}$ are the
inputs and $f(g_{1}),\ldots,f(g_{n})$ are the outputs. We have
$f(g_{1}),\ldots,f(g_{n})$ following a multivariate normal distribution,
$\mathcal{N}_{n}(\boldsymbol{\mu}_{n},\mathbf{K}_{n})$, with mean
$\boldsymbol{\mu}_{n}=\mu\boldsymbol{1}_{n}$ and covariance $\mathbf{K}_{n}$,
where $\boldsymbol{1}_{n}$ is a size-$n$ all-ones vector and
$(\mathbf{K}_{n})_{j,k}=K(g_{j},g_{k})$. The unknown parameters, including
$\mu$ and the hyperparameters associated with the kernel function, can be
estimated by likelihood-based approaches or Bayesian approaches. We refer the
details of the estimation methods to Santner et al., (2018) and Gramacy,
(2020).

Suppose $g\in V$ is an untried new function. By the property of the
conditional multivariate normal distribution, the corresponding output $f(g)$
follows a normal distribution with the mean and variance,

$\displaystyle\mathbb{E}[f(g)|\mathbf{y}_{n}]=$
$\displaystyle\mu+\mathbf{k}(g)^{T}\mathbf{K}^{-1}_{n}(\mathbf{y}_{n}-\boldsymbol{\mu}_{n}),\quad\text{and}$
(2) $\displaystyle\mathbb{V}[f(g)|\mathbf{y}_{n}]=$ $\displaystyle
K(g,g)-\mathbf{k}(g)^{T}\mathbf{K}^{-1}_{n}\mathbf{k}(g),$ (3)

where $\mathbf{y}_{n}=(y_{1},...,y_{n})$, $y_{i}=f(g_{i})$, and
$\mathbf{k}(g)=(K(g,g_{1}),...,K(g,g_{n}))^{T}$. The conditional mean of (2)
is used to predict $f(g)$, and the conditional variance of (3) can be used to
quantify the prediction uncertainty.

##  3 A New Class of Kernel Functions

In this section, a new class of kernel functions for the functional-input GP
is introduced. Based on the proposed models, the asymptotic convergence
properties of the resulting mean squared prediction errors are derived.
Section 3.1 focuses on the discussions of a linear kernel and Section 3.2
extends the discussions to a nonlinear kernel. A practical guidance on the
selection of optimal kernel is discussed in Section 3.3. For notational
simplicity, the mean in (1) is assumed to be zero in this section but the
results can be easily extended to non-zero cases.

###  3.1 Linear kernel for functional inputs

We first introduce a linear kernel for functional inputs $g_{1}$ and $g_{2}$:

$\displaystyle
K(g_{1},g_{2})=\int_{\Omega}\int_{\Omega}g_{1}(\mathbf{x})g_{2}(\mathbf{x}^{\prime})\Psi(\mathbf{x},\mathbf{x}^{\prime}){\rm{d}}\mathbf{x}{\rm{d}}\mathbf{x}^{\prime},$
(4)

where $g_{1},g_{2}\in V$ and $\Psi$ is a positive definite function defined on
$\Omega\times\Omega$. It can be shown in the following proposition that this
kernel function is semi-positive definite.

######  Proposition 3.1.

The linear kernel $K$ defined in (4) is semi-positive definite on $V\times V$.

By Mercer’s theorem (Rahman,, 2007), we have

$\displaystyle\Psi(\mathbf{x},\mathbf{x}^{\prime})=\sum_{j=1}^{\infty}\lambda_{j}\phi_{j}(\mathbf{x})\phi_{j}(\mathbf{x}^{\prime}),$
(5)

where $\mathbf{x},\mathbf{x}^{\prime}\in\Omega$, and
$\lambda_{1}\geq\lambda_{2}\geq...>0$ and $\\{\phi_{k}\\}_{k\in\mathbb{N}}$
are the eigenvalues and the orthonormal basis in $L_{2}(\Omega)$,
respectively. Given the positive definite function $\Psi$, we can construct a
GP via the Karhunen–Loève expansion:

$\displaystyle
f(g)=\sum_{j=1}^{\infty}\sqrt{\lambda_{j}}\langle\phi_{j},g\rangle_{L_{2}(\Omega)}Z_{j},$
(6)

where $Z_{j}$’s are independent standard normal random variables, and
$\langle\phi_{j},g\rangle_{L_{2}(\Omega)}$ is the inner product of $\phi_{j}$
and $g$, which is
$\langle\phi_{j},g\rangle_{L_{2}(\Omega)}=\int_{\Omega}\phi_{j}(\mathbf{x})g(\mathbf{x}){\rm{d}}\mathbf{x}$.
It can be shown that the covariance function of the constructed GP in (6) is
$K(g_{1},g_{2})$ defined in (4), i.e.,

$\displaystyle{\rm Cov}(f(g_{1}),f(g_{2}))=$
$\displaystyle\sum_{j=1}^{\infty}\lambda_{j}\langle\phi_{j},g_{1}\rangle_{L_{2}(\Omega)}\langle\phi_{j},g_{2}\rangle_{L_{2}(\Omega)}=\sum_{j=1}^{\infty}\lambda_{j}\int_{\Omega}\int_{\Omega}g_{1}(\mathbf{x})\phi_{j}(\mathbf{x})g_{2}(\mathbf{x}^{\prime})\phi_{j}(\mathbf{x}^{\prime}){\rm{d}}\mathbf{x}{\rm{d}}\mathbf{x}^{\prime}$
$\displaystyle=$
$\displaystyle\int_{\Omega}\int_{\Omega}g_{1}(\mathbf{x})g_{2}(\mathbf{x}^{\prime})\Psi(\mathbf{x},\mathbf{x}^{\prime}){\rm{d}}\mathbf{x}{\rm{d}}\mathbf{x}^{\prime}$
(7)

for any $g_{1},g_{2}\in V$.

######  Proposition 3.2.

The Gaussian process, $f(g)$, constructed as in (6) is linear, i.e., for any
$a,b\in\mathbb{R}$ and $g_{1},g_{2}\in V$, it follows that
$f(ag_{1}+bg_{2})=af(g_{1})+bf(g_{2})$.

###### Proof of Proposition 3.2.

For any $a,b\in\mathbb{R}$ and $g_{1},g_{2}\in V$, it follows that

$\displaystyle f(ag_{1}+bg_{2})=$
$\displaystyle\sum_{j=1}^{\infty}\sqrt{\lambda_{j}}\langle\phi_{j},ag_{1}+bg_{2}\rangle_{L_{2}(\Omega)}Z_{j}=\sum_{j=1}^{\infty}\sqrt{\lambda_{j}}\left(a\langle\phi_{j},g_{1}\rangle_{L_{2}(\Omega)}+b\langle\phi_{j},g_{2}\rangle_{L_{2}(\Omega)}\right)Z_{j}$
$\displaystyle=$ $\displaystyle af(g_{1})+bf(g_{2}),$

which finishes the proof. ∎

The proposed kernel function has an intuitive interpretation that connects to
Bayesian modeling. In a Bayesian linear regression, the conditional mean
function is assumed to be $f(\mathbf{x})=\mathbf{x}^{T}\mathbf{w}$, where
$\mathbf{w}$ is typically assumed to have a multivariate normal prior, i.e.,
$\mathbf{w}\sim\mathcal{N}(0,\Sigma_{\mathbf{w}})$. Hence, for any two points,
$\mathbf{x}$ and $\mathbf{x}^{\prime}$, the covariance of $f(\mathbf{x})$ and
$f(\mathbf{x}^{\prime})$ is
$\text{Cov}(f(\mathbf{x}^{\prime}),f(\mathbf{x}^{\prime}))=\mathbf{x}^{T}\Sigma_{\mathbf{w}}\mathbf{x}^{\prime}$,
which can be interpreted as a weighted inner product of $\mathbf{x}$ and
$\mathbf{x}^{\prime}$. The proposed model (6) can be viewed as an analogy to
the Bayesian linear model with functional inputs and the covariance (3.1) can
be viewed as a weighted inner product of the two functions $g_{1}$ and
$g_{2}$.

To understand the performance of the proposed predictor of (2) with the kernel
function defined in (4), we first characterize the mean squared prediction
error by the following theorem. Denote the reproducing kernel Hilbert space
(RKHS) associated with a kernel $\Psi$ as $\mathcal{N}_{\Psi}(\Omega)$.

######  Theorem 3.3.

Let $\hat{f}(g)=\mathbb{E}[f(g)|\mathbf{y}_{n}]$ as in (2). For any continuous
function $g\in V\subset L_{2}(\Omega)$, define a linear operator on
$L_{2}(\Omega)$:

$\displaystyle\mathcal{T}g(\mathbf{x})=\int_{\Omega}g(\mathbf{x}^{\prime})\Psi(\mathbf{x},\mathbf{x}^{\prime}){\rm{d}}\mathbf{x}^{\prime}.$

The mean squared prediction error (MSPE) of $\hat{f}(g)$ can be written as

$\displaystyle\mathbb{E}\left(f(g)-\hat{f}(g)\right)^{2}=\min_{(u_{1},...,u_{n})\in\mathbb{R}^{n}}\left\|\mathcal{T}g-\sum_{j=1}^{n}u_{j}\mathcal{T}g_{j}\right\|_{\mathcal{N}_{\Psi}(\Omega)}^{2},$
(8)

where $\left\|\cdot\right\|_{\mathcal{N}_{\Psi}(\Omega)}$ is the RKHS norm of
$\mathcal{N}_{\Psi}(\Omega)$.

By Proposition 10.28 in Wendland, (2004), it is straightforward to shown that
$\mathcal{T}g\in\mathcal{N}_{\Psi}(\Omega)$ and therefore the right-hand side
of (8) exists. According to Theorem 3.3, the MSPE can be represented as the
distance between $\mathcal{T}g$ and its projection on the linear space spanned
by $\\{\mathcal{T}g_{1},...,\mathcal{T}g_{n}\\}$. This distance can be reduced
if $g_{j}$’s can be designed so that the spanned space can well approximate
the space $V$. We highlight some designs of $g_{j}$’s in the following two
corollaries where the convergence rates of MSPE can be explicitly discussed.

In the following corollaries, the kernel function $\Psi$ is assumed to be
Matérn (Stein,, 1999):

$\displaystyle\Psi(\mathbf{x},\mathbf{x}^{\prime})=\psi(\|\Theta(\mathbf{x}-\mathbf{x}^{\prime})\|_{2})$
(9)

with

$\displaystyle\psi(r)=\frac{\sigma^{2}}{\Gamma(\nu)2^{\nu-1}}(2\sqrt{\nu}r)^{\nu}B_{\nu}(2\sqrt{\nu}r),$
(10)

where $\Theta$ is a scale parameter, which is a $d\times d$ positive diagonal
matrix, $\|\cdot\|_{2}$ denotes the Euclidean norm, $\sigma^{2}$ is a positive
scalar, $B_{\nu}$ is the modified Bessel function of the second kind, and
$\nu$ represents a smoothness parameter. Without loss of generality, we assume
that $\Theta$ is an identity matrix and $\sigma^{2}=1$ for the theoretical
developments in this section. More detailed discussions of these parameters,
including $\Theta,\sigma^{2}$ and $\nu$, are given in Section 4.

######  Corollary 3.4.

Suppose $g_{j}$, $j=1,\ldots,n$ are the first $n$ eigenfunctions of $\Psi$,
i.e, $g_{j}=\phi_{j}$. For $g\in V\subset L_{2}(\Omega)$, there exists a
constant $C_{1}>0$ such that

$\displaystyle\mathbb{E}\left(f(g)-\hat{f}(g)\right)^{2}\leq
C_{1}\|g\|_{L_{2}(\Omega)}^{2}n^{-\frac{2\nu}{d}}.$ (11)

Furthermore, if $g\in\mathcal{N}_{\Psi}(\Omega)$, then there exists a constant
$C_{2}>0$ such that

$\displaystyle\mathbb{E}\left(f(g)-\hat{f}(g)\right)^{2}\leq
C_{2}\|g\|_{\mathcal{N}_{\Psi}(\Omega)}^{2}n^{-\frac{4\nu}{d}}.$ (12)

Corollary 3.4 shows that if we can design the input functions to be
eigenfunctions of $\Psi$, the convergence rate of MSPE is polynomial. If the
functional space is further assumed to be the RKHS associated with the kernel
$\Psi$ (i.e. $g\in\mathcal{N}_{\Psi}(\Omega)$), which is smaller than
$L_{2}(\Omega)$, the convergence rate becomes faster as indicated by (12).

Instead of selecting the input functions to be eigenfunctions, an alternative
is to design the input functions through a set of knots in $\Omega$, i.e.,
$\mathbf{X}_{n}\equiv\\{\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\\}$, where
$\mathbf{x}_{j}\in\Omega$ for $j=1,\ldots,n$, and the convergence rate is
derived in the following corollary.

######  Corollary 3.5.

  * (1)

Suppose $g_{j}(\mathbf{x})=\Psi(\mathbf{x},\mathbf{x}_{j})$, where
$\mathbf{x},\mathbf{x}_{j}\in\Omega$ for $j=1,\ldots,n$. For
$g\in\mathcal{N}_{\Psi}(\Omega)$, there exists a constant $C_{3}>0$ such that

$\displaystyle\mathbb{E}\left(f(g)-\hat{f}(g)\right)^{2}\leq
C_{3}\|g\|_{\mathcal{N}_{\Psi}(\Omega)}^{2}h_{\mathbf{X}_{n},\Omega}^{2\nu},$

where $h_{\mathbf{X}_{n},\Omega}$ is the fill distance of $\mathbf{X}_{n}$,
i.e.,

$\displaystyle
h_{\mathbf{X}_{n},\Omega}:=\sup_{\mathbf{x}\in\Omega}\inf_{\mathbf{x}_{j}\in\mathbf{X}_{n}}\|\mathbf{x}-\mathbf{x}_{j}\|_{2}.$

  * (2)

For a quasi-uniform design, there exists a positive constant $C$ such that
$h_{\mathbf{X}_{n},\Omega}\leq Cn^{-1/d}$ (Wendland,, 2004; Müller,, 2009).
Therefore, there exists a constant $C_{4}>0$ such that

$\displaystyle\mathbb{E}\left(f(g)-\hat{f}(g)\right)^{2}\leq
C_{4}\|g\|_{\mathcal{N}_{\Psi}(\Omega)}^{2}n^{-\frac{2\nu}{d}}.$ (13)

As compared with the results in Corollary 3.4, the convergence rate of quasi-
uniform designs as shown in (13) is slower than the choice of eigenfunctions
as shown in (12). Despite a slower rate of convergence, designing functional
inputs by $\Psi(\mathbf{x}_{j},\cdot)$ with space-filling $\mathbf{x}_{j}$’s
can be relatively easier in practice than finding eigenfunctions of $\Psi$.

The proposed linear kernel can be naturally modified to accommodate the
potential non-linearity in $f$ by enlarging the feature space using a pre-
sepecified nonlinear transformation $\mathcal{M}$ on $g$, i.e.,
$\mathcal{M}:V\rightarrow V_{1}$, where $V_{1}$ is a function class. The
resulting kernel function can be written as

$\displaystyle K(g_{1},g_{2})=\int_{\Omega}\int_{\Omega}\mathcal{M}\circ
g_{1}(\mathbf{x})\mathcal{M}\circ
g_{2}(\mathbf{x}^{\prime})\Psi(\mathbf{x},\mathbf{x}^{\prime}){\rm{d}}\mathbf{x}{\rm{d}}\mathbf{x}^{\prime},$

and the corresponding GP can be constructed by

$\displaystyle
f(g)=\sum_{j=1}^{\infty}\sqrt{\lambda_{j}}\langle\phi_{j},\mathcal{M}\circ
g\rangle_{L_{2}(\Omega)}Z_{j}.$ (14)

The convergence results of MSPE can be extended to (14). There are many
possible ways to define $\mathcal{M}$ so that the feature space can be
enlarged; however, the flexibility comes with a higher estimation complexity.
In the next section, we propose an alternative to address the non-linearity
through a kernel function, which is computationally more efficient.

###  3.2 Nonlinear kernel for functional inputs

In this section, we introduce a new type of kernel function for functional
inputs that takes into account the non-linearity via a radial basis function.
Let $\psi(r):\mathbb{R}^{+}\rightarrow\mathbb{R}$ be a radial basis function
whose corresponding kernel in $\mathbb{R}^{d}$ is strictly positive definite
for any $d\geq 1$. Note that the radial basis function of (10), whose
corresponding kernel is Matérn, satisfies this condition. Define $K:V\times
V\rightarrow\mathbb{R}$ as

$\displaystyle K(g_{1},g_{2})=\psi(\gamma\|g_{1}-g_{2}\|_{L_{2}(\Omega)}),$
(15)

where $\|\cdot\|_{L_{2}(\Omega)}$ is the $L_{2}$-norm of a function, defined
by $\|g\|_{L_{2}(\Omega)}=(\langle g,g\rangle_{L_{2}(\Omega)})^{1/2}$, and
$\gamma$ is a parameter that controls the decay of the kernel function with
respect to the $L_{2}$-norm.

######  Proposition 3.6.

The function $K$ defined in (15) is positive definite on $V\times V$.

Assume that there exists a probability measure $P$ on $V$ such that
$\int_{V}g(t)^{2}dP(g)<\infty$, for $t\in\Omega$ (Ritter,, 2007). Based on the
positive definite function $K$ as in (15), we can construct a GP via the
Karhunen–Loève expansion:

$\displaystyle f(g)=\sum_{j=1}^{\infty}\sqrt{\lambda_{j}}\varphi_{j}(g)Z_{j},$
(16)

where $\varphi_{j}$’s are the orthonormal basis obtained by a generalized
version of Mercer’s theorem,
$K(g_{1},g_{2})=\sum_{j=1}^{\infty}\lambda_{j}\varphi_{j}(g_{1})\varphi_{j}(g_{2})$
(Steinwart and Scovel,, 2012) with respect to the probability measure $P$.
Based on the proposed nonlinear kernel, the convergence rate of the MSPE is
studied in the following theorem.

######  Theorem 3.7.

Suppose that $\Phi$ is a Matérn kernel function with smoothness $\nu_{1}$, and
$\psi$ is the radial basis function of (10), whose corresponding kernel is
Matérn with smoothness $\nu$. Let $\tau=\min(\nu,1)$. For any $n>4$, there
exist $n$ input functions such that for any $g\in\mathcal{N}_{\Phi}(\Omega)$
with $\|g\|_{\mathcal{N}_{\Phi}(\Omega)}\leq 1$, the MSPE can be bounded by

$\displaystyle\mathbb{E}\left(f(g)-\hat{f}(g)\right)^{2}\leq C_{5}(\log
n)^{-\frac{(\nu_{1}+d/2)\tau}{d}}\log\log n.$ (17)

Based on (17), it appears that the convergence rate is slower than the
conventional GP where the inputs are defined in the Euclidean space. Although
potentially there are rooms to further sharpen this rate, a slower rate of
convergence for functional inputs is not surprising because the input space is
much larger than the Euclidean space.

###  3.3 Selection of kernels

Based on Sections 3.1 and 3.2, the linear kernel of (4) results in a less
flexible model leading to a lower prediction variance but higher bias, while
the nonlinear kernel of (15) results in a more flexible model leading to a
higher variance but lower bias (Hastie et al.,, 2009). To find an optimal
kernel function that balances the bias–variance trade-off, the idea of cross-
validation is adopted that allows us to select the kernel by minimizing the
estimated prediction error.

Although cross-validation methods are typically expensive to implement in many
situations, the leave-one-out cross-validation (LOOCV) of GP models can be
expressed in a closed form, which makes computation less demanding (Zhang and
Wang,, 2010; Rasmussen and Williams,, 2006; Currin et al.,, 1988).
Specifically, denote $\tilde{y}_{i}$ as the prediction mean based on all data
except $i$th observation and $y_{i}$ as the real output of $i$th observation.
One can show that, for a kernel candidate $K$, which can be either the linear
kernel (4) or the nonlinear kernel (15), the LOOCV error is

$\frac{1}{n}\sum^{n}_{i=1}(y_{i}-\tilde{y}_{i})^{2}=\frac{1}{n}\|\boldsymbol{\Lambda}_{n}^{-1}\mathbf{K}_{n}^{-1}(\mathbf{y}_{n}-\mu\mathbf{1}_{n})\|^{2}_{2},$
(18)

where $\boldsymbol{\Lambda}_{n}$ is a diagonal matrix with the element
$(\boldsymbol{\Lambda}_{n})_{j,j}=(\mathbf{K}_{n}^{-1})_{j,j}$. Thus, between
the linear and nonlinear kernels, the optimal one can be selected by
minimizing the LOOCV error.

##  4 Numerical Study

In Section 4.1, the sample paths of the functional-input GP with different
settings of the parameters are explored. In Section 4.2, numerical experiments
are conducted to examine the emulation performance of the proposed method.

###  4.1 Sample path

####  4.1.1 Linear kernel

The focus of this section is to study how the unknown parameters in the
proposed linear kernel (4) affect the generated sample paths. We focus on the
Matérn kernel function, which has the form of (9). Three types of unknown
parameters are studied, including the $d$ positive diagonal elements of the
diagonal matrix $\Theta$, denoted by $\boldsymbol{\theta}$, the positive
scalar $\sigma^{2}$, and the smoothness parameter $\nu$.

The sample paths are generated by the input functions $g(x)=\sin(\alpha x)$,
where $x\in\Omega=[0,2\pi]$ and $\alpha\in[0,1]$. The value $\alpha$ indicates
the frequency of the periodic function and the RKHS norm of $g$, i.e.,
$\|g\|_{\mathcal{N}_{\Psi}(\Omega)}$, increases monotonically with respect to
$\alpha$. As a result, this input function creates an analogy to the sample
paths in conventional GP by studying the paths as a function of $\alpha$ with
different parameter settings. In Figure 2, the sample paths are demonstrated
with different settings of the three types of parameters. The first row
illustrates the sample paths with three different settings of $\nu$, given
$\theta=1$ and $\sigma^{2}=1$. It appears that the smoothness of the resulting
sample paths is not significantly affected by the setting of $\nu$, which
typically controls the smoothness in conventional GPs. This is mainly because,
unlike the conventional Matérn kernel function (Stein,, 1999), the derivative
of (4) with respect to the input $g$ is not directly related to the parameter
$\nu$. The middle panels in Figure 2 demonstrate the sample paths with
different settings of $\theta$, given $\nu=2.5$ and $\sigma^{2}=1$. It shows
that, as $\theta$ increases, the number of the local maxima and minima
increases, which agrees with the observations in conventional GPs. Lastly, the
bottom three panels show the sample paths with different settings of
$\sigma^{2}$, given $\nu=2.5$ and $\theta=1$. Similar to conventional GPs,
$\sigma^{2}$ controls the amplitude of the resulting sample paths.

####  4.1.2 Nonlinear kernel

Based on the nonlinear kernel of (15) with the Matérn kernel function defined
in (10), the sample paths are studied with respect to different settings of
the $\gamma$, $\nu$ and $\sigma^{2}$. As shown in Figure 3, the results appear
to be consistent with the observations in conventional GPs where $\nu$
controls the smoothness of the function, $\theta$ controls the number of the
local maxima and minima, and $\sigma^{2}$ controls the amplitude of the
functions.

Figure 2: Sample paths of linear kernels. Top panel shows the effect of
varying the parameter $\nu$ with the fixed $\theta=1$ and $\sigma^{2}=1$,
middle panel shows the effect of varying the parameter $\theta$ with the fixed
$\nu=2.5$ and $\sigma^{2}=1$, and the bottom panel shows the effect of varying
the parameter $\sigma^{2}$ with the fixed $\nu=2.5$ and $\theta=1$. Figure 3:
Sample paths of nonlinear kernels. Top panel shows the effect of varying the
parameter $\nu$ with the fixed $\gamma=0.01$ and $\sigma^{2}=1$, middle panel
shows the effect of varying the parameter $\gamma$ with the fixed $\nu=2.5$
and $\sigma^{2}=1$, and the bottom panel shows the effect of varying the
parameter $\sigma^{2}$ with the fixed $\nu=2.5$ and $\gamma=0.01$.

###  4.2 Prediction performance

The prediction performance of the proposed method is examined by synthetic
examples in this section. Numerical results are simulated from the three
synthetic examples, including a linear operator,
$f_{1}(g)=\int_{\Omega}\int_{\Omega}g(\mathbf{x}){\rm{d}}x_{1}{\rm{d}}x_{2}$,
and two nonlinear ones,
$f_{2}(g)=\int_{\Omega}\int_{\Omega}g(\mathbf{x})^{2}{\rm{d}}x_{1}{\rm{d}}x_{2}$
and
$f_{3}(g)=\int_{\Omega}\int_{\Omega}\sin(g(\mathbf{x})){\rm{d}}x_{1}{\rm{d}}x_{2}$,
where $\mathbf{x}=(x_{1},x_{2})\in\Omega\equiv[0,1]^{2}$ and
$g(\mathbf{x}):[0,1]^{2}\rightarrow\mathbb{R}$. Consider eight functional
inputs, which are shown in the first row of Table 1, and assume that these
eight functional inputs are implemented for each of the synthetic examples and
their outputs are given in Table 1.

Three types of functional inputs are considered for predictions:
$g_{9}(\mathbf{x})=1+\sin(\alpha_{1}x_{1}+\alpha_{2}x_{2}),g_{10}(\mathbf{x})=\beta+x_{1}^{2}+x_{2}^{3}$,
and $g_{11}(\mathbf{x})=\exp\\{-\kappa x_{1}x_{2}\\}$, where
$\alpha_{1},\alpha_{2},\beta,\kappa\in[0,1]$. Based on 100 random samples of
$\alpha_{1},\alpha_{2},\beta$ and $\kappa$ from $[0,1]$, the prediction
performance is evaluated by averaging the mean absolute percentage errors
(MAPEs), where
$\text{MAPE}=\frac{1}{3}\sum^{11}_{j=9}\left|\frac{f(g_{j})-\hat{f}(g_{j})}{f(g_{j})}\right|\times
100\%.$

$g(\mathbf{x})$ |  $x_{1}+x_{2}$ |  $x_{1}^{2}$ |  $x_{2}^{2}$ |  $1+x_{1}$ |  $1+x_{2}$ |  $1+x_{1}x_{2}$ |  $\sin(x_{1})$ |  $\cos(x_{1}+x_{2})$
---|---|---|---|---|---|---|---|---
$f_{1}(g)$ |  1 |  0.33 |  0.33 |  1.5 |  1.5 |  1.25 |  0.46 |  0.50
$f_{2}(g)$ |  1.17 |  0.20 |  0.20 |  2.33 |  2.33 |  1.61 |  0.27 |  0.35
$f_{3}(g)$ |  0.77 |  0.31 |  0.31 |  0.96 |  0.96 |  0.93 |  0.43 |  0.45
Table 1: Training data set for the numerical study, where
$f_{1}(g)=\int_{\Omega}\int_{\Omega}g(\mathbf{x}){\rm{d}}x_{1}{\rm{d}}x_{2},f_{2}(g)=\int_{\Omega}\int_{\Omega}g(\mathbf{x})^{2}{\rm{d}}x_{1}{\rm{d}}x_{2}$,
and
$f_{3}(g)=\int_{\Omega}\int_{\Omega}\sin(g(\mathbf{x})){\rm{d}}x_{1}{\rm{d}}x_{2}$.

The proposed method is performed with the Matérn kernel function with the
smoothness parameter $\nu=5/2$, which leads to a simplified form of (10):

$\psi(r)=\left(1+\sqrt{5}r+\frac{5}{3}r^{2}\right)\exp\left(-5r\right).$ (19)

Other parameters, including $\boldsymbol{\theta},\tau$ and $\gamma$, are
estimated via maximum likelihood estimation. Both of the linear kernel (4) and
the nonlinear kernel (15) are performed for the proposed functional input GP
and their LOOCV errors are reported in Table 2. According to Section 3.3,
LOOCV is then used to identify the optimal kernel. By minimizing the LOOCV
errors, the linear kernel is identified as the optimal choice for the linear
synthetic example, $f_{1}(g)$, and the nonlinear kernel is identified as the
optimal choice for the nonlinear synthetic examples, $f_{2}(g)$ and
$f_{3}(g)$. For the three synthetic simulations, their MAPEs are summarized in
Table 2. It appears that the optimal kernels selected by LOOCV are consistent
with the selections based on minimizing the MAPEs, which shows that LOOCV is a
reasonable indicator of the optimal kernel when the ground truth is unknown.
Furthermore, based on the optimal kernels, the MAPEs (written in boldface) for
the three test functions are less than $6\%$, indicating a high accuracy in
predictions.

Measurements | Kernel | $f_{1}(g)=\int_{\Omega}\int_{\Omega}g$ | $f_{2}(g)=\int_{\Omega}\int_{\Omega}g^{2}$ | $f_{3}(g)=\int_{\Omega}\int_{\Omega}\sin(g)$
---|---|---|---|---
LOOCV | linear | $\boldsymbol{1.5\times 10^{-9}}$ | 0.1751 | 0.0308
nonlinear | $3.8\times 10^{-6}$ | 0.0210 | 0.0014
MAPE (%) | linear | 0.0006 | 36.88 | 16.30
nonlinear | 0.0219 | 5.70 | 3.21
Table 2: The mean absolute percentage errors (MAPEs) and the leave-one-out
cross-validation errors for three testing functions, in which the errors
corresponding to the optimal kernel are boldfaced.

##  5 Applications to Inverse Scattering Problem

In this section, we revisit the inverse scattering problem shown in Figure 4.
Let $D\subset\mathbb{R}^{2}$ denote an inhomogeneous isotropic scattering
region of interest, and the functional input $g$, the support of which is $D$,
is related to the refractive index for the region $D$ of the unknown
scatterer. Given a set of finite element simulations as the training data, the
goal of inverse scattering is to recover the functional input from a given
far-field pattern. To achieve this goal, an important task is to construct a
surrogate model for functional inputs.

In this study, 10 functional inputs, including
$1+x_{1},1-x_{1},1+x_{1}x_{2},1-x_{1}x_{2},1+x_{2},1-x_{2},1+x_{1}^{2},1-x_{1}^{2},1+x_{2}^{2}$,
and $1-x_{2}^{2}$, are conducted in the training set and the corresponding
far-field patterns are shown in Figure 4. The pre-processing of using
principal component analysis (PCA) is first applied to reduce the
dimensionality of the output images. The first three principal components,
denoted by $\mathbf{u}_{l}\in\mathbb{R}^{1024},l=1,2,3$, are shown in Figure
5, which explain more than 99.9% variability of the data. Therefore, given the
functional input $g_{i}$ for $i=1,\ldots,10$, the output of far-field images
can be approximated by $\sum^{3}_{l=1}f_{l}(g_{i})\mathbf{u}_{l}$, where
$f_{1}(g_{i}),f_{2}(g_{i}),f_{3}(g_{i})$ are the first three PC scores.

Figure 4: Training data. Figure 5: Principle components, which explain more
than 99.9% variations of the data.

After the dimension reduction, the three-dimensional outputs
$f_{1}(g),f_{2}(g)$ and $f_{3}(g)$ are assumed to be mutually independent and
follow the functional-input GP as the surrogate model. For any untried
functional input $g\in V$, by the results of (2) and (3), the far-field
pattern can by predicted by the following normal distribution:

$\mathcal{N}\left(\sum^{3}_{l=1}k_{l}(g)^{T}\mathbf{K}^{-1}_{n,l}(\mathbf{f}_{l}-\mu_{l}\mathbf{1}_{n})\mathbf{u}_{l},\sum^{3}_{l=1}(K_{l}(g,g)-k_{l}(g)^{T}\mathbf{K}^{-1}_{n,l}k_{l}(g))\mathbf{u}^{2}_{l}\right),$

where $\mathbf{f}_{l}=(f_{l}(g_{1}),\ldots,f_{l}(g_{n}))$,
$k_{l}(g)=(K_{l}(g,g_{1}),\ldots,K_{l}(g,g_{n}))^{T},(\mathbf{K}_{n,l})_{i,j}=K_{l}(g_{i},g_{j})$,
and $K_{l}$ is the kernel function with the hyperparameters estimated based on
$\mathbf{f}_{l}$.

The linear kernel of (4) and the nonlinear kernel of (15) are both performed.
The optimal kernel is selected by comparing the LOOCV errors in predicting the
far-field pattern. The LOOCV error based on the linear kernel is $9.6\times
10^{-8}$, which is slightly smaller than that of the nonlinear kernel,
$9.9\times 10^{-8}$. Therefore, we apply the linear kernel and examine its
prediction performance for two test functions, $g(\mathbf{x})=1-\sin(x_{2})$
and $g(\mathbf{x})=1$. The images of the true far-field patterns and their
predictions are illustrated in Figure 6. Compared to the ground truth, the
predictions can capture the underlying structures reasonably well with some
discrepancies appearing on the lower right corner when predicting the first
test function. The MAPEs in predicting the far-field patterns for these two
test functions are $1.10\%$ and $0.03\%$, respectively.

Figure 6: Prediction on the validation functions. The left two panels are the
true output and prediction for $g(x_{1},x_{2})=1-\sin(x_{2})$, and the right
two panels are the true output and prediction for $g(x_{1},x_{2})=1$.

##  6 Concluding Remarks

Although GP surrogates are widely used in the analysis of complex system as an
alternative to the direct analysis using computer experiments, most of the
existing work are not applicable to the problems with functional inputs. To
address this issue, two new types of kernel functions are introduced for
functional-input GPs, including a linear and nonlinear kernel. Theoretical
properties of the proposed surrogates, such as the convergence rate of the
mean squared prediction error, are discussed. Numerical studies and the
application to surrogate modeling in an inverse scattering problem show high
prediction accuracy of the proposed method.

There are extensive studies in experimental design for conventional GP
surrogates, but the study of optimal designs for GPs with functional inputs
remains scarce. In this paper, it is shown that space-fillingness is a
desirable property in controlling the convergence rate of MSPE. An interesting
topic for future research is to explore the construction procedure for
efficient space-filling designs with functional inputs. Apart from
experimental designs, another important future work is to explore systematic
approaches to efficiently identify the functional input given an observed far-
field pattern, which is the ultimate goal of inverse scattering problems.
Based on the proposed GP surrogate, we will explore a Bayesian inverse
framework and a calibration procedure that integrates computer experiments and
real observations.

Supplemental Materials Additional supporting materials can be found in
Supplemental Materials, including the theoretical proofs of Propositions 3.1
and 3.6, Theorems 3.3 and 3.7, and Corollaries 3.4 and 3.5, and the data, R
code, and workflow for reproducing the results in Sections 4 and 5.

## References

  * Cakoni et al., (2016) Cakoni, F., Colton, D., and Haddar, H. (2016).  Inverse Scattering Theory and Transmission Eigenvalues,CBMS-NSF Regional Conference Series in Applied Mathematics.  Philadelphia: SIAM.
  * Chen et al., (2011) Chen, D., Hall, P., and Müller, H.-G. (2011).  Single and multiple index functional regression models with nonparametric link.  The Annals of Statistics, 39(3):1720–1747.
  * Chen et al., (2021) Chen, J., Mak, S., Joseph, V. R., and Zhang, C. (2021).  Function-on-function kriging, with applications to three-dimensional printing of aortic tissues.  Technometrics, 63(3):384–395.
  * Currin et al., (1988) Currin, C., Mitchell, T., Morris, M., and Ylvisaker, D. (1988).  A bayesian approach to the design and analysis of computer experiments.  Technical report, Oak Ridge National Lab., TN (USA).
  * Fan and Zhang, (1999) Fan, J. and Zhang, W. (1999).  Statistical estimation in varying coefficient models.  The Annals of Statistics, 27(5):1491–1518.
  * Gramacy, (2020) Gramacy, R. B. (2020).  Surrogates: Gaussian Process Modeling, Design, and Optimization for the Applied Sciences.  CRC Press.
  * Hall and Horowitz, (2007) Hall, P. and Horowitz, J. L. (2007).  Methodology and convergence rates for functional linear regression.  The Annals of Statistics, 35(1):70–91.
  * Hastie et al., (2009) Hastie, T., Tibshirani, R., and Friedman, J. (2009).  The Elements of Statistical Learning: Data Mining, Inference, and Prediction.  Springer, second edition.
  * James, (2002) James, G. M. (2002).  Generalized linear models with functional predictors.  Journal of the Royal Statistical Society: Series B, 64(3):411–432.
  * Malfait and Ramsay, (2003) Malfait, N. and Ramsay, J. O. (2003).  The historical functional linear model.  Canadian Journal of Statistics, 31(2):115–128.
  * Morris, (2012) Morris, M. D. (2012).  Gaussian surrogates for computer models with time-varying inputs and outputs.  Technometrics, 54(1):42–50.
  * Mueller and Siltanen, (2020) Mueller, J. L. and Siltanen, S. (2020).  The D-bar method for electrical impedance tomography–demystified.  Inverse Problems, 36(4):093001.
  * Müller, (2009) Müller, S. (2009).  Komplexität und Stabilität von kernbasierten Rekonstruktionsmethoden.  PhD thesis, Niedersächsische Staats-und Universitätsbibliothek Göttingen.
  * Nguyen and Peraire, (2015) Nguyen, N. C. and Peraire, J. (2015).  Gaussian functional regression for linear partial differential equations.  Computer Methods in Applied Mechanics and Engineering, 287(15 April):69–89.
  * Rahman, (2007) Rahman, M. (2007).  Integral Equations and Their Applications.  WIT press.
  * Rasmussen and Williams, (2006) Rasmussen, C. E. and Williams, C. K. I. (2006).  Gaussian Processes for Machine Learning.  The MIT Press.
  * Ritter, (2007) Ritter, K. (2007).  Average-Case Analysis of Numerical Problems.  Springer.
  * Santner et al., (2018) Santner, T. J., Williams, B. J., and Notz, W. I. (2018).  The Design and Analysis of Computer Experiments.  Springer New York, second edition.
  * Stein, (1999) Stein, M. L. (1999).  Interpolation of Spatial Data: Some Theory for Kriging.  Springer Science & Business Media.
  * Steinwart and Scovel, (2012) Steinwart, I. and Scovel, C. (2012).  Mercer’s theorem on general domains: On the interaction between measures, kernels, and rkhss.  Constructive Approximation, 35(3):363–417.
  * Tan, (2019) Tan, M. H. (2019).  Gaussian process modeling of finite element models with functional inputs.  SIAM/ASA Journal on Uncertainty Quantification, 7(4):1133–1161.
  * Wang et al., (2010) Wang, S., Qian, L., and Carroll, R. J. (2010).  Generalized empirical likelihood methods for analyzing longitudinal data.  Biometrika, 97(1):79–93.
  * Wendland, (2004) Wendland, H. (2004).  Scattered Data Approximation.  Cambridge University Press.
  * Zhang and Wang, (2010) Zhang, H. and Wang, Y. (2010).  Kriging and cross-validation for massive spatial data.  Environmetrics: The official journal of the International Environmetrics Society, 21(3-4):290–304.

Supplementary Materials for “Functional-Input Gaussian Processes with
Applications to Inverse Scattering Problems”

##  S1 Proof of Proposition 3.1

For any non-zero $(\alpha_{1},...,\alpha_{n})$ and $(g_{1},...,g_{n})$, it
follows the quadratic form

$\displaystyle\sum_{j=1}^{n}\sum_{k=1}^{n}\alpha_{j}\alpha_{k}K(g_{j},g_{k})=\sum_{j=1}^{n}\sum_{k=1}^{n}\alpha_{j}\alpha_{k}\int_{\Omega}\int_{\Omega}g_{j}(\mathbf{x})g_{k}(\mathbf{x}^{\prime})\Psi(\mathbf{x},\mathbf{x}^{\prime}){\rm{d}}\mathbf{x}{\rm{d}}\mathbf{x}^{\prime}$
$\displaystyle=$
$\displaystyle\int_{\Omega}\int_{\Omega}\sum_{j=1}^{n}\sum_{k=1}^{n}\alpha_{j}\alpha_{k}g_{j}(\mathbf{x})g_{k}(\mathbf{x}^{\prime})\Psi(\mathbf{x},\mathbf{x}^{\prime}){\rm{d}}\mathbf{x}{\rm{d}}\mathbf{x}^{\prime}$
$\displaystyle=$
$\displaystyle\int_{\Omega}\int_{\Omega}\left(\sum_{j=1}^{n}\alpha_{j}g_{j}(\mathbf{x})\right)^{2}\Psi(\mathbf{x},\mathbf{x}^{\prime}){\rm{d}}\mathbf{x}{\rm{d}}\mathbf{x}^{\prime}\geq
0,$

and the quadratic form is strictly greater than zero if $g_{1},...,g_{n}$ are
linearly independent. This finishes the proof.

##  S2 Proof of Theorem 3.3

We need the following lemma, which is Proposition 10.28 of Wendland, (2004).

######  Lemma S2.1.

Suppose $\Psi$ is a symmetric and positive definite kernel on $\Omega$. Then
the integral operator $\mathcal{T}$ maps $L_{2}(\Omega)$ continuously into the
reproducing kernel Hilbert space $\mathcal{N}_{\Psi}(\Omega)$. It is the
adjoint of the embedding operator of the reproducing kernel Hilbert space
$\mathcal{N}_{\Psi}(\Omega)$ into $L_{2}(\Omega)$, i.e., it satisfies

$\displaystyle\langle g,v\rangle_{L_{2}(\Omega)}=\langle
g,\mathcal{T}v\rangle_{\mathcal{N}_{\Psi}(\Omega)},g\in\mathcal{N}_{\Psi}(\Omega),v\in
L_{2}(\Omega).$

Now we are ready to prove Theorem 3.3. For any
$\mathbf{u}=(u_{1},...,u_{n})^{T}\in\mathbb{R}^{n}$, direct computation shows

$\displaystyle\mathbb{E}\left(f(g)-\sum_{j=1}^{n}u_{j}f(g_{j})\right)^{2}=K(g,g)-2\sum_{j=1}^{n}u_{j}K(g,g_{j})+\sum_{j=1}^{n}\sum_{l=1}^{n}u_{j}u_{l}K(g_{j},g_{l})$
$\displaystyle=$
$\displaystyle\int_{\Omega}\int_{\Omega}g(\mathbf{x})g(\mathbf{x}^{\prime})\Psi(\mathbf{x},\mathbf{x}^{\prime}){\rm{d}}\mathbf{x}{\rm{d}}\mathbf{x}^{\prime}-2\sum_{j=1}^{n}u_{j}\int_{\Omega}\int_{\Omega}g(\mathbf{x})g_{j}(\mathbf{x}^{\prime})\Psi(\mathbf{x},\mathbf{x}^{\prime}){\rm{d}}\mathbf{x}{\rm{d}}\mathbf{x}^{\prime}$
$\displaystyle+\sum_{j=1}^{n}\sum_{l=1}^{n}u_{j}u_{l}\int_{\Omega}\int_{\Omega}g_{j}(\mathbf{x})g_{l}(\mathbf{x}^{\prime})\Psi(\mathbf{x},\mathbf{x}^{\prime}){\rm{d}}\mathbf{x}{\rm{d}}\mathbf{x}^{\prime}$
$\displaystyle=$ $\displaystyle\langle
g,\mathcal{T}g\rangle_{L_{2}(\Omega)}-2\sum_{j=1}^{n}u_{j}\langle
g,\mathcal{T}g_{j}\rangle_{L_{2}(\Omega)}+\sum_{j=1}^{n}\sum_{l=1}^{n}u_{j}u_{l}\langle
g_{j},\mathcal{T}g_{l}\rangle_{L_{2}(\Omega)}$ $\displaystyle=$
$\displaystyle\langle\mathcal{T}g,\mathcal{T}g\rangle_{\mathcal{N}_{\Psi}(\Omega)}-2\sum_{j=1}^{n}u_{j}\langle\mathcal{T}g,\mathcal{T}g_{j}\rangle_{\mathcal{N}_{\Psi}(\Omega)}+\sum_{j=1}^{n}\sum_{l=1}^{n}u_{j}u_{l}\langle\mathcal{T}g_{j},\mathcal{T}g_{l}\rangle_{\mathcal{N}_{\Psi}(\Omega)}$
$\displaystyle=$
$\displaystyle\left\|\mathcal{T}g-\sum_{j=1}^{n}u_{j}\mathcal{T}g_{j}\right\|_{\mathcal{N}_{\Psi}(\Omega)}^{2},$
(S2.1)

where the third equality is by Lemma S2.1.

Since $\hat{\mathbf{u}}:=\mathbf{K}_{n}^{-1}\mathbf{k}(g)$ minimizes the MSPE,
it also minimizes
$\left\|\mathcal{T}g-\sum_{j=1}^{n}u_{j}\mathcal{T}g_{j}\right\|_{\mathcal{N}_{\Psi}(\Omega)}^{2}$,
where $\mathbf{K}_{n}$ and $\mathbf{k}(g)$ are as in (2). By taking minimum on
both sides of (S2), we obtain

$\displaystyle\mathbb{E}\left(f(g)-\hat{f}(g)\right)^{2}=\min_{\mathbf{u}\in\mathbb{R}^{n}}\left\|\mathcal{T}g-\sum_{j=1}^{n}u_{j}\mathcal{T}g_{j}\right\|_{\mathcal{N}_{\Psi}(\Omega)}^{2},$
(S2.2)

which finishes the proof.

##  S3 Proof of Corollary 3.4

We need the following lemma, which states the asymptotic rates of the
eigenvalues of $K$. Lemma S3.1 is implied by the proof of Lemma 18 of Tuo and
Wang, (2020). In the rest of the Supplementary Materials, we will use the
following notation. For two positive sequences $a_{n}$ and $b_{n}$, we write
$a_{n}\asymp b_{n}$ if, for some constants $C,C^{\prime}>0$, $C\leq
a_{n}/b_{n}\leq C^{\prime}$. For notational simplicity, we will use
$C,C^{\prime},C_{1},C_{2},...$ to denote the constants, of which the values
can change from line to line.

######  Lemma S3.1.

Let $\Psi$ be a Matérn kernel function defined in (9), and
$\lambda_{1}\geq\lambda_{2}\geq...>0$ be its eigenvalues. Then,
$\lambda_{k}\asymp k^{-2\nu/d}$.

###### Proof of Corollary 3.4.

By (5), we have
$\mathcal{T}g(\mathbf{x})=\int_{\Omega}g(\mathbf{x}^{\prime})\Psi(\mathbf{x},\mathbf{x}^{\prime}){\rm
d}\mathbf{x}^{\prime}=\sum_{j=1}^{\infty}\lambda_{j}\langle
g,\phi_{j}\rangle_{L_{2}(\Omega)}\phi_{j}(\mathbf{x}).$ Therefore,

$\displaystyle\left\|\mathcal{T}g-\sum_{j=1}^{n}u_{j}\mathcal{T}g_{j}\right\|_{\mathcal{N}_{\Psi}(\Omega)}^{2}$
$\displaystyle=\left\|\sum_{j=1}^{\infty}\lambda_{j}\langle
g-\sum_{k=1}^{n}u_{k}g_{k},\phi_{j}\rangle_{L_{2}(\Omega)}\phi_{j}\right\|_{\mathcal{N}_{\Psi}(\Omega)}^{2}$
$\displaystyle=\sum_{j=1}^{\infty}\lambda_{j}\langle
g-\sum_{k=1}^{n}u_{k}g_{k},\phi_{j}\rangle_{L_{2}(\Omega)}^{2},$ (S3.3)

where the last equality holds by Theorem 10.29 of Wendland, (2004). Recall
that $g_{k}=\phi_{k}$ for $k=1,...,n$. Now take $u_{k}=\langle
g,\phi_{k}\rangle_{L_{2}(\Omega)}$ for $k=1,...,n$. It follows from Theorem
3.3 and (S3) that

$\displaystyle\mathbb{E}\left(f(g)-\hat{f}(g)\right)^{2}=\min_{\mathbf{u}\in\mathbb{R}^{n}}\left\|\mathcal{T}g-\sum_{j=1}^{n}u_{j}\mathcal{T}g_{j}\right\|_{\mathcal{N}_{\Psi}(\Omega)}^{2}$
$\displaystyle\leq$ $\displaystyle\sum_{j=1}^{\infty}\lambda_{j}\langle
g-\sum_{k=1}^{n}u_{k}g_{k},\phi_{j}\rangle_{L_{2}(\Omega)}^{2}=\sum_{j=n+1}^{\infty}\lambda_{j}\langle
g,\phi_{j}\rangle_{L_{2}(\Omega)}^{2}.$ (S3.4)

This indicates that the MSPE depends on the tail behavior of the summation
$\sum_{j=1}^{\infty}\lambda_{j}\langle g,\phi_{j}\rangle_{L_{2}(\Omega)}^{2}$.
Because $\Psi$ is a Matérn kernel function defined in (9), Lemma S3.1 implies
$\lambda_{j}\asymp j^{-\frac{2\nu}{d}}$. Then, an explicit convergence rate
can be obtained via

$\displaystyle\mathbb{E}\left(f(g)-\hat{f}(g)\right)^{2}\leq\sum_{j=n+1}^{\infty}\lambda_{j}\langle
g,\phi_{j}\rangle_{L_{2}(\Omega)}^{2}\leq\lambda_{n}\sum_{j=n+1}^{\infty}\langle
g,\phi_{j}\rangle_{L_{2}(\Omega)}^{2}\leq
C_{1}\|g\|_{L_{2}}^{2}n^{-\frac{2\nu}{d}},$

where we utilize $\sum_{j=n+1}^{\infty}\langle
g,\phi_{j}\rangle_{L_{2}(\Omega)}^{2}\leq\sum_{j=1}^{\infty}\langle
g,\phi_{j}\rangle_{L_{2}(\Omega)}^{2}=\|g\|_{L_{2}}^{2}$. This finishes the
proof of (11).

Next, we consider the case $g\in\mathcal{N}_{\Psi}(\Omega)$. Theorem 10.29 of
Wendland, (2004) yields that

$\displaystyle\|g\|_{\mathcal{N}_{\Psi}(\Omega)}^{2}=\sum_{j=1}^{\infty}\frac{\langle
g,\phi_{j}\rangle_{L_{2}(\Omega)}^{2}}{\lambda_{j}}<\infty.$

Then an alternative way to bound (S3) is by

$\displaystyle\mathbb{E}\left(f(g)-\hat{f}(g)\right)^{2}\leq\sum_{j=n+1}^{\infty}\lambda_{j}\langle
g,\phi_{j}\rangle_{L_{2}(\Omega)}^{2}=\sum_{j=n+1}^{\infty}\lambda_{j}^{2}\frac{\langle
g,\phi_{j}\rangle_{L_{2}(\Omega)}^{2}}{\lambda_{j}}$ $\displaystyle\leq$
$\displaystyle\lambda_{n}^{2}\sum_{j=n+1}^{\infty}\frac{\langle
g,\phi_{j}\rangle_{L_{2}(\Omega)}^{2}}{\lambda_{j}}\leq\lambda_{n}^{2}\|g\|_{\mathcal{N}_{\Psi}(\Omega)}^{2}\leq
C_{2}\|g\|_{\mathcal{N}_{\Psi}(\Omega)}^{2}n^{-\frac{4\nu}{d}}.$ (S3.5)

This finishes the proof of (12), and thus the proof of Corollary 3.4. ∎

##  S4 Proof of Corollary 3.5

The following lemma is used in the proof.

######  Lemma S4.1 (Wu and Schaback,, 1993, Theorem 5.14).

Let $\Omega$ be compact and convex with a positive Lebesgue measure;
$\Psi(\mathbf{x},\mathbf{x}^{\prime})$ be a Matérn kernel given by (9) with
the smoothness parameter $\nu$. Then there exist constants $c,c_{0}$ depending
only on $\Omega,\nu$ and the scale parameter $\gamma$ in (9), such that
$\Psi(\mathbf{x},\mathbf{x})-\mathbf{r}_{M}(\mathbf{x})^{T}\mathbf{R}_{M}^{-1}\mathbf{r}_{M}(\mathbf{x})\leq
ch_{\mathbf{X}_{n},\Omega}^{2\nu}$ provided that
$h_{\mathbf{X}_{n},\Omega}\leq c_{0}$, where
$\mathbf{R}_{M}=(\Psi(\mathbf{x}_{j},\mathbf{x}_{k}))_{jk}$ and
$\mathbf{r}_{M}(\mathbf{x})=(\Psi(\mathbf{x},\mathbf{x}_{1}),...,\Psi(\mathbf{x},\mathbf{x}_{n}))^{T}$.

###### Proof of Corollary 3.5.

For any $\mathbf{u}=(u_{1},...,u_{n})^{T}$, by (S3), we have

$\displaystyle\left\|\mathcal{T}g-\sum_{j=1}^{n}u_{j}\mathcal{T}g_{j}\right\|_{\mathcal{N}_{\Psi}(\Omega)}^{2}=\sum_{j=1}^{\infty}\lambda_{j}\langle
g-\sum_{k=1}^{n}u_{k}g_{k},\phi_{j}\rangle_{L_{2}(\Omega)}^{2}$
$\displaystyle=$ $\displaystyle\sum_{j=1}^{\infty}\lambda_{j}\langle
g-\sum_{k=1}^{n}u_{k}\Psi(\mathbf{x}_{k},\cdot),\phi_{j}\rangle_{L_{2}(\Omega)}^{2}\leq
C\sum_{j=1}^{\infty}\langle
g-\sum_{k=1}^{n}u_{k}\Psi(\mathbf{x}_{k},\cdot),\phi_{j}\rangle_{L_{2}(\Omega)}^{2}$
$\displaystyle=$ $\displaystyle
C\left\|g-\sum_{k=1}^{n}u_{k}\Psi(\mathbf{x}_{k},\cdot)\right\|_{L_{2}(\Omega)}^{2}\leq
C_{1}\sup_{\mathbf{x}\in\Omega}\left|g(\mathbf{x})-\sum_{k=1}^{n}u_{k}\Psi(\mathbf{x}_{k},\mathbf{x})\right|^{2},$
(S4.6)

where the last equality holds because $\phi_{j}$’s are orthogonal basis in
$L_{2}(\Omega)$. Therefore, we can take
$\mathbf{u}=(u_{1},...,u_{n})^{T}=\mathbf{R}^{-1}\mathbf{Y}$, where
$\mathbf{R}=(\Psi(\mathbf{x}_{j},\mathbf{x}_{k}))_{jk}$ and
$\mathbf{Y}=(g(\mathbf{x}_{1}),...,g(\mathbf{x}_{n}))^{T}$. Then the term
$\sup_{\mathbf{x}\in\Omega}\left|g(\mathbf{x})-\sum_{k=1}^{n}u_{k}\Psi(\mathbf{x}_{k},\mathbf{x})\right|$
becomes the prediction error of the radial basis function interpolation
(Wendland,, 2004), which is well understood in the literature.

For the completeness of this proof, we include the proof of an upper bound
here. Let
$\mathbf{v}=(v_{1},...,v_{n})^{T}=\mathbf{R}^{-1}\mathbf{r}(\mathbf{x})$,
where
$\mathbf{r}(\mathbf{x})=(\Psi(\mathbf{x},\mathbf{x}_{1}),...,\Psi(\mathbf{x},\mathbf{x}_{n}))^{T}$.
For any $\mathbf{x}\in\Omega$, the reproducing property implies that

$\displaystyle\left|g(\mathbf{x})-\sum_{k=1}^{n}u_{k}\Psi(\mathbf{x}_{k},\mathbf{x})\right|^{2}=\left|\langle
g,\Psi(\mathbf{x},\cdot)\rangle_{\mathcal{N}_{\Psi}(\Omega)}-\sum_{k=1}^{n}v_{k}(\mathbf{x})g(\mathbf{x}_{k})\right|^{2}$
$\displaystyle=$ $\displaystyle\left|\langle
g,\Psi(\mathbf{x},\cdot)\rangle_{\mathcal{N}_{\Psi}(\Omega)}-\sum_{k=1}^{n}v_{k}(\mathbf{x})\langle
g,\Psi(\mathbf{x}_{k},\cdot)\rangle_{\mathcal{N}_{\Psi}(\Omega)}\right|^{2}$
$\displaystyle\leq$
$\displaystyle\|g\|_{\mathcal{N}_{\Psi}(\Omega)}^{2}\left\|\Psi(\mathbf{x},\cdot)-\sum_{k=1}^{n}v_{k}(\mathbf{x})\Psi(\mathbf{x}_{k},\cdot)\right\|_{\mathcal{N}_{\Psi}(\Omega)}^{2}$
$\displaystyle=$
$\displaystyle\|g\|_{\mathcal{N}_{\Psi}(\Omega)}^{2}(\Psi(\mathbf{x},\mathbf{x})-\mathbf{r}(\mathbf{x})^{T}\mathbf{R}^{-1}\mathbf{r}(\mathbf{x})),$
(S4.7)

where the inequality is by the Cauchy-Schwarz inequality. A bound on
$\Psi(\mathbf{x},\mathbf{x})-\mathbf{r}(\mathbf{x})^{T}\mathbf{R}^{-1}\mathbf{r}(\mathbf{x})$
can be obtained via Lemma S4.1, which gives

$\displaystyle\Psi(\mathbf{x},\mathbf{x})-\mathbf{r}(\mathbf{x})^{T}\mathbf{R}^{-1}\mathbf{r}(\mathbf{x})\leq
C_{2}h_{\mathbf{X}_{n},\Omega}^{2\nu},$ (S4.8)

where $h_{\mathbf{X}_{n},\Omega}$ is the fill distance of $\mathbf{X}_{n}$.
Combining (8) with (S4), (S4) and (S4.8), we obtain that

$\displaystyle\mathbb{E}\left(f(g)-\hat{f}(g)\right)^{2}\leq
C_{3}h_{\mathbf{X}_{n},\Omega}^{2\nu},$

which finishes the proof. ∎

##  S5 Proof of Proposition 3.6

Without loss of generality, we assume $\gamma=1$. It suffices to show that
$\mathbf{R}=(K(g_{j},g_{k}))_{jk}$ is positive definite for any
$g_{1},...,g_{n}\in V$, which can be done by showing that there exist $n$
points $\mathbf{a}_{1},...,\mathbf{a}_{n}\in\mathbb{R}^{d}$ with $d<\infty$
such that
$\|\mathbf{a}_{j}-\mathbf{a}_{k}\|_{2}=\|g_{j}-g_{k}\|_{L_{2}(\Omega)}$. Let
$V_{n}={\rm span}(\\{g_{1},...,g_{n}\\})$, which is the linear space spanned
by $g_{1},...,g_{n}$. Clearly, $V_{n}$ is a finite dimensional space. Let
$\phi_{1},...,\phi_{m}$ be an orthogonal basis of $V_{n}$ with respect to
$L_{2}(\Omega)$, with $m\leq n$, which can be found via the Gram–Schmidt
process. Given this basis, we can write each $g_{j}$ as

$\displaystyle g_{j}=\sum_{k=1}^{m}a_{jk}\phi_{k},$

and let $\mathbf{a}_{j}=(a_{j1},...,a_{jm})^{T}\in\mathbb{R}^{m}$. Then it can
be checked that
$\|\mathbf{a}_{j}-\mathbf{a}_{k}\|_{2}=\|g_{j}-g_{k}\|_{L_{2}(\Omega)}$. Since
$\mathbf{R}=(K(g_{j}-g_{k}))_{jk}=(\psi(\|g_{j}-g_{k}\|_{L_{2}(\Omega)}))_{jk}=(\psi(\|\mathbf{a}_{j}-\mathbf{a}_{k}\|_{2}))_{jk}$,
which is positive definite, we finish the proof.

##  S6 Proof of Theorem 3.7

Before introducing the convergence results of MSPE of (16), we first provide a
characterization on the function class $V$. Since $\Phi$ is a positive
definite function, we can apply Mercer’s theorem to $\Phi$ and obtain

$\displaystyle\Phi(\mathbf{x},\mathbf{x}^{\prime})=\sum_{j=1}^{\infty}\lambda_{\Phi,j}\varphi_{j}(\mathbf{x})\varphi_{j}(\mathbf{x}^{\prime}),\quad\mathbf{x},\mathbf{x}^{\prime}\in\Omega,$
(S6.9)

where $\lambda_{\Phi,1}\geq\lambda_{\Phi,2}\geq...>0$ are eigenvalues and
$\\{\varphi_{k}\\}_{k\in\mathbb{N}}$ are eigenfunctions, and the summation is
uniformly and absolutely convergent. Because $g\in
V\subset\mathcal{N}_{\Phi}(\Omega)$, by Theorem 10.29 of Wendland, (2004), the
summation

$\displaystyle\|g\|_{\mathcal{N}_{\Phi}(\Omega)}^{2}=\sum_{j=1}^{\infty}\lambda_{\Phi,j}^{-1}\langle
g,\phi_{j}\rangle_{L_{2}(\Omega)}^{2}\leq 1.$

######  Definition S6.1.

A set $\Omega\subset{\mathbb{R}^{d}}$ is said to satisfy an interior cone
condition if there exists an angle $\alpha\in(0,\pi/2)$ and a radius
$\mathfrak{r}>0$ such that for every $\mathbf{x}\in\Omega$, a unit vector
$\xi(\mathbf{x})$ exists such that the cone

$C(\mathbf{x},\xi(\mathbf{x}),\alpha,\mathfrak{r}):=\left\\{\mathbf{x}+\eta\tilde{\mathbf{x}}:\tilde{\mathbf{x}}\in\mathbb{R}^{d},\|\tilde{\mathbf{x}}\|=1,\tilde{\mathbf{x}}^{T}\xi(\mathbf{x})\geq\cos(\alpha),\eta\in[0,\mathfrak{r}]\right\\}$

is contained in $\Omega$.

We need the following lemma, which is Theorem 11.8 of Wendland, (2004); also
see Theorem 3.14 of Wendland, (2004). Lemma S6.1 ensures the existence of the
local polynomial reproduction.

######  Lemma S6.1 (Local polynomial reproduction).

Let $l\in\mathbb{N}_{0}$ and $\pi_{l}(\mathbb{R}^{d})$ be the set of
$d$-variate polynomials with absolute degree no more than $l$. Suppose
$\Omega\subset\mathbb{R}^{d}$ is bounded and satisfies an interior cone
condition. Define

$\displaystyle
C_{1}=2,C_{2}=\frac{16(1+\sin(\alpha))^{2}l^{2}}{3\sin^{2}(\alpha)},c_{0}=\frac{\mathfrak{r}}{C_{2}},$

where $\alpha$ and $\mathfrak{r}$ are as in Definition S6.1. Then for all
$\mathbf{X}=\\{\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\\}\subset\Omega$ with
$h_{\mathbf{X},\Omega}\leq c_{0}$ and every $\mathbf{x}\in\Omega$ there exist
numbers $\tilde{u}_{1}(\mathbf{x}),\ldots,\tilde{u}_{n}(\mathbf{x})$ with

(1) $\sum_{j=1}^{n}\tilde{u}_{j}(\mathbf{x})p(\mathbf{x}_{j})=p(\mathbf{x})$
for all $p\in\pi_{l}(\mathbb{R}^{d})$,

(2) $\sum_{j=1}^{n}|\tilde{u}_{j}(\mathbf{x})|\leq C_{1}$,

(3) $\tilde{u}_{j}(\mathbf{x})=0$, if
$\|\mathbf{x}-\mathbf{x}_{j}\|>C_{2}h_{\mathbf{X},\Omega}$.

Lemma S6.2 comes from the proof of Theorem 11.9 and (11.6) of Wendland,
(2004), which provides an upper bound of the MSPE using $\tilde{u}_{j}$.

######  Lemma S6.2.

Suppose $\Phi=r(\|\cdot\|_{2})\in C^{k}(\mathbb{R}^{d})$ is positive definite.
Let $\Omega$ be a compact region satisfying an interior cone condition. Then
for $h_{\mathbf{X},\Omega}\leq h_{0}$,

$\displaystyle\Phi(\mathbf{x},\mathbf{x})-2\sum_{j=1}^{n}\tilde{u}_{j}\Phi(\mathbf{x},\mathbf{x}_{j})+\sum_{j=1}^{n}\sum_{k=1}^{n}\tilde{u}_{j}\tilde{u}_{k}\Phi(\mathbf{x}_{j},\mathbf{x}_{k})\leq(1+C_{1})^{2}\max_{0\leq
s\leq 2C_{2}h_{\mathbf{X},\Omega}}|\phi(s)-p(s^{2})|,$

where $p\in\pi_{\lfloor l/2\rfloor}(\mathbb{R})$.

###### Proof of Theorem 3.7.

Define a map $h:V\rightarrow W$ between the function class $V$ and the set
$W=\\{\mathbf{a}=(a_{1},...,a_{n},...)^{T}:\sum_{j=1}^{\infty}\lambda_{\Phi,j}^{-1}a_{j}^{2}\leq
1\\}\subset l_{2}(\mathbb{R}^{\infty})$ as

$\displaystyle h(g)=(\langle g,\phi_{1}\rangle_{L_{2}(\Omega)},...,\langle
g,\phi_{n}\rangle_{L_{2}(\Omega)},...)^{T}.$

It can be verified that $\|g\|_{L_{2}(\Omega)}=\|h(g)\|_{2}$. Therefore, we
can define a new positive definite function $\Psi_{1}$ on $W$ which satisfies

$\displaystyle\Psi_{1}(\mathbf{a},\mathbf{a}^{\prime})=\psi(\|\mathbf{a}-\mathbf{a}^{\prime}\|_{2})=\psi(\|g-g^{\prime}\|_{L_{2}(\Omega)})=\Psi(g,g^{\prime}),\forall
g,g^{\prime}\in V,$ (S6.10)

where $\mathbf{a}=h(g)$ and $\mathbf{a}^{\prime}=h(g^{\prime})$. Let
$\mathbf{a}^{(j)}=h(g_{j})$, and $\mathbf{a}=h(g)$. For any
$\mathbf{u}=(u_{1},...,u_{n})^{T}\in\mathbb{R}^{n}$, direct computation shows

$\displaystyle\mathbb{E}\left(f(g)-\sum_{j=1}^{n}u_{j}f(g_{j})\right)^{2}$
$\displaystyle=$ $\displaystyle
K(g,g)-2\sum_{j=1}^{n}u_{j}K(g,g_{j})+\sum_{j=1}^{n}\sum_{k=1}^{n}u_{j}u_{k}K(g_{j},g_{k})$
$\displaystyle=$
$\displaystyle\psi(0)-2\sum_{j=1}^{n}u_{j}\psi(\|\mathbf{a}-\mathbf{a}^{(j)}\|_{2})+\sum_{j=1}^{n}\sum_{k=1}^{n}u_{j}u_{l}\psi(\|\mathbf{a}^{(j)}-\mathbf{a}^{(k)}\|_{2}).$
(S6.11)

Let $\mathbf{b}_{j}=(a_{1}^{(j)},...,a_{m}^{(j)})^{T}$,
$\mathbf{b}=(a_{1},...,a_{m})^{T}$,
$\mathbf{a}_{c}^{(j)}=(a_{m+1}^{(j)},a_{m+2}^{(j)},...)$, and
$\mathbf{a}_{c}=(a_{m+1},a_{m+2},...)$, where $m$ will be determined later.
Then $\mathbf{a}^{(j)}=(\mathbf{b}_{j}^{T},(\mathbf{a}_{c}^{(j)})^{T})^{T}$
and $\mathbf{a}=(\mathbf{b}^{T},\mathbf{a}_{c}^{T})^{T}$. Applying Lemma S6.1
to (S6), we obtain that for some $\tilde{u}_{j}$,

$\displaystyle\sum_{j=1}^{n}\tilde{u}_{j}p(\mathbf{b}_{j})=p(\mathbf{b}),\mbox{
for all }p\in\pi_{l}(\mathbb{R}^{m}),$
$\displaystyle\sum_{j=1}^{n}|\tilde{u}_{j}(\mathbf{b})|\leq C_{1},\mbox{ and
}$ $\displaystyle\tilde{u}_{j}(\mathbf{b})=0,\mbox{ if
}\|\mathbf{b}-\mathbf{b}_{j}\|_{2}>C_{2}h_{\mathbf{B},\Omega},$ (S6.12)

when $h_{\mathbf{B},\Omega}\leq c_{0}$, where
$\mathbf{B}=\\{\mathbf{b}_{1},...,\mathbf{b}_{n}\\}$. Note that the
“constants” $C_{1},C_{2}$ and $c_{0}$ are related to the interior cone
condition and $l$. We need to compute these “constants” before providing an
upper bound. Here we use quote marks to denote that these constants are
related to the dimension $m$ of $\mathbf{b}$ and the degree $l$ of polynomials
$p$ in (S6.12). The constant $C_{1}=2$, thus we only need to consider $C_{2}$
and $c_{0}$.

Since $\mathbf{a},\mathbf{a}^{(j)}\in W$, we have that
$\mathbf{a}_{k}\leq\sqrt{\lambda_{\Phi,j}}$ and
$\mathbf{a}_{k}^{(j)}\leq\sqrt{\lambda_{\Phi,j}}$. Define a set
$V_{2}=\bigtimes_{k=1}^{m}[0,\sqrt{\lambda_{\Phi,k}}]$, it can be verified
that $\mathbf{b},\mathbf{b}_{j}\in V_{2}$. Set $\alpha=\pi/6$ and
$\mathfrak{r}=\sqrt{\lambda_{\Phi,m}}/2$. It can be verified that the interior
cone condition is satisfied. Then $C_{2}=48l^{2}$ and
$c_{0}=\frac{\sqrt{\lambda_{\Phi,m}}}{96l^{2}}$.

With $\tilde{u}_{j}$ defined in (S6.12), by (S6), it can be checked that

$\displaystyle\mathbb{E}\left(f(g)-\hat{f}(g)\right)^{2}$ $\displaystyle\leq$
$\displaystyle\psi(0)-2\sum_{j=1}^{n}\tilde{u}_{j}\psi(\|\mathbf{a}-\mathbf{a}^{(j)}\|_{2})+\sum_{j=1}^{n}\sum_{l=1}^{n}\tilde{u}_{j}\tilde{u}_{l}\psi(\|\mathbf{a}^{(j)}-\mathbf{a}^{(l)}\|_{2})$
$\displaystyle=$
$\displaystyle\left(\psi(0)-2\sum_{j=1}^{n}\tilde{u}_{j}\psi(\|\mathbf{b}-\mathbf{b}_{j}\|_{2})+\sum_{j=1}^{n}\sum_{k=1}^{n}\tilde{u}_{j}\tilde{u}_{k}\psi(\|\mathbf{b}_{j}-\mathbf{b}_{k}\|_{2})\right)$
$\displaystyle+\left(-2\sum_{j=1}^{n}\tilde{u}_{j}\left(\psi(\|\mathbf{a}-\mathbf{a}^{(j)}\|_{2})-\psi(\|\mathbf{b}-\mathbf{b}_{j}\|_{2})\right)+\sum_{j=1}^{n}\sum_{k=1}^{n}\tilde{u}_{j}\tilde{u}_{k}\left(\psi(\|\mathbf{a}^{(j)}-\mathbf{a}^{(l)}\|_{2})-\psi(\|\mathbf{b}_{j}-\mathbf{b}_{k}\|_{2})\right)\right)$
$\displaystyle=$ $\displaystyle I_{1}+I_{2}.$ (S6.13)

The first term can be bounded by Lemma S6.2, which gives us

$\displaystyle I_{1}\leq$ $\displaystyle 9\max_{0\leq s\leq
2C_{2}h_{\mathbf{B},\Omega}}|\psi(s)-p(s^{2})|=9\max_{0\leq
s\leq\sqrt{\lambda_{\Phi,m}}}|\psi(s)-p(s^{2})|,$ (S6.14)

for some $p\in\pi_{\lfloor l/2\rfloor}(\mathbb{R})$, provided
$h_{\mathbf{B},\Omega}\leq c_{0}$. Because
$\|\mathbf{b}_{j}-\mathbf{b}\|_{2}\leq\|\mathbf{a}-\mathbf{a}^{(j)}\|_{2}=\|g-g_{j}\|_{L_{2}(\Omega)}$,
it can be seen that $h_{\mathbf{B},\Omega}\leq h_{G_{n},V}$, and thus
$h_{\mathbf{B},\Omega}\leq c_{0}$ holds.

Next, we consider bounding $I_{1}$ with Matérn kernel functions. Lemma S3.1
implies that $\lambda_{\Phi,j}\asymp j^{-\frac{2\nu_{1}}{d}}$. By the
expansion of modified Bessel function (Abramowitz and Stegun,, 1948), $\psi$
can be written as

$\displaystyle\psi(r)=\sum_{k=0}^{\lfloor\nu\rfloor}c_{k}r^{2k}+c_{\psi}(r),$

where

$\displaystyle c_{\psi}(r)=\left\\{\begin{array}[]{cc}cr^{2\nu}\log
r+O(r^{2\nu})&\nu=1,2,...\\\ cr^{2\nu}+O(r^{2(\lfloor\nu\rfloor+1)})&\mbox{
otherwise.}\end{array}\right.$

Therefore, we can take $p(s^{2})=-\sum_{k=0}^{\lfloor\nu\rfloor}c_{k}s^{2k}$
and obtain that

$\displaystyle\max_{0\leq
s\leq\sqrt{\lambda_{\Phi,m}}}|\psi(s)-p(s^{2})|\leq\left\\{\begin{array}[]{ll}C_{2}\lambda_{\Phi,m}^{\nu}\log(\lambda_{\Phi,m}^{-1})\leq
C_{3}m^{-\frac{2\nu\nu_{1}}{d}}\log m,&\nu=2,...\\\
C_{4}\lambda_{\Phi,m}^{\nu}\leq
C_{5}m^{-\frac{2\nu\nu_{1}}{d}},&\mbox{otherwise.}\end{array}\right.$ (S6.17)

By (S6.14),

$\displaystyle I_{1}\leq 9\max_{0\leq
s\leq\sqrt{\lambda_{\Phi,m}}}|\psi(s)-p(s^{2})|\leq\left\\{\begin{array}[]{ll}C_{6}m^{-\frac{2\nu\nu_{1}}{d}}\log
m,&\nu=2,...\\\
C_{7}m^{-\frac{2\nu\nu_{1}}{d}},&\mbox{otherwise.}\end{array}\right.$ (S6.20)

It remains to bound $I_{2}$. For Matérn kernel functions, it can be checked
that for all $s_{1},s_{2}\in[0,s]$,

$\displaystyle|\psi(s_{1})-\psi(s_{2})|\leq C_{8}|s_{1}-s_{2}|^{2\tau},$

where $\tau=\min(\nu,1)$. Therefore, we can rewrite (S6) as

$\displaystyle|I_{2}|\leq$ $\displaystyle
2C_{8}\sum_{j=1}^{n}\tilde{u}_{j}|\|\mathbf{a}-\mathbf{a}^{(j)}\|_{2}-\|\mathbf{b}-\mathbf{b}_{j}\|_{2}|^{2\tau}+C_{8}\sum_{j=1}^{n}\sum_{k=1}^{n}\tilde{u}_{j}\tilde{u}_{k}|\|\mathbf{a}^{(j)}-\mathbf{a}^{(k)}\|_{2}-\|\mathbf{b}_{j}-\mathbf{b}_{k}\|_{2}|^{2\tau}$
$\displaystyle\leq$ $\displaystyle
2C_{8}\sum_{j=1}^{n}\left|\tilde{u}_{j}\right|\|\mathbf{a}-\mathbf{a}^{(j)}-(\mathbf{b}-\mathbf{b}_{j})\|_{2}^{2\tau}+C_{8}\sum_{j=1}^{n}\sum_{k=1}^{n}\left|\tilde{u}_{j}\right|\left|\tilde{u}_{l}\right|\|\mathbf{a}^{(j)}-\mathbf{a}^{(l)}-(\mathbf{b}_{j}-\mathbf{b}_{k})\|_{2}^{2\tau}$
$\displaystyle\leq$ $\displaystyle
C_{8}(4C_{1}+2C_{1}^{2})\lambda_{\Phi,m+1}^{\tau}\leq
C_{8}(4C_{1}+2C_{1}^{2})m^{-\frac{2\tau\nu_{1}}{d}}.$ (S6.21)

Because $m^{-\frac{2\nu\nu_{1}}{d}}\log m\leq m^{-\frac{2\nu_{1}\tau}{d}}\log
m$ and $\log m>1$, combining (S6.17), (S6) and (S6) leads to

$\displaystyle\mathbb{E}\left(f(g)-\hat{f}(g)\right)^{2}\leq
C_{9}m^{-\frac{2\nu_{1}\tau}{d}}\log m.$ (S6.22)

The last step is to compute $m$ such that there exist $n$ functions,
$h_{G_{n},V}\leq C_{0}m^{-\frac{2\nu_{1}}{d}}$. Since it is known that a unit
ball of $\mathcal{N}_{\Phi}(\Omega)$ has a covering number

$\displaystyle N(\delta,V,\|\cdot\|_{L_{\infty}})\leq
C_{10}\exp(C_{11}\delta^{-\frac{d}{\nu_{1}+d/2}}).$

Thus, in order to make $h_{G_{n},V}\leq C_{0}m^{-\frac{2\nu_{1}}{d}}$, we set
$\delta=C_{0}m^{-\frac{2\nu_{1}}{d}}$. Thus, as long as $n\geq
C_{12}\exp(C_{13}m^{\frac{2\nu_{1}}{\nu_{1}+d/2}})$, $h_{G_{n},V}\leq
C_{0}m^{-\frac{2\nu_{1}}{d}}$ holds. This implies $m\leq C_{14}(\log
n)^{\frac{\nu_{1}+d/2}{2\nu_{1}}}$. Plugging $m\leq C_{14}(\log
n)^{\frac{\nu_{1}+d/2}{2\nu_{1}}}$ in (S6.20), (S6), and (S6), we finish the
proof. ∎

## References

  * Tuo and Wang, (2020) Tuo, R. and Wang, W. (2020).  Kriging prediction with isotropic Matern correlations: robustness and experimental designs.  Journal of Machine Learning Research, 21(187):1–38.
  * Wendland, (2004) Wendland, H. (2004).  Scattered Data Approximation.  Cambridge University Press.
  * Wu and Schaback, (1993) Wu, Z. and Schaback, R. (1993).  Local error estimates for radial basis function interpolation of scattered data.  IMA Journal of Numerical Analysis, 13(1):13–27.
  * Abramowitz and Stegun, (1948) Abramowitz, M. and Stegun, I. A. (1948).  Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, volume 55.  US Government Printing Office.