# Optimal plug-in Gaussian processes for modelling derivatives

Zejian Liu<EMAIL_ADDRESS>Meng Li<EMAIL_ADDRESS>
###### Abstract

Derivatives are a key nonparametric functional in wide-ranging applications
where the rate of change of an unknown function is of interest. In the
Bayesian paradigm, Gaussian processes (GPs) are routinely used as a flexible
prior for unknown functions, and are arguably one of the most popular tools in
many areas. However, little is known about the optimal modelling strategy and
theoretical properties when using GPs for derivatives. In this article, we
study a plug-in strategy by differentiating the posterior distribution with GP
priors for derivatives of any order. This practically appealing plug-in GP
method has been previously perceived as suboptimal and degraded, but this is
not necessarily the case. We provide posterior contraction rates for plug-in
GPs and establish that they remarkably adapt to derivative orders. We show
that the posterior measure of the regression function and its derivatives,
with the same choice of hyperparameter that does not depend on the order of
derivatives, converges at the minimax optimal rate up to a logarithmic factor
for functions in certain classes. This to the best of our knowledge provides
the first positive result for plug-in GPs in the context of inferring
derivative functionals, and leads to a practically simple nonparametric
Bayesian method with guided hyperparameter tuning for simultaneously
estimating the regression function and its derivatives. Simulations show
competitive finite sample performance of the plug-in GP method. A climate
change application on analyzing the global sea-level rise is discussed.

##  1 Introduction

Consider the nonparametric regression model

$Y_{i}=f(X_{i})+\varepsilon_{i},\quad\varepsilon_{i}\sim N(0,\sigma^{2}),$ (1)

where the data $\mathbb{D}_{n}=\\{X_{i},Y_{i}\\}_{i=1}^{n}$ are independent
and identically distributed samples from a distribution $\mathbb{P}_{0}$ on
$\mathcal{X}\times\mathbb{R}$ that is determined by $\mathbb{P}_{X}$, $f_{0}$,
and $\sigma^{2}$, which are respectively the marginal distribution of $X_{i}$,
the true regression function, and the noise variance that is possibly unknown.
Let $p_{X}$ denote the density of $\mathbb{P}_{X}$ with respect to the
Lebesgue measure $\mu$. Here $\mathcal{X}\subset\mathbb{R}^{p}$ is a compact
metric space for $p\geq 1$.

We are interested in the inference on the derivative functions of $f$.
Derivatives emerge as a key nonparametric quantity when the rate of change of
an unknown surface is of interest. Examples include surface roughness for
digital terrain models, temperature or rainfall slope in meteorology, and
pollution curvature for environmental data. The importance of derivatives,
either as a nonparametric functional or localized characteristic of $f$, can
be also found in efficient modelling of functional data (Dai et al.,, 2018),
shape constrained function estimation (Riihimäki and Vehtari,, 2010; Wang and
Berger,, 2016), and the detection of stationary points (Yu et al.,, 2022).

Gaussian processes (GPs) are a popular nonparametric Bayesian method in many
areas such as spatially correlated data analysis (Stein,, 2012; Gelfand et
al.,, 2003; Banerjee et al.,, 2003), functional data analysis (Shi and Choi,,
2011), and machine learning (Rasmussen and Williams,, 2006); see also the
excellent review article by Gelfand and Schliep, (2016) which elaborates on
the instrumental role GPs have played as a key ingredient in an extensive list
of twenty years of modelling work. GPs not only provide a flexible process for
unknown functions but also serve as a building block in hierarchical models
for broader applications.

For function derivatives, the so-called plug-in strategy that directly
differentiates the posterior distribution of GP priors is practically
appealing, as it would allow users to employ the same prior no matter the
inference goal is on the regression function or its derivatives. However, this
plug-in estimator has been perceived as suboptimal and degraded for a decade
(Stein,, 2012; Holsclaw et al.,, 2013) based on heuristics, while a
theoretical understanding is lacking partly owing to technical challenges
posed by the irregularity and nonparametric nature of derivative functionals
(see Section 2.2 for more details). As a result, there is limited study of
plug-in GPs ever since, and substantially more complicated methods that hamper
easy implementation and often restrict to one particular derivative order are
pursued.

In this article, we study the plug-in strategy with GPs for derivative
functionals by characterizing large sample properties of the plug-in posterior
measure with GP priors, and obtain the first positive result. We show that the
plug-in posterior distribution, with the same choice of hyperparameter in the
GP prior, concentrates at the derivative functionals of any order at a nearly
minimax rate in specific examples, thus achieving a remarkable plug-in
property for nonparametric functionals that gains increasing attention
recently (Yoo and Ghosal,, 2016; Liu and Li,, 2022). It is known that many
commonly used nonparametric methods such as smoothing splines and local
polynomials do not enjoy this property when estimating derivatives (Wahba and
Wang,, 1990; Charnigo et al.,, 2011), and the only nonparametric Bayesian
method with established plug-in property, to the best of our knowledge, is
random series priors with B-splines (Yoo and Ghosal,, 2016).

In recent years, the nonparametric Bayesian literature has seen remarkable
adapativity of GP priors in various regression settings (van der Vaart and van
Zanten,, 2009; Bhattacharya et al.,, 2014; Jiang and Tokdar,, 2021). Our
findings contribute to this growing literature and indicate that the widely
used GP priors offer more than inferring regression functions. In particular,
the established theory reassures the use of plug-in GPs for optimal modelling
of derivatives, and further sheds lights into hyperparameter tuning in the
presence of varying derivative orders. That is, the Bayes procedure using GP
priors automatically adapts to the order of derivative, leading to a
practically simple nonparametric Bayesian method with guided hyperparameter
tuning for simultaneously estimating the regression function and its
derivatives. These theoretical guarantees are complemented by competitive
finite sample performance using simulations, as well as a climate change
application to analyzing the global sea-level rise.

The following notation is used throughout this paper. We write
$X=(X_{1}^{T},\ldots,X_{n}^{T})^{T}\in\mathbb{R}^{n\times p}$ and
$Y=(Y_{1},\ldots,Y_{n})^{T}\in\mathbb{R}^{n}$. Let $\|\cdot\|$ be the
Euclidean norm; for $f,g:\mathcal{X}\rightarrow\mathbb{R}$, let
$\|f\|_{\infty}$ be the $L_{\infty}$ (supremum) norm,
$\|f\|_{2}=(\int_{\mathcal{X}}f^{2}d\mathbb{P}_{X})^{1/2}$ the $L_{2}$ norm
with respect to the covariate distribution $\mathbb{P}_{X}$, and
$\left<f,g\right>_{2}=(\int_{\mathcal{X}}fgd\mathbb{P}_{X})^{1/2}$ the inner
product. The corresponding $L_{2}$ space relative to $\mathbb{P}_{X}$ is
denoted by $L^{2}_{p_{X}}(\mathcal{X})$; we write $L^{2}(\mathcal{X})$ as the
$L_{2}$ space with respect to the Lebesgue measure $\mu$. Denote the space of
all essentially bounded functions by $L^{\infty}(\mathcal{X})$. Let
$\mathbb{N}$ be the set of all positive integers and write
$\mathbb{N}_{0}=\mathbb{N}\cup\\{0\\}$. We let $C(\mathcal{X})$ and
$C(\mathcal{X},\mathcal{X})$ denote the space of continuous functions and
continuous bivariate functions. In one-dimensional case, for
$\Omega\subset\mathbb{R}$, a function $f:\Omega\rightarrow\mathbb{R}$ and
$k\in\mathbb{N}$, we use $f^{(k)}$ to denote its $k$-th derivative as long as
it exists and $f^{(0)}=f$. Let
$C^{m}(\Omega)=\\{f:\Omega\rightarrow\mathbb{R}\mid f^{(k)}\in C(\Omega)\ {\rm
for\ all\ }1\leq k\leq m\\}$ denote the space of $m$-times continuously
differentiable functions and
$C^{2m}(\Omega,\Omega)=\\{K:\Omega\times\Omega\rightarrow\mathbb{R}\mid\partial^{k}_{x}\partial^{k}_{x^{\prime}}K(x,x^{\prime})\in
C(\Omega,\Omega)\ {\rm for\ all\ }1\leq k\leq m\\}$ denote the space of
$m$-times continuously differentiable bivariate functions, where
$\partial^{k}_{x}=\partial^{k}/\partial x^{k}$. For two sequences $a_{n}$ and
$b_{n}$, we write $a_{n}\lesssim b_{n}$ if $a_{n}\leq Cb_{n}$ for a universal
constant $C>0$, and $a_{n}\asymp b_{n}$ if $a_{n}\lesssim b_{n}$ and
$b_{n}\lesssim a_{n}$.

##  2 Main results

###  2.1 Plug-in Gaussian process for derivative functionals

We assign a Gaussian process prior $\Pi$ on the regression function such that
$f\sim\mathrm{GP}(0,\sigma^{2}(n\lambda)^{-1}K)$. Here
$K(\cdot,\cdot):\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ is a
continuous, symmetric and positive definite kernel function, and $\lambda>0$
is a regularization parameter that possibly depends on the sample size $n$.
The rescaling factor $(n\lambda)^{-1}$ in the covariance kernel connects Bayes
estimator with kernel ridge regression (Wahba,, 1990; Cucker and Zhou,, 2007);
see also Theorem 11.61 in Ghosal and van der Vaart, (2017) for more discussion
on this connection.

It is not difficult to derive that the posterior distribution
$\Pi_{n}(\cdot\mid\mathbb{D}_{n})$ is also a GP:
$f|\mathbb{D}_{n}\sim\mathrm{GP}(\hat{f}_{n},\hat{V}_{n})$, where the
posterior mean $\hat{f}_{n}$ and posterior covariance $\hat{V}_{n}$ are given
by

$\displaystyle\hat{f}_{n}(x)$ $\displaystyle=K(x,X)[K(X,X)+n\lambda
I_{n}]^{-1}Y,$ (2) $\displaystyle\hat{V}_{n}(x,x^{\prime})$
$\displaystyle=\sigma^{2}(n\lambda)^{-1}\\{K(x,x^{\prime})-K(x,X)[K(X,X)+n\lambda
I_{n}]^{-1}K(X,x^{\prime})\\},$ (3)

for any $x,x^{\prime}\in\mathcal{X}$. Here $K(X,X)$ is the $n$ by $n$ matrix
$(K(X_{i},X_{j}))_{i,j=1}^{n}$, $K(x,X)$ is the 1 by $n$ vector
$(K(x,X_{i}))_{i=1}^{n}$, and $I_{n}$ is the $n$ by $n$ identity matrix.

We now define the plug-in Gaussian process for differential operators. For
simplicity we focus on the one-dimensional case where $\mathcal{X}=[0,1]$
throughout the paper, but remark that the studied plug-in strategy can be
extended to multivariate cases straightforwardly despite more complicated
notation for high-order derivatives.

For any $k\in\mathbb{N}$, define the $k$-th differential operator
$D^{k}:C^{k}[0,1]\rightarrow C[0,1]$ by $D^{k}(f)=f^{(k)}$. If $K\in
C^{2k}([0,1],[0,1])$, then the posterior distribution of the derivative
$f^{(k)}|\mathbb{D}_{n}$, denoted by $\Pi_{n,k}(\cdot\mid\mathbb{D}_{n})$, is
also a Gaussian process since differentiation is a linear operator. In
particular, $f^{(k)}|\mathbb{D}_{n}\sim{\rm
GP}(\hat{f}^{(k)}_{n},\tilde{V}^{k}_{n})$, where

$\displaystyle\hat{f}^{(k)}_{n}(x)$ $\displaystyle=K_{k0}(x,X)[K(X,X)+n\lambda
I_{n}]^{-1}Y,$ (4) $\displaystyle\tilde{V}^{k}_{n}(x,x^{\prime})$
$\displaystyle=\sigma^{2}(n\lambda)^{-1}\left\\{K_{kk}(x,x^{\prime})-K_{k0}(x,X)[K(X,X)+n\lambda
I_{n}]^{-1}K_{0k}(X,x^{\prime})\right\\},$ (5)

with $K_{k0}(x,X)=(\partial^{k}_{x}K(x,X_{i}))_{i=1}^{n}$ and
$K_{kk}(x,x^{\prime})=\partial^{k}_{x}\partial^{k}_{x^{\prime}}K(x,x^{\prime})$.
Then the nonparametric plug-in procedure for $D^{k}$ refers to using the plug-
in posterior measure $\Pi_{n,k}(\cdot\mid\mathbb{D}_{n})$ for inference on
$D^{k}(f)$.

The plug-in posterior measure $\Pi_{n,k}(\cdot\mid\mathbb{D}_{n})$ has a
closed-form expression with given $\lambda$ and $\sigma^{2}$, substantially
facilitating its implementation in practice. The plug-in strategy is
practically appealing but has been perceived as suboptimal for a decade
(Stein,, 2012; Holsclaw et al.,, 2013) based on heuristics. To the contrary,
we will establish optimality of plug-in GPs and uncover its adaptivity to
derivative orders. Before we move on to studying large sample properties of
$\Pi_{n,k}(\cdot\mid\mathbb{D}_{n})$, in the next section we first take a
detour to present technical challenges when studying derivative functionals
that hamper theoretical development for this problem.

###  2.2 Nonparametric plug-in property and technical challenges

The technical challenges to study derivative functionals are twofold: the
irregularity of function derivatives at fixed points and nonparametric
extension of such derivatives.

The first challenge is related to the “plug-in property” in the literature.
The plug-in property proposed by Bickel and Ritov, (2003) refers to the
phenomenon that a rate-optimal nonparametric estimator also efficiently
estimates some bounded linear functionals. A parallel concept has been studied
in the Bayesian paradigm relying on posterior distributions and posterior
contraction rates (Rivoirard and Rousseau,, 2012; Castillo and Nickl,, 2013;
Castillo and Rousseau,, 2015).

However, function derivatives may not fall into the classical plug-in property
framework. To see this, let $D_{t}=f^{\prime}(t)$ be a functional which maps
$f$ to its derivative at any fixed point $t\in[0,1]$. While it is easy to see
that $D_{t}$ is a linear functional, the following Proposition 1 (Conway,,
1994, page 13) shows that $D_{t}$ is not bounded.

######  Proposition 1.

Let $t\in[0,1]$ and define $D_{t}:C^{1}[0,1]\rightarrow\mathbb{R}$ by
$D_{t}(f)=f^{\prime}(t)$. Then, there is no bounded linear functional on
$L^{2}[0,1]$ that agrees with $D_{t}$ on $C^{1}[0,1]$.

Therefore, it appears difficult to analyze function derivatives evaluated at a
fixed point, as existing work on the plug-in property typically assumes the
functional to be bounded (Bickel and Ritov,, 2003; Castillo and Nickl,, 2013;
Castillo and Rousseau,, 2015).

The second challenge is linked to that the differential operator $D^{k}$
points to function-valued functionals, or nonparametric functionals, as
opposed to real-valued functionals studied in the classical plug-in property
literature. Hence, one needs to analyze the function-valued functionals
uniformly for all points in the support. To distinguish the plug-in property
for nonparametric functionals from its traditional counterpart for real-valued
functionals, we term this property as nonparametric plug-in property.

We overcome these challenges by resting on an operator-theoretic framework
(Smale and Zhou,, 2005, 2007), the equivalent kernel technique, and a recent
non-asymptotic analysis of nonparametric quantities (Liu and Li,, 2022), and
show that GPs enjoy the nonparametric plug-in property for differential
operators.

###  2.3 Posterior contraction for function derivatives

Throughout this article, we assume the true regression function $f_{0}\in
C^{k}[0,1]$ and the covariance kernel $K\in C^{2k}([0,1],[0,1])$. Let
$\\{\mu_{i}\\}_{i=1}^{\infty}$ and $\\{\phi_{i}\\}_{i=1}^{\infty}$ be the
eigenvalues and eigenfunctions of the kernel $K$ such that
$K(x,x^{\prime})=\sum_{i=1}^{\infty}\mu_{i}\phi_{i}(x)\phi_{i}(x^{\prime})$
for any $x,x^{\prime}\in[0,1]$, where the eigenvalues satisfy
$\mu_{1}\geq\mu_{2}\geq\cdots>0$ and $\mu_{i}\downarrow 0$, and eigenfunctions
form an orthonormal basis of $L_{p_{X}}^{2}[0,1]$. The existence of such
eigendecomposition is ensured by Mercer’s theorem. It can also be seen that
$\phi_{i}\in C^{k}[0,1]$ for all $i\in\mathbb{N}$ as $K\in
C^{2k}([0,1],[0,1])$.

We make the following assumptions on the eigenfunctions of the covariance
kernel.

######  Assumption (A).

There exists $C_{k,\phi}>0$ such that $\|\phi^{(k)}_{i}\|_{\infty}\leq
C_{k,\phi}i^{k}$ for any $i\in\mathbb{N}$.

######  Assumption (B).

There exists $L_{k,\phi}>0$ such that
$|\phi^{(k)}_{i}(x)-\phi^{(k)}_{i}(x^{\prime})|\leq
L_{k,\phi}i^{k+1}|x-x^{\prime}|$ for all $i\in\mathbb{N}$ and any
$x,x^{\prime}\in[0,1]$.

We will make extensive use of the so-called equivalent kernel $\tilde{K}$
(Rasmussen and Williams,, 2006, Chapter 7), which shares the same
eigenfunctions with $K$ with altered eigenvalues
$\nu_{i}={\mu_{i}}/{(\lambda+\mu_{i})}$ for $i\in\mathbb{N}$, i.e.,
$\tilde{K}(x,x^{\prime})=\sum_{i=1}^{\infty}\nu_{i}\phi_{i}(x)\phi_{i}(x^{\prime})$.
Note that $\tilde{K}$ is also a continuous, symmetric, and positive definite
kernel.

Under Assumption (A), we define an $m$-th order analog of effective dimension
of the kernel $K$ with respect to $p_{X}$ for $0\leq m\leq k$:

$\tilde{\kappa}_{m}^{2}=\sup_{x,x^{\prime}\in[0,1]}\partial_{x}^{m}\partial_{x^{\prime}}^{m}\tilde{K}(x,x^{\prime})=\sup_{x\in[0,1]}\sum_{i=1}^{\infty}\frac{\mu_{i}}{\lambda+\mu_{i}}\phi^{(m)}_{i}(x)^{2}.$
(6)

When $m=0$, the definition above recovers the effective dimension of $K$ in
the literature (Zhang,, 2005), namely,
$\tilde{\kappa}^{2}=\tilde{\kappa}_{0}^{2}=\sup_{x\in[0,1]}\tilde{K}(x,x)\lesssim\sum_{i=1}^{\infty}\frac{\mu_{i}}{\lambda+\mu_{i}}$.
We also define $\hat{\kappa}_{m}^{2}$ as a useful upper bound of
$\tilde{\kappa}_{m}^{2}$ up to a constant multiplier under Assumption (A) for
$0\leq m\leq k$:

$\hat{\kappa}_{m}^{2}=\sum_{i=1}^{\infty}\frac{i^{2m}\mu_{i}}{\lambda+\mu_{i}}\gtrsim\tilde{\kappa}_{m}^{2}.$
(7)

Note that $\tilde{\kappa}_{m}^{2}$ and $\hat{\kappa}_{m}^{2}$ depend on the
regularization parameter $\lambda$.

Let $\epsilon_{n}$ be a sequence such that $\epsilon_{n}\rightarrow 0$, and
$n\epsilon_{n}^{2}\rightarrow\infty$. For $\epsilon_{n}$ to become a posterior
contraction rate for $\Pi_{n,k}(\cdot\mid\mathbb{D}_{n})$ at $f_{0}^{(k)}$, we
would additionally need the following assumptions.

######  Assumption (C).

The regularization parameter $\lambda$ is chosen such that
$\tilde{\kappa}^{2}=o(\sqrt{n/\log n})$,
$\tilde{\kappa}_{k}^{2}=O(n\epsilon_{n}^{2}/\log n)$, and
$\hat{\kappa}_{k+1}^{2}=O(n)$.

######  Assumption (D).

$\epsilon_{n}$ is a non-asymptotic convergence rate of $\hat{f}_{n}^{(k)}$
under the $L_{2}$ or $L_{\infty}$ norm, i.e.,
$\|\hat{f}_{n}^{(k)}-f_{0}^{(k)}\|_{p}\lesssim\epsilon_{n}$ for $p=2$ or
$\infty$ with $\mathbb{P}_{0}^{(n)}$-probability tending to 1.

Assumptions (A), (B) and (C) are related to the eigendecomposition of the
covariance kernel, while Assumption (D) is on the convergence rate of the
derivatives of the posterior mean, which are well studied in nonparametric
statistics. Assumption (C) can be verified through direct calculations based
on the decay rate of eigenvalues; we provide the rates of
$\tilde{\kappa}^{2}$, $\tilde{\kappa}_{k}^{2}$ and $\hat{\kappa}_{k+1}^{2}$
for specific kernels later in the paper in Lemma 8 and Lemma 9. In view of the
connection between posterior mean and kernel ridge regression (KRR) estimator,
we can take advantage of the rich literature on KRR and kernel learning theory
to verify Assumption (D). Along this line, Section 2.5 will give the
convergence rates for special examples in Lemma 10 and Lemma 11.

The following theorem obtains the posterior contraction rates using plug-in
GPs for function derivatives.

######  Theorem 2.

Suppose that $f_{0}\in C^{k}[0,1]$, $K\in C^{2k}([0,1],[0,1])$ and
$\epsilon_{n}$ is a sequence such that $\epsilon_{n}\rightarrow 0$, and
$n\epsilon_{n}^{2}\rightarrow\infty$. Under Assumptions (A), (B), (C) and (D),
$\epsilon_{n}$ is a posterior contraction rate of
$\Pi_{n,k}(\cdot\mid\mathbb{D}_{n})$ at $f_{0}^{(k)}$, i.e., for any
$M_{n}\rightarrow\infty$,
$\Pi_{n,k}\left(\|f^{(k)}-f_{0}^{(k)}\|_{p}>M_{n}\epsilon_{n}\,\big{|}\,\mathbb{D}_{n}\right)\rightarrow
0$ in $\mathbb{P}_{0}^{(n)}$-probability.

Theorem 2 provides an approach for establishing posterior contraction rates.
Computing the contraction rate now boils down to analyzing the
eigendecomposition of covariance kernel (Assumptions (A), (B) and (C)) and the
convergence rate of the derivatives of the posterior mean (Assumption (D)). In
Section 2.5, we will apply Theorem 2 to derive the minimax optimal contraction
rate for estimating the derivatives of analytic-type functions (Theorem 4) and
Hölder smooth functions (Theorem 5). Although our focus is on function
derivatives, setting $k=0$ in Theorem 2 also provides contraction rates for
estimating the regression function.

###  2.4 Unknown error variance

In practice, the error variance $\sigma^{2}$ in the regression model is
unknown and needs to be estimated. Let $\sigma^{2}_{0}$ be its true value. van
der Vaart and van Zanten, (2009) proposed a fully Bayesian scheme by endowing
the standard error $\sigma$ with a hyperprior, which is supported on a compact
interval $[a,b]\subset(0,\infty)$ that contains $\sigma_{0}$ with a Lebesgue
density bounded away from zero. This approach has been followed by many others
(Bhattacharya et al.,, 2014; Li and Dunson,, 2020). de Jonge and van Zanten,
(2013) showed a Bernstein-von Mises theorem for the marginal posterior of
$\sigma$, where the prior for $\sigma$ is relaxed to be supported on
$(0,\infty)$. In other words, it is possible to simultaneously estimate the
regression function at an optimal nonparametric rate and the standard
deviation at a parametric rate.

Here we consider an empirical Bayes approach, which is widely used in practice
and eliminates the need of jointly sampling $f$ and $\sigma$. Under model (1),
the marginal likelihood is

$Y|X\sim N(0,\sigma^{2}(n\lambda)^{-1}K(X,X)+\sigma^{2}I_{n}).$ (8)

We estimate $\sigma^{2}$ by its maximum marginal likelihood estimator (MMLE)

$\hat{\sigma}^{2}_{n}=\lambda Y^{T}[K(X,X)+n\lambda I_{n}]^{-1}Y.$ (9)

Then we endow a new prior
$f\sim\mathrm{GP}(0,\hat{\sigma}^{2}_{n}(n\lambda)^{-1}K)$ by substituting
$\sigma^{2}$ in the original GP prior with $\hat{\sigma}^{2}_{n}$. The induced
posterior measure for $\Pi_{n,k}(\cdot\mid\mathbb{D}_{n})$ is denoted by
$\Pi_{n,k,\text{EB}}(\cdot\mid\mathbb{D}_{n})$.

The next theorem shows that the established theory holds under the empirical
Bayes scheme.

######  Theorem 3.

If $\lambda$ is chosen such that $n\lambda\rightarrow\infty$, then Theorem 2
holds under the empirical Bayes scheme, i.e., $\epsilon_{n}$ is a posterior
contraction rate of $\Pi_{n,k,\rm{EB}}(\cdot\mid\mathbb{D}_{n})$ at
$f_{0}^{(k)}$ under the same conditions.

###  2.5 Minimax optimality using special examples

We present two concrete examples as a direct application of Theorem 2, where
we use kernels with exponentially decaying and polynomially decaying
eigenvalues, respectively, for estimating analytic-type and Hölder smooth
functions. In this section, we consider a uniform sampling process for $p_{X}$
on $[0,1]$, and the corresponding probability measure $\mathbb{P}_{X}$ becomes
the Lebesgue measure.

We first consider the analytic-type function class $A^{\gamma_{0}}[0,1]$ for
the true regression function $f_{0}$:

$A^{\gamma_{0}}[0,1]=\left\\{f\in
L^{2}[0,1]:\|f\|^{2}_{A^{\gamma_{0}}[0,1]}=\sum_{i=1}^{\infty}e^{\gamma_{0}i}|f_{i}|<\infty,f_{i}=\left<f,\phi_{i}\right>_{2}\right\\},$
(10)

where the eigenfunctions are the Fourier basis functions

$\psi_{1}(x)=1,\ \psi_{2i}(x)=\sqrt{2}\cos(2\pi ix),\
\psi_{2i+1}=\sqrt{2}\sin(2\pi ix),\quad i\in\mathbb{N}.$ (11)

Such eigenfunctions easily satisfy Assumption (A) with
$C_{k,\psi}=\sqrt{2}(2\pi)^{k}$ and Assumption (B) with
$L_{k,\psi}=\sqrt{2}(2\pi)^{k+1}$. We use kernels with exponentially decaying
eigenvalues, i.e., the covariance kernel $K$ has an eigendecomposition
relative to $\mathbb{P}_{X}$ such that the eigenvalues satisfy

$\mu_{i}\asymp e^{-2\gamma i},\quad i\in\mathbb{N},$ (12)

for some $\gamma>0$. We denote such kernels by $K_{\gamma}$, which encompasses
the well-known squared exponential kernel in an approximate sense (Rasmussen
and Williams,, 2006; Pati and Bhattacharya,, 2015).

By verifying Assumptions (C) and (D) using our non-asymptotic analysis in
Section 3, Theorem 2 yields the nonparametric plug-in property for estimating
the derivatives of analytic-type functions: the posterior distribution
contracts at a nearly parametric rate under the $L_{2}$ norm.

######  Theorem 4.

Suppose $f_{0}\in A^{\gamma_{0}}[0,1]$ for $\gamma_{0}>\gamma>1/2$. If
$K_{\gamma}$ is used in the GP prior with the regularization parameter
$\lambda\asymp\log n/n$, then for any $k\in\mathbb{N}_{0}$, the posterior
distribution $\Pi_{n,k}(\cdot\mid\mathbb{D}_{n})$ contracts at $f_{0}^{(k)}$
at a nearly parametric rate $\epsilon_{n}=(\log n)^{k+1}/\sqrt{n}$ under the
$L_{2}$ norm.

Our second example considers the following Hölder space $H^{\alpha}[0,1]$ for
$f_{0}$:

$H^{\alpha}[0,1]=\left\\{f\in
L^{2}[0,1]:\|f\|^{2}_{H^{\alpha}[0,1]}=\sum_{i=1}^{\infty}i^{\alpha}|f_{i}|<\infty,f_{i}=\left<f,\psi_{i}\right>_{2}\right\\},$
(13)

for some $\alpha>0$, with the eigenfunctions being the Fourier basis
$\\{\psi_{i}\\}_{i=1}^{\infty}$. For any $f\in H^{\alpha}[0,1]$, $f$ has
continuous derivatives up to order $\lfloor\alpha\rfloor$ and the
$\lfloor\alpha\rfloor$th derivative is Lipschitz continuous of order
$\alpha-\lfloor\alpha\rfloor$.

We consider kernels with polynomially decaying eigenvalues $K_{\alpha}$, where

$\mu_{i}\asymp i^{-2\alpha},\quad i\in\mathbb{N}.$ (14)

Verifying Assumptions (C) and (D), and invoking Theorem 2, we obtain the
nonparametric plug-in property for estimating the derivatives of functions in
the Hölder class: the posterior distribution contracts at a nearly minimax
optimal rate under the $L_{2}$ norm (Stone,, 1982).

######  Theorem 5.

Suppose $f_{0}\in H^{\alpha}[0,1]$ for $\alpha>3/2$. If $K_{\alpha}$ is used
in the GP prior with the regularization parameter $\lambda\asymp({\log
n}/{n})^{\frac{2\alpha}{2\alpha+1}}$, then for any $k<\alpha-3/2$ and
$k\in\mathbb{N}_{0}$, the posterior distribution
$\Pi_{n,k}(\cdot\mid\mathbb{D}_{n})$ contracts at $f_{0}^{(k)}$ at the nearly
minimax optimal rate $\epsilon_{n}=(\log n/n)^{\frac{\alpha-k}{2\alpha+1}}$
under the $L_{2}$ norm.

Remark. Given the smoothness level of the regression function, the rate-
optimal estimation of $f_{0}$ and its derivatives are achieved under the same
choice of the regularization parameter $\lambda$ that does not depend on the
derivative order (Theorem 4 and Theorem 5). Therefore, the plug-in GP prior
automatically adapts to the order of the derivative to be estimated, leading
to easy parameter tuning. This nonparametric plug-in property in the context
of derivative estimation is not seen in some popular nonparametric methods,
including smoothing splines and local polynomial regression; see Liu and Li,
(2022) for more discussion. In the Bayesian paradigm, Yoo and Ghosal, (2016)
established this property for B-splines, and our work reassuringly suggests
that GP priors attain the same remarkable adaptivity. We compare the finite-
sample performance of these two adaptive methods in Section 4.

By verifying the condition $n\lambda\rightarrow\infty$ in Theorem 3, the
established results of posterior contraction in the above two examples hold
under the empirical Bayes scheme.

######  Corollary 6.

Theorem 4 and Theorem 5 hold under the empirical Bayes scheme with the
regularization parameter $\lambda\asymp\log n/n$ and $\lambda\asymp({\log
n}/{n})^{\frac{2\alpha}{2\alpha+1}}$, respectively.

##  3 Non-asymptotic analysis of key nonparametric quantities

In this section, we present the error rates for a list of key quantities in GP
regression that are useful for deriving our nonparametric plug-in theory and
verifying theorem assumptions. We defer the preliminaries for reproducing
kernel Hilbert space and other auxiliary technical results to Appendix B.1.

We first study the rates for the posterior variances of derivatives of
Gaussian processes, i.e., the variance of
$\Pi_{n,k}(\cdot\mid\mathbb{D}_{n})$; this error bound will be used in
establishing Theorem 2. Here we consider a general kernel $K\in
C^{2k}([0,1],[0,1])$. The posterior covariance of the $k$-th derivative of
$f$, denoted by $\tilde{V}^{k}_{n}(x,x^{\prime})$, is given in (5). By
comparing (3) and (5), we can equivalently rewrite
$\tilde{V}^{k}_{n}(x,x^{\prime})$ as

$\tilde{V}^{k}_{n}(x,x^{\prime})=\partial^{k}_{x}\partial^{k}_{x^{\prime}}\hat{V}_{n}(x,x^{\prime}).$
(15)

Therefore, the posterior covariance of derivative of GP is exactly the mixed
derivative of the posterior covariance of the original GP. This is expected as
the differential operator is linear. We write
$\tilde{V}^{k}_{n}(x)=\tilde{V}^{k}_{n}(x,x)$ for the posterior variance. We
caution that, however, $\tilde{V}^{k}_{n}(x)$ may not be obtained by taking
the derivatives of $\hat{V}_{n}(x)$ since differentiation and evaluation at
$x=x^{\prime}$ may not be exchangeable. The next lemma provides a non-
asymptotic error bound for $\tilde{V}^{k}_{n}(x)$.

######  Lemma 7.

Suppose $K\in C^{2k}([0,1],[0,1])$. Under Assumption (A), by choosing
$\lambda$ such that $\tilde{\kappa}^{2}=o(\sqrt{n/\log n})$, it holds with
$\mathbb{P}_{0}^{(n)}$-probability at least $1-n^{-10}$ that

$\|\tilde{V}^{k}_{n}\|_{\infty}\leq\frac{2\sigma^{2}\tilde{\kappa}_{k}^{2}}{n}.$
(16)

Next, we present the rates for the higher-order analog of the effective
dimension under specific examples, which help verify Assumption (C). We
consider $K_{\gamma}$ and $K_{\alpha}$ defined in Section 2.5 for the GP prior
and a uniform sampling process. We denote the equivalent kernel of
$K_{\gamma}$ and $K_{\alpha}$ by $\tilde{K}_{\gamma}$ and
$\tilde{K}_{\alpha}$, respectively.

Let the higher-order analog of the effective dimension for $K_{\gamma}$ and
$K_{\alpha}$ be $\tilde{\kappa}_{\gamma,m}^{2}$ and
$\tilde{\kappa}_{\alpha,m}^{2}$, respectively, where $m\in\mathbb{N}_{0}$ and
the subscripts $\gamma$ and $\alpha$ emphasize the use of specific kernels
compared to the general definition in (6). Note that we allow $m=0$, which
corresponds to $\tilde{\kappa}_{\gamma,0}^{2}=\tilde{\kappa}_{\gamma}^{2}$ and
$\tilde{\kappa}_{\alpha,0}^{2}=\tilde{\kappa}_{\alpha}^{2}$. The same
convention in notation applies to $\hat{\kappa}_{\gamma,m}^{2}$ and
$\hat{\kappa}_{\alpha,m}^{2}$. Lemma 8 and Lemma 9 provide the exact order of
$\tilde{\kappa}_{\gamma,m}^{2}$, $\hat{\kappa}_{\gamma,m}^{2}$ and
$\tilde{\kappa}_{\alpha,m}^{2}$, $\hat{\kappa}_{\alpha,m}^{2}$ with respect to
the regularization parameter $\lambda$.

######  Lemma 8.

$\tilde{K}_{\gamma}\in C^{2m}([0,1]\times[0,1])$ and
$\tilde{\kappa}_{\gamma,m}^{2}\asymp\hat{\kappa}_{\gamma,m}^{2}\asymp(-\log\lambda)^{2m+1}$
for any $\gamma>0$, $0<\lambda<1$ and $m\in\mathbb{N}_{0}$.

######  Lemma 9.

$\tilde{K}_{\alpha}\in C^{2m}([0,1]\times[0,1])$ and
$\tilde{\kappa}_{\alpha,m}^{2}\asymp\hat{\kappa}_{\alpha,m}^{2}\asymp\lambda^{-\frac{2m+1}{2\alpha}}$
for any $\alpha>m+1/2$ and $m\in\mathbb{N}_{0}$.

Finally, we derive non-asymptotic convergence rates of $\hat{f}_{n}^{(k)}$ for
analytic-type functions and Hölder smooth functions considered in Section 2.5.
These results facilitate justifying Assumption (D).

######  Lemma 10.

Suppose $f_{0}\in A^{\gamma_{0}}[0,1]$ for $\gamma_{0}>\gamma>1/2$. If
$K_{\gamma}$ is used in the GP prior, then for any $k\in\mathbb{N}_{0}$ it
holds with $\mathbb{P}_{0}^{(n)}$-probability at least $1-n^{-10}$ that

$\|\hat{f}_{n}^{(k)}-f_{0}^{(k)}\|_{2}\lesssim\frac{(\log
n)^{k+1}}{\sqrt{n}},$ (17)

with the corresponding choice of regularization parameter $\lambda\asymp 1/n$.

######  Lemma 11.

Suppose $f_{0}\in H^{\alpha}[0,1]$ for $\alpha>k+1/2$ and
$k\in\mathbb{N}_{0}$. If $K_{\alpha}$ is used in the GP prior, then it holds
with $\mathbb{P}_{0}^{(n)}$-probability at least $1-n^{-10}$ that

$\|\hat{f}_{n}^{(k)}-f_{0}^{(k)}\|_{2}\lesssim\left(\frac{\log
n}{n}\right)^{\frac{\alpha-k}{2\alpha+1}},$ (18)

with the corresponding choice of regularization parameter $\lambda\asymp({\log
n}/{n})^{\frac{2\alpha}{2\alpha+1}}$.

##  4 Simulation

We carry out a simulation study to assess the finite sample performance of the
proposed nonparametric plug-in procedure for function derivatives. We consider
the true function $f_{0}(x)=\sqrt{2}\sum_{i=1}^{\infty}i^{-4}\sin
i\cos[(i-1/2)\pi x]$, $x\in[0,1]$, which has Hölder smoothness level
$\alpha=3$. We simulate $n$ observations from the regression model
$Y_{i}=f_{0}(X_{i})+\varepsilon_{i}$ with $\varepsilon_{i}\sim N(0,0.1)$ and
$X_{i}\sim\text{Unif}[0,1]$. We consider three sample sizes 100, 500, and
1000, and replicate the simulation 100 times.

For Gaussian process priors, we consider three commonly used covariance
kernels: the Matérn kernel, squared exponential (SE) kernel, and second-order
Sobolev kernel, which are given by

$\displaystyle K_{\text{Mat},\nu}(x,x^{\prime})$
$\displaystyle=\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\sqrt{2\nu}|x-x^{\prime}|\right)^{\nu}B_{\nu}\left(\sqrt{2\nu}|x-x^{\prime}|\right),$
(19) $\displaystyle K_{\text{SE}}(x,x^{\prime})$
$\displaystyle=\exp(-(x-x^{\prime})^{2}),$ (20) $\displaystyle
K_{\text{Sob}}(x,x^{\prime})$
$\displaystyle=1+xx^{\prime}+\min\\{x,x^{\prime}\\}^{2}(3\max\\{x,x^{\prime}\\}-\min\\{x,x^{\prime}\\})/6.$
(21)

Here $B_{\nu}(\cdot)$ is the modified Bessel function of the second kind with
$\nu$ being the smoothness parameter to be determined. For the Matérn kernel,
it is well known that the eigenvalues of $K_{\text{Mat},\nu}$ decay at a
polynomial rate, that is, $\mu_{i}\asymp i^{-2(\nu+1/2)}$ for
$i\in\mathbb{N}$.

We compare various Gaussian process priors with a random series prior using
B-splines. B-splines are widely used in nonparametric regression (James et
al.,, 2009; Wang et al.,, 2022). In the context of estimating function
derivatives, the B-spline prior with normal basis coefficients has been
recently shown to enjoy the plug-in property (cf. Theorem 4.2 in Yoo and
Ghosal, (2016)), which to the best of our knowledge constitutes the only
Bayesian method in the existing literature with this property. The
implementation of this B-spline prior follows Yoo and Ghosal, (2016). In
particular, for any $x\in[0,1]$, let $b_{J,4}(x)=(B_{j,4}(x))_{j=1}^{J}$ be a
B-spline of order $4$ and degrees of freedom $J$ with uniform knots. The prior
on $f$ is given as $f(x)=b_{J,4}(x)^{T}\beta$ with each entry of $\beta$
following $N(0,\sigma^{2})$ independently. The unknown variance $\sigma^{2}$
is estimated by its MMLE
$\hat{\sigma}_{n}^{2}=n^{-1}Y^{T}(BB^{T}+I_{n})^{-1}Y$, where
$B=(b_{J,4}(X_{1}),\ldots,b_{J,4}(X_{n}))^{T}$. The number of interior knots
$N=J-4$ is selected from $\\{1,2,\ldots,10\\}$ using leave-one-out cross
validation.

For Gaussian process priors, we adopt the same strategy of leave-one-out cross
validation to select the degree of freedom $\nu$ in $K_{\text{Mat},\nu}$ for a
fair comparison. In particular, we consider $\nu$ from the set {2, 2.5, 3,
3.5, …, 10}. The regularization parameter $\lambda$ and unknown $\sigma^{2}$
are determined by empirical Bayes through maximizing the marginal likelihood.
In addition to the estimation using each kernel, we also automatically select
the best kernel among the three via leave-one-out cross validation.

For each method, we evaluate the posterior mean $\hat{f}_{n}$ and
$\hat{f}_{n}^{\prime}$ at 100 equally spaced points in $[0,1]$, and calculate
the root mean square error (RMSE) between the estimates and the true
functions:

$\mathrm{RMSE}=\sqrt{\frac{1}{100}\sum_{t=0}^{99}\\{\hat{s}(t/99)-s(t/99)\\}^{2}},$
(22)

where $\hat{s}$ is the estimated function ($\hat{f}_{n}$ or
$\hat{f}_{n}^{\prime}$) and $s$ is the true function ($f_{0}$ or
$f_{0}^{\prime}$).

Table 1 reports the average RMSE of all methods over 100 repetitions for
$f_{0}$ and $f_{0}^{\prime}$. Clearly the RMSE of all methods steadily
decreases as the sample size $n$ increases. The squared exponential kernel and
Matérn kernel are the two leading approaches for all sample sizes and for both
$f_{0}$ and $f_{0}^{\prime}$. While the difference between various methods for
$f_{0}$ tends to vanish when the sample size increases to $n=1000$, the
performance gap in estimating $f_{0}^{\prime}$ is more profound. In
particular, compared to the squared exponential kernel, the Sobolev kernel
increases the average RMSE by nearly 60% from 1.53 to 2.41, while the increase
becomes more than twofold for the B-spline method (from 1.53 to 3.47). We
notice considerably large RMSEs for B-splines in a proportion of simulations,
so we also calculate the median RMSEs for better robustness, which are 0.08,
0.03, 0.02, 0.51, 0.31 and 0.24 for the six scenarios. This improves the
summarized RMSEs to close to the Sobolev kernel for both $f_{0}$ and
$f_{0}^{\prime}$ when $n=1000$. RMSEs of GP priors with the selected kernel
via cross validation do not significantly deviate from the best performing
kernel relative to the standard errors, suggesting that it can be a useful
strategy to choose kernels if desired.

Table 1: RMSE of estimating $f_{0}$ and $f_{0}^{\prime}$, averaged over 100 repetitions. The first four rows are the plug-in GP prior with various kernels (Matérn kernel, squared exponential kernel, second-order Sobolev kernel, and the selected kernel via cross validation). The last row is the random series prior using B-splines. Standard errors are provided in parentheses.  | $f_{0}$ | $f_{0}^{\prime}$
---|---|---
| $n=100$ | 500 | 1000 | $n=100$ | 500 | 1000
Matérn | 0.05 (0.02) | 0.02 (0.01) | 0.02 (0.01) | 0.27 (0.11) | 0.17 (0.04) | 0.19 (0.11)
SE | 0.05 (0.02) | 0.02 (0.01) | 0.02 (0.01) | 0.23 (0.10) | 0.15 (0.04) | 0.15 (0.03)
Sobolev | 0.05 (0.02) | 0.03 (0.01) | 0.03 (0.01) | 0.31 (0.07) | 0.25 (0.04) | 0.24 (0.02)
CV | 0.05 (0.02) | 0.02 (0.01) | 0.02 (0.01) | 0.30 (0.08) | 0.17 (0.05) | 0.16 (0.13)
B-splines | 0.08 (0.03) | 0.03 (0.01) | 0.02 (0.01) | 0.57 (0.36) | 0.49 (0.58) | 0.33 (0.28)

We next choose two representative examples from the 100 repetitions to
visualize the estimates of $f_{0}$ and $f_{0}^{\prime}$ in Figure 1. The
dotted line stands for the posterior mean $\hat{f}_{n}$ and dashes lines for
the 95% simultaneous $L_{\infty}$ credible bands, where the radius is
estimated by the 95% quantile of the posterior samples of
$\|f-\hat{f}_{n}\|_{\infty}$ and
$\|f^{\prime}-\hat{f}_{n}^{\prime}\|_{\infty}$. For the B-spline prior, we use
the default setting as in Yoo and Ghosal, (2016) by specifying the inflation
factor $\rho=0.5$. Note that these credible bands have fixed width by
construction.

|  |  Matérn |  SE |  Sobolev |  B-splines
---|---|---|---|---|---
Example 1 |  $f_{0}$ |  |  |  |
|  $f_{0}^{\prime}$ |  |  |  |
Example 2 |  $f_{0}$ |  |  |  |
|  $f_{0}^{\prime}$ |  |  |  |
Figure 1: Visualization of estimates of $f_{0}$ and $f_{0}^{\prime}$ in two
examples with $n=1000$. Dots: posterior mean $\hat{f}_{n}$; Solid: true
function $f_{0}$ or $f_{0}^{\prime}$; Dashes: 95% $L_{\infty}$-credible bands.

The first and third rows of Figure 1 show that the four methods lead to
comparable estimation and uncertainty quantification when estimating $f_{0}$.
However, the deviation between different methods is considerably widened for
the estimation of $f_{0}^{\prime}$. Matérn and squared exponential kernels
constantly give the most accurate point estimation, and their credible bands
cover the ground truth with reasonable width, indicating the effectiveness of
the nonparametric plug-in procedure using GP priors. There is a tendency for
the second-order Sobolev kernel to fail to capture $f_{0}^{\prime}$ around the
left endpoint, particularly in Example 1. The performance of the B-spline
method continues to exhibit sensitivity to the choice of $N$, selected by
leave-one-out cross validation. In Example 1, the selected $N$ is 1, and the
B-spline prior yields comparable credible bands of $f_{0}^{\prime}$ than GP
priors with slightly altered estimation near zero; in Example 2 shown in the
fourth row, the B-spline method with $N=5$ gives a point estimate that is
substantially worse than the other three GP methods, and the associated
credible bands are off the chart. We remark that the performance of B-splines
might be substantially improved had the number of knots been selected by a
different tuning method or with a different simulation setting. While leave-
one-out cross validation may be appropriate for selecting $N$ when estimating
$f_{0}$, as observed in Yoo and Ghosal, (2016), our results suggest that
adjustments or alternative strategies seem to be needed when the objective is
to make inference on $f_{0}^{\prime}$. For GP priors, our numerical results
suggest choosing $\lambda$ using empirical Bayes seems to be a reasonable
strategy for both $f_{0}$ and $f_{0}^{\prime}$, which is in line with the
nonparametric plug-in property of GP priors.

In another simulation reported in Appendix D, we compare our plug-in GP
estimators with the inverse method proposed in Holsclaw et al., (2013), where
the authors perceived the plug-in estimator as suboptimal. We consider the
regression function $f_{0}(x)=x\sin(x)/10$ given by one simulated example in
Holsclaw et al., (2013); the implementation of the inverse method follows the
authors’ setup and uses their published code online. We generate $n=100$ and
500 data points on a regular grid in $[0,10]$, and generate noise
$\varepsilon_{i}\sim N(0,\sigma^{2})$ with $\sigma=0.1$, 0.2 and 0.3. The goal
is to estimate $f_{0}^{\prime}$. We have found that the inverse method has a
noticeable frequency to produce a zero estimate and the estimates have
extremely high variability across replications, partly due to complicated and
poor hyperparameters tuning. We compare its performance with and without these
simulations in which the inverse method is not stable. The proposed plug-in GP
method continues to give competitive performance, and we have found no
scenario in the considered settings with varying sample size and noise
standard deviation where the inverse method has significantly smaller RMSE
relative to standard errors. In contrast, plug-in GPs are significantly better
when the sample size is $n=500$ and noise standard derivation $\sigma=0.1$. We
also note that the inverse method is not only computationally intensive but
also restricts to one particular derivative order, and its generalization to
other derivative orders is nontrivial. The proposed plug-in GP is instead
computationally efficient and applicable for general derivative orders.

In addition to the empirical Bayes approach to choose $\sigma^{2}$ and
$\lambda$, in Appendix D we also implement a fully Bayesian alternative where
inverse Gamma and Gamma priors are assigned to $\sigma^{2}$ and $\lambda$,
respectively. We do not observe any significant difference between the results
of the two treatments under the considered simulation setting.

##  5 Real data application

We apply the proposed plug-in GP method to analyze the rate of global sea rise
using global mean sea-level (GMSL) records from coastal and island tide-gauge
measurements from the late 19th to early 21st century. This global tide-gauge
record dataset has been described in detail and previously analyzed in Church
and White, (2006, 2011). Our analysis uses three variables: time ($x$) in
years AD, GMSL ($y$) in millimeters, and one-sigma sea-level observational
error ($\sigma_{y}$) in millimetres. The sample size is $n=130$. We are
interested in estimating the rate of GMSL rise, denoted by $f^{\prime}$ with
$f$ being the regression function of $y$ on $x$.

As a descriptive summary, the total GMSL rise from January 1880 to December
2009 is about 210 mm over the 130 years. The least-squares linear trend of
sea-level rise from 1900 AD to 2009 AD is 1.7 mm/yr $\pm$ 0.3 mm/yr (Church
and White,, 2006). Still using linear regression models, Church and White,
(2011) reduced the uncertainty and obtained an estimate of 1.7 mm/yr $\pm$ 0.2
mm/yr for the same period. Nonparametric inference on $f^{\prime}(x)$ as an
unknown function allows us to depict flexible rate changes that are possibly
time-varying.

We assume a nonparametric regression model for the observed data
$y_{i}=f(x_{i})+\varepsilon_{i}$, for $i=1880,1881,\ldots,2009$. We consider
two error structures: (i) a homogeneous error structure with
$\varepsilon_{i}\sim N(0,\sigma^{2})$, and (ii) a heterogeneous structure
incorporating the observational error $\sigma^{2}_{y_{i}}$, that is,
$\varepsilon_{i}\sim N(0,\sigma^{2}_{y_{i}}+\sigma^{2})$. Our proposed method
can be applied to these two models straightforwardly. We assign a Gaussian
process prior $f\sim\mathrm{GP}(0,\sigma^{2}(n\lambda)^{-1}K)$, and use leave-
one-out cross validation to choose the kernel function from the Matérn kernel,
squared exponential kernel, and Sobolev kernel. For both error structures, the
Sobolev kernel is selected. We optimize $\sigma^{2}$ and $\lambda$ jointly by
maximizing the marginal likelihood function, noting that for model (i) the
closed-form expression for $\sigma^{2}$ is available as discussed in Section
2.4. Then the plug-in posterior distribution of $f^{\prime}$ is analytically
tractable and would require no further tuning or posterior sampling, based on
which one can make inference on the rate of global sea-level rise.

For both error structures, Figure 2 presents the rate estimation using the
posterior mean function along with 90% pointwise and simultaneous credible
bands for uncertainty quantification. Model (a) reveals that the rate of GMSL
rise is not a constant over time, and instead has accelerated much from 1.28
mm/yr in 1880 AD to 1.82 mm/yr in 2010 AD. Model (b) shows a similar
acceleration pattern but with a narrower overall range and credible band,
suggesting reduced uncertainty by incorporating the observational variance
$\sigma_{y_{i}}^{2}$ in the model. Our findings of the nonlinear, increasing
rate of global sea-level rise have also been observed in earlier literature
such as Cahill et al., (2015) based on integrated Gaussian processes, which
necessitate complex parameter tuning and computationally intensive posterior
sampling. Our uncertainty quantification indicates that the rate 1.7 mm/yr in
a least-squares line roughly falls into the 90% simultaneous credible
intervals from approximately 1955 AD to 1985 AD. Overall, our analysis
captures the acceleration in the rate of global sea rise that would otherwise
be missed by least-squares fitting, indicates the least-squares slope 1.7mm/yr
may be only representative for a limited time period, and provides a time-
varying description of global sea-level rise for the past 130 years.

(a) $\varepsilon_{i}\sim N(0,\sigma^{2})$.

(b) $\varepsilon_{i}\sim N(0,\sigma^{2}_{y_{i}}+\sigma^{2})$.

Figure 2: Rate of global sea-level rise calculated ($f^{\prime}$) under each
error structure. Shading denotes 90% pointwise (dark) and simultaneous
credible bands (light) for the rate process. The reference rate 1.7 mm/yr is
the least-squares slope in a linear model assuming a constant rate of sea-
level rise.

## Acknowledgements

We would like to thank William Yoo for providing R code to implement the
B-spline prior in the simulation section.

##  Appendix A Proofs in Section 2

###### Proof of Theorem 2.

Consider the centered Gaussian process:
$(f^{(k)}-\hat{f}^{(k)}_{n})|\mathbb{D}_{n}\sim{\rm GP}(0,\tilde{V}_{n}^{k})$.
According to Borell-TIS inequality (cf. Proposition A.2.1 in van der Vaart and
Wellner, (1996)), we have for any $j\in\mathbb{N}$,

$\Pi_{n,k}\left(\left|\|f^{(k)}-\hat{f}^{(k)}_{n}\|_{\infty}-\mathbb{E}\|f^{(k)}-\hat{f}^{(k)}_{n}\|_{\infty}\right|>j\epsilon_{n}\,\big{|}\,\mathbb{D}_{n}\right)\leq
2\exp\left(-j^{2}\epsilon_{n}^{2}/2\|\tilde{V}^{k}_{n}\|_{\infty}\right).$
(23)

Note that

$\displaystyle\Pi_{n,k}\left(\left|\|f^{(k)}-\hat{f}^{(k)}_{n}\|_{\infty}-\mathbb{E}\|f^{(k)}-\hat{f}^{(k)}_{n}\|_{\infty}\right|>j\epsilon_{n}\,\big{|}\,\mathbb{D}_{n}\right)$
(24)
$\displaystyle\qquad\qquad\geq\Pi_{n,k}\left(\|f^{(k)}-\hat{f}^{(k)}_{n}\|_{\infty}-\mathbb{E}\|f^{(k)}-\hat{f}^{(k)}_{n}\|_{\infty}>j\epsilon_{n}\,\big{|}\,\mathbb{D}_{n}\right).$
(25)

Thus, by Lemma 7, it holds with $\mathbb{P}_{0}^{(n)}$-probability at least
$1-n^{-10}$ that

$\|\tilde{V}^{k}_{n}\|_{\infty}\leq\frac{2\sigma^{2}\tilde{\kappa}_{k}^{2}}{n},$
(26)

leading to

$\Pi_{n,k}\left(\|f^{(k)}-\hat{f}^{(k)}_{n}\|_{\infty}-\mathbb{E}\|f^{(k)}-\hat{f}^{(k)}_{n}\|_{\infty}>j\epsilon_{n}\,\big{|}\,\mathbb{D}_{n}\right)\leq
2\exp\left(-n\epsilon_{n}^{2}j^{2}/4\sigma^{2}\tilde{\kappa}_{k}^{2}\right).$
(27)

By Lemma 12, there exists $C_{1}>0$ such that with probability at least
$1-n^{-10}$,

$\mathbb{E}\|f^{(k)}-\hat{f}^{(k)}_{n}\|_{\infty}\leq
C_{1}\tilde{\kappa}_{k}\sqrt{\frac{\log n}{n}}.$ (28)

We next consider the two cases when $p=\infty$ and $p=2$ separately. If
$p=\infty$, by Assumption (D), with $\mathbb{P}_{0}^{(n)}$-probability tending
to 1 we have $\|\hat{f}_{n}^{(k)}-f_{0}^{(k)}\|_{\infty}\lesssim\epsilon_{n}$
and thus

$\|f^{(k)}-\hat{f}_{n}^{(k)}\|_{\infty}\geq\|f^{(k)}-f_{0}^{(k)}\|_{\infty}-\|\hat{f}_{n}^{(k)}-f_{0}^{(k)}\|_{\infty}\geq\|f^{(k)}-f_{0}^{(k)}\|_{\infty}-C\epsilon_{n}$
(29)

for some $C>0$, which implies

$\Pi_{n,k}\left(\|f^{(k)}-f_{0}^{(k)}\|_{\infty}>C\epsilon_{n}+C_{1}\tilde{\kappa}_{k}\sqrt{\frac{\log
n}{n}}+j\epsilon_{n}\,\big{|}\,\mathbb{D}_{n}\right)\leq
2\exp\left(-n\epsilon_{n}^{2}j^{2}/4\sigma^{2}\tilde{\kappa}_{k}^{2}\right).$
(30)

According to Assumption (C), $\tilde{\kappa}_{k}^{2}=O(n\epsilon_{n}^{2}/\log
n)$ and $\tilde{\kappa}_{k}\sqrt{\frac{\log n}{n}}\leq C_{2}\epsilon_{n}$ for
some $C_{2}>0$, there exists a $J\in\mathbb{N}$ such that $C+C_{1}C_{2}\leq
J$. Thus, for any $j\geq J$, it holds with $\mathbb{P}_{0}^{(n)}$-probability
tending to 1 that

$\Pi_{n,k}\left(\|f^{(k)}-f_{0}^{(k)}\|_{\infty}>2j\epsilon_{n}\,\big{|}\,\mathbb{D}_{n}\right)\leq
2\exp\left(-n\epsilon_{n}^{2}j^{2}/4\sigma^{2}\tilde{\kappa}_{k}^{2}\right),$
(31)

which implies that

$\Pi_{n,k}\left(\|f^{(k)}-f_{0}^{(k)}\|_{\infty}>j\epsilon_{n}\,\big{|}\,\mathbb{D}_{n}\right)=\exp\left(-Kn\epsilon_{n}^{2}j^{2}/\tilde{\kappa}_{k}^{2}\right)$
(32)

for some $K>0$. Let $A_{n}$ denote the event that the preceding display holds,
which satisfies $\mathbb{P}_{0}^{(n)}(A_{n}^{c})\rightarrow 0$. Hence, for any
$M_{n}\rightarrow\infty$,

$\displaystyle\mathbb{P}_{0}^{(n)}\Pi_{n,k}\left(\|f^{(k)}-f_{0}^{(k)}\|_{\infty}>M_{n}\epsilon_{n}\,\big{|}\,\mathbb{D}_{n}\right)$
(33)
$\displaystyle\qquad\qquad\leq\mathbb{P}_{0}^{(n)}\Pi_{n,k}\left(\|f^{(k)}-f_{0}^{(k)}\|_{\infty}>M_{n}\epsilon_{n}\,\big{|}\,\mathbb{D}_{n}\right)\mathbbm{1}(A_{n})+\mathbb{P}_{0}^{(n)}(A_{n}^{c})\rightarrow
0.$ (34)

Therefore, $\epsilon_{n}$ is a contraction rate under the $L_{\infty}$ norm.

Now we consider the case of $p=2$. Since
$\|\hat{f}_{n}^{(k)}-f_{0}^{(k)}\|_{2}\lesssim\epsilon_{n}$ with
$\mathbb{P}_{0}^{(n)}$-probability tending to 1, we have

$\|f^{(k)}-\hat{f}_{n}^{(k)}\|_{\infty}\geq\|f^{(k)}-\hat{f}_{n}^{(k)}\|_{2}\geq\|f^{(k)}-f_{0}^{(k)}\|_{2}-\|\hat{f}_{n}^{(k)}-f_{0}^{(k)}\|_{2}\geq\|f^{(k)}-f_{0}^{(k)}\|_{2}-C\epsilon_{n}.$
(35)

Comparing the preceding display with (29) and following the same arguments, we
have with $\mathbb{P}_{0}^{(n)}$-probability tending to 1 that

$\Pi_{n,k}\left(\|f^{(k)}-f_{0}^{(k)}\|_{2}>j\epsilon_{n}\,\big{|}\,\mathbb{D}_{n}\right)\leq\exp\left(-Kn\epsilon_{n}^{2}j^{2}/\tilde{\kappa}_{k}^{2}\right).$
(36)

Then a similar argument as in the case of $p=\infty$ yields that
$\epsilon_{n}$ is also a posterior contraction rate. This completes the proof.
∎

###### Proof of Theorem 3.

Let $u_{1}\geq u_{2}\geq\ldots\geq u_{n}$ denote the eigenvalues of $K(X,X)$.
Since $K(X,X)$ is non-negative definite, we have $u_{i}\geq 0$ for $1\leq
i\leq n$. Note that $\sup_{x\in\mathcal{X}}K(x,x)<\infty$ as $K$ is a
continuous bivariate function on a compact support
$\mathcal{X}\times\mathcal{X}$. Then we have
$\sum_{i=1}^{n}u_{i}=\mathrm{tr}(K(X,X))\leq
n\cdot\sup_{x\in\mathcal{X}}K(x,x)\leq n\|K\|_{\infty}$. Let
$\bm{f}_{0}=(f_{0}(X_{1}),\ldots,f_{0}(X_{n}))^{T}$. The MMLE
$\hat{\sigma}_{n}^{2}$ is a quadratic form in $Y$. In view of the well-known
formula for the expectation of quadratic forms (cf. Theorem 11.19 in Schott,
(2016)), we obtain

$\mathbb{E}(\hat{\sigma}^{2}_{n}|X)=\lambda\sigma_{0}^{2}\mathrm{tr}([K(X,X)+n\lambda
I_{n}]^{-1})+\lambda\bm{f}_{0}^{T}[K(X,X)+n\lambda I_{n}]^{-1}\bm{f}_{0}.$
(37)

Therefore,

$\displaystyle|\mathbb{E}(\hat{\sigma}^{2}_{n}|X)-\sigma_{0}^{2}|$
$\displaystyle\leq\left|\lambda\sigma_{0}^{2}\mathrm{tr}([K(X,X)+n\lambda
I_{n}]^{-1})-\sigma_{0}^{2}\right|+\lambda\bm{f}^{T}_{0}[K(X,X)+n\lambda
I_{n}]^{-1}\bm{f}_{0}$ (38)
$\displaystyle\leq\left|n^{-1}\sigma_{0}^{2}\mathrm{tr}([(n\lambda)^{-1}K(X,X)+I_{n}]^{-1})-\sigma_{0}^{2}\right|+\lambda\bm{f}^{T}_{0}[K(X,X)+n\lambda
I_{n}]^{-1}\bm{f}_{0}.$ (39)

It follows that the first term is bounded by

$\displaystyle\left|n^{-1}\sigma_{0}^{2}\mathrm{tr}([(n\lambda)^{-1}K(X,X)+I_{n}]^{-1})-\sigma_{0}^{2}\right|$
(40)
$\displaystyle\qquad\qquad=n^{-1}\sigma_{0}^{2}\mathrm{tr}(I_{n}-[(n\lambda)^{-1}K(X,X)+I_{n}]^{-1})=n^{-1}\sigma_{0}^{2}\sum_{i=1}^{n}\left(1-\frac{1}{u_{i}/n\lambda+1}\right)$
$\displaystyle\qquad\qquad=n^{-1}\sigma_{0}^{2}\sum_{i=1}^{n}\frac{u_{i}}{n\lambda+u_{i}}\leq
n^{-1}\sigma_{0}^{2}\sum_{i=1}^{n}\frac{u_{i}}{n\lambda}\leq\sigma_{0}^{2}\|K\|_{\infty}(n\lambda)^{-1}.$

Let $\lambda_{\max}(A)$ denote the largest eigenvalue of a matrix $A$. We have
$\lambda_{\max}([K(X,X)+n\lambda
I_{n}]^{-1})=(u_{n}+n\lambda)^{-1}\leq(n\lambda)^{-1}$. Hence,

$\lambda\bm{f}^{T}_{0}[K(X,X)+n\lambda
I_{n}]^{-1}\bm{f}_{0}\leq\lambda\cdot\lambda_{\max}([K(X,X)+n\lambda
I_{n}]^{-1})\|\bm{f}_{0}\|_{2}^{2}\leq\|f\|_{\infty}^{2}n^{-1}.$ (41)

Combining (40) and (41) gives

$|\mathbb{E}(\hat{\sigma}^{2}_{n}|X)-\sigma_{0}^{2})|\leq\sigma_{0}^{2}\|K\|_{\infty}(n\lambda)^{-1}+\|f_{0}\|_{\infty}^{2}n^{-1}.$
(42)

We now bound the variance of $\hat{\sigma}^{2}_{n}$. Using the variance
formula for quadratic forms (cf. Theorem 11.23 in Schott, (2016)), we have

$\displaystyle\mathrm{Var}(\hat{\sigma}^{2}_{n}|X)$
$\displaystyle=2\lambda^{2}\sigma_{0}^{4}\mathrm{tr}([K(X,X)+n\lambda
I_{n}]^{-2})+4\lambda^{2}\sigma_{0}^{2}\bm{f}^{T}_{0}[K(X,X)+n\lambda
I_{n}]^{-2}\bm{f}_{0}$ (43) $\displaystyle\leq 2\lambda^{2}\sigma_{0}^{4}\cdot
n(n\lambda)^{-2}+4\lambda^{2}\sigma_{0}^{2}\cdot\lambda_{\max}([K(X,X)+n\lambda
I_{n}]^{-2})\|\bm{f}_{0}\|_{2}^{2}$ (44) $\displaystyle\leq
2\sigma_{0}^{4}n^{-1}+4\sigma_{0}^{2}\|f_{0}\|_{\infty}^{2}n^{-1}.$ (45)

Therefore, for any $X\in\mathcal{X}^{n}$,

$\displaystyle\mathbb{E}[(\hat{\sigma}^{2}_{n}-\sigma_{0}^{2})^{2}|X]$
$\displaystyle=\mathrm{Var}(\hat{\sigma}^{2}_{n}|X)+[\mathbb{E}(\hat{\sigma}^{2}_{n}|X)-\sigma_{0}^{2})]^{2}$
(46)
$\displaystyle\leq\sigma_{0}^{2}\|K\|_{\infty}(n\lambda)^{-1}+\|f_{0}\|_{\infty}^{2}n^{-1}+2\sigma_{0}^{4}n^{-1}+4\sigma_{0}^{2}\|f_{0}\|_{\infty}^{2}n^{-1},$
(47)

which implies that

$\mathbb{E}[(\hat{\sigma}^{2}_{n}-\sigma_{0}^{2})^{2}]=\mathbb{E}[\mathbb{E}[(\hat{\sigma}^{2}_{n}-\sigma_{0}^{2})^{2}|X]]=o(1).$
(48)

It thus follows that $\hat{\sigma}^{2}_{n}$ converges to $\sigma_{0}^{2}$ in
$\mathbb{P}_{0}^{(n)}$-probability by applying Chebyshev’s inequality.

Now we prove Theorem 2 under the empirical Bayes scheme. Consider a shrinking
neighborhood $\mathcal{B}_{n}=(\sigma_{0}^{2}-r_{n},\sigma_{0}^{2}+r_{n})$ of
$\sigma_{0}^{2}$ with $r_{n}=o(1)$, which satisfies that
$\mathbb{P}_{0}^{(n)}(\hat{\sigma}^{2}_{n}\in\mathcal{B}_{n})\rightarrow 1$
according to the arguments above. Conditional on $\mathcal{B}_{n}$, Lemma 7
gives

$\|\tilde{V}^{k}_{n}\|_{\infty}\leq\frac{2(\sigma_{0}^{2}+o(1))\tilde{\kappa}_{k}^{2}}{n}.$
(49)

Then, all the established inequalities in the proof of Theorem 2 hold
uniformly over $\sigma^{2}\in\mathcal{B}_{n}$. It follows that

$\sup_{\sigma^{2}\in\mathcal{B}_{n}}\Pi_{n,k}\left(f:\|f^{(k)}-f_{0}^{(k)}\|_{p}\geq
M_{n}\epsilon_{n}\big{|}\,\mathbb{D}_{n},\sigma^{2}\right)\rightarrow 0$ (50)

in $\mathbb{P}_{0}^{(n)}$-probability for $p=2,\infty$, which directly implies
that

$\Pi_{n,k,\text{EB}}\left(f:\|f^{(k)}-f_{0}^{(k)}\|_{p}\geq
M_{n}\epsilon_{n}\big{|}\,\mathbb{D}_{n}\right)\rightarrow 0$ (51)

in $\mathbb{P}_{0}^{(n)}$-probability for $p=2,\infty$. This completes the
proof. ∎

###### Proof of Theorem 4.

Note that the Fourier basis $\\{\psi_{i}\\}_{i=1}^{\infty}$ satisfies
Assumption (A) with $C_{k,\psi}=\sqrt{2}(2\pi)^{k}$ and Assumption (B) with
$L_{k,\psi}=\sqrt{2}(2\pi)^{k+1}$. According to Lemma 8, we have
$\tilde{\kappa}_{\gamma,m}^{2}\asymp(-\log\lambda)^{2m+1}$ for any
$m\in\mathbb{N}_{0}$. Since $\lambda\asymp 1/n$, we have
$\tilde{\kappa}_{\gamma}^{2}\asymp(-\log\lambda)\asymp\log n=o(\sqrt{n/\log
n})$. Besides, $\tilde{\kappa}_{\gamma,k}^{2}\asymp(\log
n)^{2k+1}=O(n\epsilon_{n}^{2}/\log n)$. It also follows that
$\hat{\kappa}_{\gamma,k+1}^{2}\asymp(-\log\lambda)^{2k+3}\asymp(\log
n)^{2k+3}=O(n)$ for any $k\in\mathbb{N}_{0}$. This verifies Assumption (C).
Finally, Assumption (D) is given by Lemma 10, which shows a convergence rate
of $\hat{f}_{n}^{(k)}$ under the $L_{2}$ norm is $\epsilon_{n}=(\log
n)^{k+1}/\sqrt{n}$ when $\gamma_{0}>\gamma>1/2$. Invoking Theorem 2,
$\epsilon_{n}$ is a contraction rate of the posterior distribution
$\Pi_{n,k}(\cdot\mid\mathbb{D}_{n})$. This completes the proof. ∎

###### Proof of Theorem 5.

The arguments are similar to those used to prove Theorem 4, and we only note
the key differences below. We first verify Assumption (C). According to Lemma
9, we have
$\tilde{\kappa}_{\alpha,m}^{2}\asymp\lambda^{-\frac{2m+1}{2\alpha}}$ for any
$m<\alpha-1/2$ and $m\in\mathbb{N}$. Since $\lambda\asymp(\log
n/n)^{\frac{2\alpha}{2\alpha+1}}$, we have
$\tilde{\kappa}_{\alpha}^{2}\asymp\lambda^{-\frac{1}{2\alpha}}\asymp({n}/{\log
n})^{\frac{1}{2\alpha+1}}=o(\sqrt{n/\log n})$. Besides,
$\tilde{\kappa}_{\alpha,k}^{2}\asymp(n/\log
n)^{\frac{2k+1}{2\alpha+1}}=O(n\epsilon_{n}^{2}/\log n)$. It also follows that
$\hat{\kappa}_{\alpha,k+1}^{2}\asymp\lambda^{-\frac{2k+3}{2\alpha}}\asymp\left(n/\log
n\right)^{\frac{2k+3}{2\alpha+1}}=O(n)$ for any $k<\alpha-3/2$ and
$k\in\mathbb{N}_{0}$. Assumption (D) follows from Lemma 11. This completes the
proof. ∎

##  Appendix B Preliminaries and Proofs in Section 3

###  B.1 Preliminaries for reproducing kernel Hilbert space

Before presenting the proofs, we review preliminaries for reproducing kernel
Hilbert space (RKHS). These background materials can be found in many existing
works such as Rasmussen and Williams, (2006); Berlinet and Thomas-Agnan,
(2011), and are included here for completeness and introducing notation.

Let $\mathbb{H}$ be the reproducing kernel Hilbert space induced by the kernel
$K$. For any $f=\sum_{i=1}^{\infty}f_{i}\phi_{i}\in\mathbb{H}$ and
$g=\sum_{i=1}^{\infty}g_{i}\phi_{i}\in\mathbb{H}$, the inner product is
$\left<f,g\right>_{\mathbb{H}}=\sum_{i=1}^{\infty}f_{i}g_{i}/\mu_{i}$. For any
$f\in L^{2}_{p_{X}}(\mathcal{X})$, we define an integral operator
$L_{K}:L^{2}_{p_{X}}(\mathcal{X})\rightarrow\mathbb{H}$:

$L_{K}(f)(x)=\int_{\mathcal{X}}K(x,x^{\prime})f(x^{\prime})d\mathbb{P}_{X}(x^{\prime}),\quad
x\in\mathcal{X}.$ (52)

The compact, positive definite and self-adjoint operator $L_{K}$ and the
kernel $K$ can be uniquely determined by each other. $L_{K}$ can also be
defined on $\mathbb{H}$ or $L^{2}_{p_{X}}(\mathcal{X})$.

We introduce the sample analog
$L_{K,X}:L^{2}_{p_{X}}(\mathcal{X})\rightarrow\mathbb{H}$ by

$L_{K,X}(f)=\frac{1}{n}\sum_{i=1}^{n}f(X_{i})K_{X_{i}},$ (53)

where $K_{X_{i}}(\cdot)=K(X_{i},\cdot)$. Note that the operator $L_{K,X}$ and
the matrix $K(X,X)/n$ share the same eigenvalues. We remark that $L_{K}$ and
$L_{K,X}$ can also be defined on subspaces of $L^{2}_{p_{X}}(\mathcal{X})$
such as $C(\mathcal{X})$, $C^{k}(\mathcal{X})$, and $H^{\alpha}(\mathcal{X})$.

Let $I$ be the identity operator. We approximate $f_{0}$ by a function in
$\mathbb{H}$ defined as follows:

$f_{\lambda}=(L_{K}+\lambda I)^{-1}L_{K}f_{0},$ (54)

which minimizes $\|f-f_{0}\|_{2}^{2}+\lambda\|f\|^{2}_{\mathbb{H}}$ subject to
$f\in\mathbb{H}.$

Let $\tilde{\mathbb{H}}$ be the RKHS induced by the equivalent kernel
$\tilde{K}$, which is equivalent to $\mathbb{H}$ as a function space, but with
a different inner product

$\left<f,g\right>_{\tilde{\mathbb{H}}}=\left<f,g\right>_{2}+\lambda\left<f,g\right>_{\mathbb{H}}.$
(55)

We call the corresponding norm $\|\cdot\|_{\tilde{\mathbb{H}}}$ the equivalent
RKHS norm. The equivalent RKHS norm upper bounds the $L_{\infty}$ norm by
noting that

$f(x)=\left<f,\tilde{K}_{x}\right>_{\tilde{\mathbb{H}}}\leq\|f\|_{\tilde{\mathbb{H}}}\|K_{x}\|_{\tilde{\mathbb{H}}}=\sqrt{\tilde{K}(x,x)}\|f\|_{\tilde{\mathbb{H}}},$
(56)

which follows from Cauchy-Schwarz inequality. Taking the supremum on both
sides, we obtain that for any $f\in\tilde{\mathbb{H}}$, it holds

$\|f\|_{\infty}\leq\tilde{\kappa}\|f\|_{\tilde{\mathbb{H}}}.$ (57)

The equivalent kernel $\tilde{K}$ provides an alternative interpretation of
the proximate function $f_{\lambda}$. Let the integral operator
$L_{\tilde{K}}$ be the counterpart of $L_{K}$ induced by $\tilde{K}$, then we
have $f_{\lambda}=L_{\tilde{K}}f_{0}$.

###  B.2 Proofs in Section 3

###### Proof of Lemma 7.

Let $K_{0k,x^{\prime}}(\cdot)=\partial_{x^{\prime}}^{k}K(\cdot,x^{\prime})$
and

$\hat{K}_{0k,x^{\prime}}(\cdot)=K(\cdot,X)[K(X,X)+n\lambda
I_{n}]^{-1}K_{0k}(X,x^{\prime}).$

It is easy to see that $\hat{K}_{0k,x^{\prime}}$ is the solution to a noise-
free KRR with observations $K_{0k,x^{\prime}}$. Moreover, we have

$\begin{split}|\sigma^{-2}n\lambda\tilde{V}^{k}_{n}(x^{\prime})|&=\left|\partial^{k}_{x}(\hat{K}_{0k,x^{\prime}}-K_{0k,x^{\prime}})(x)|_{x=x^{\prime}}\right|\\\
&\leq\|\partial^{k}_{x}(\hat{K}_{0k,x^{\prime}}-K_{0k,x^{\prime}})\|_{\infty}\leq\tilde{\kappa}_{k}\|\hat{K}_{0k,x^{\prime}}-K_{0k,x^{\prime}}\|_{\tilde{\mathbb{H}}}.\end{split}$
(58)

According to Lemma 15, it holds with $\mathbb{P}_{0}^{(n)}$-probability at
least $1-n^{-10}$ that

$\|\hat{K}_{0k,x^{\prime}}-K_{0k,x^{\prime}}\|_{\tilde{\mathbb{H}}}\leq
2\|L_{\tilde{K}}K_{0k,x^{\prime}}-K_{0k,x^{\prime}}\|_{\tilde{\mathbb{H}}}\leq
2\lambda\tilde{\kappa}_{k}.$ (59)

∎

###### Proof of Lemma 8.

According to the proof of Lemma 11 in Liu and Li, (2022), we have

$\tilde{\kappa}_{\gamma,m}^{2}\asymp\sum_{i=1}^{\infty}\frac{i^{2m}\mu_{i}}{\lambda+\mu_{i}}=\hat{\kappa}_{\gamma,m}^{2}\asymp\sum_{i=1}^{\infty}\frac{e^{-2\gamma
i}i^{2m}}{\lambda+e^{-2\gamma i}}\asymp\int_{0}^{\infty}\frac{x^{2m}}{\lambda
e^{2\gamma x}+1}dx.$ (60)

Letting $e^{-\eta}=\lambda\in(0,1)$ where $\eta>0$, it follows that

$\tilde{\kappa}_{\gamma,m}^{2}\asymp(2\gamma)^{-2m-1}\int_{0}^{\infty}\frac{t^{2m}}{\lambda
e^{t}+1}dt=(2\gamma)^{-2m-1}(2m)!F_{2m}(\eta),$ (61)

where $F_{2m}$ is the Fermi-Dirac integral of order $2m$ (McDougall and
Stoner,, 1938).

According to Equation (7.3) in McDougall and Stoner, (1938), as
$\eta\rightarrow\infty$ (or $\lambda\rightarrow 0$), it holds that

$F_{2m}(\eta)\rightarrow\frac{1}{2m+1}\eta^{2m+1}=\frac{1}{2m+1}(-\log\lambda)^{2m+1}.$
(62)

Therefore,

$\tilde{\kappa}_{\gamma,m}^{2}\asymp\hat{\kappa}_{\gamma,m}^{2}\asymp(-\log\lambda)^{2m+1}.$
(63)

The differentiability of $\tilde{K}_{\gamma}$ directly follows from the
boundedness of $\tilde{\kappa}^{2}_{\gamma,m}$ for sufficiently small
$\lambda$. ∎

###### Proof of Lemma 9.

According to the proof of Lemma 11 in Liu and Li, (2022), it holds that

$\tilde{\kappa}_{\gamma,m}^{2}\asymp\sum_{i=1}^{\infty}\frac{i^{2m}\mu_{i}}{\lambda+\mu_{i}}=\hat{\kappa}_{\gamma,m}^{2}\asymp\sum_{i=1}^{\infty}\frac{i^{-2\alpha}i^{2m}}{\lambda+i^{-2\alpha}}\asymp\lambda^{-\frac{2m+1}{2\alpha}}.$
(64)

∎

###### Proof of Lemma 10.

According to the proof of Theorem 2 in Liu and Li, (2022) and substituting
$\delta=n^{-10}$, it holds with $\mathbb{P}_{0}^{(n)}$-probability at least
$1-n^{-10}$ that

$\|\hat{f}_{n}-f_{\lambda}\|_{\tilde{\mathbb{H}}}\leq\frac{\tilde{\kappa}^{-1}C(n,\tilde{\kappa})}{1-C(n,\tilde{\kappa})}\|f_{\lambda}-f_{0}\|_{\infty}+\frac{1}{1-C(n,\tilde{\kappa})}\frac{C_{1}\tilde{\kappa}\sigma\sqrt{10\log(3n)}}{\sqrt{n}},$
(65)

where $C_{1}>0$ does not depend on $K$ or $n$ and
$C(n,\tilde{\kappa})=\frac{\tilde{\kappa}^{2}\sqrt{10\log(3n)}}{\sqrt{n}}\left(4+\frac{4\tilde{\kappa}\sqrt{10\log(3n)}}{3\sqrt{n}}\right)$.
Applying Lemma 13 to the preceding inequality yields that

$\displaystyle\|\hat{f}_{n}^{(k)}-f_{\lambda}^{(k)}\|_{2}$
$\displaystyle\leq\tilde{\kappa}_{\gamma}^{-1}\tilde{\kappa}_{\gamma,k}\|\hat{f}_{n}-f_{\lambda}\|_{\tilde{\mathbb{H}}_{\gamma}}$
(66)
$\displaystyle\leq\frac{\tilde{\kappa}_{\gamma,k}\tilde{\kappa}_{\gamma}^{-2}C(n,\tilde{\kappa}_{\gamma})}{1-C(n,\tilde{\kappa}_{\gamma})}\|f_{\lambda}-f_{0}\|_{\infty}+\frac{1}{1-C(n,\tilde{\kappa}_{\gamma})}\frac{C_{1}\tilde{\kappa}_{\gamma,k}\sigma\sqrt{10\log(3n)}}{\sqrt{n}}.$
(67)

Note that
$\tilde{\kappa}_{\gamma}^{-2}C(n,\tilde{\kappa}_{\gamma})\asymp\sqrt{\log
n/n}$. By choosing $\lambda$ such that
$\tilde{\kappa}_{\gamma}^{2}=o(\sqrt{n/\log n})$,
$\tilde{\kappa}_{\gamma,k}=o(\sqrt{n/\log n})$ and
$C(n,\tilde{\kappa}_{\gamma})\leq 1/2$, we arrive at

$\|\hat{f}_{n}^{(k)}-f_{\lambda}^{(k)}\|_{2}\lesssim
2\|f_{\lambda}-f_{0}\|_{\infty}+\frac{2C_{1}\tilde{\kappa}_{\gamma,k}\sigma\sqrt{10\log(3n)}}{\sqrt{n}}.$
(68)

According to Lemma 14, we obtain

$\displaystyle\|\hat{f}_{n}^{(k)}-f_{0}^{(k)}\|_{2}$
$\displaystyle\lesssim\|f_{\lambda}^{(k)}-f_{0}^{(k)}\|_{2}+2\|f_{\lambda}-f_{0}\|_{\infty}+\frac{2C_{1}\tilde{\kappa}_{\gamma,k}\sigma\sqrt{10\log(3n)}}{\sqrt{n}}$
(69)
$\displaystyle\leq\|f_{\lambda}^{(k)}-f_{0}^{(k)}\|_{\infty}+2\|f_{\lambda}-f_{0}\|_{\infty}+\frac{2C_{1}\tilde{\kappa}_{\gamma,k}\sigma\sqrt{10\log(3n)}}{\sqrt{n}}$
(70)
$\displaystyle\lesssim\lambda^{\frac{1}{2}}+\lambda^{\frac{1}{2}}+(-\log\lambda)^{m+\frac{1}{2}}\sqrt{\frac{\log
n}{n}}.$ (71)

We take $\lambda\asymp\log n/n$, which satisfies
$\tilde{\kappa}_{\gamma}^{2}\asymp\log(n/\log n)=o(\sqrt{n/\log n})$ and
$\tilde{\kappa}_{\gamma,k}\asymp(\log(n/\log
n))^{k+\frac{1}{2}}=o(\sqrt{n/\log n})$. The proof is completed by
substituting $\lambda$. ∎

###### Proof of Lemma 11.

See the proof of Theorem 14 in Liu and Li, (2022).

∎

##  Appendix C Auxiliary technical results and proofs

######  Lemma 12.

Under the conditions of Theorem 2, let $\tilde{\rho}$ be the intrinsic psuedo-
metric of the centered Gaussian process (73). Then, it holds with
$\mathbb{P}_{0}^{(n)}$-probability at least $1-n^{-10}$ that

$\mathbb{E}\|f^{(k)}-\hat{f}^{(k)}_{n}\|_{\infty}\lesssim\tilde{\kappa}_{k}\sqrt{\frac{\log
n}{n}}.$ (72)

###### Proof of Lemma 12.

Define a centered GP posterior $\tilde{Z}^{k}$ as

$\tilde{Z}^{k}=(f^{(k)}-\hat{f}^{(k)}_{n})|\mathbb{D}_{n}\sim{\rm
GP}(0,\tilde{V}_{n}^{k}),$ (73)

with psuedo-metric $\tilde{\rho}(x,x^{\prime})=\sqrt{{\rm
Var}(\tilde{Z}^{k}(x)-\tilde{Z}^{k}(x^{\prime}))}$ for $x,x^{\prime}\in[0,1]$.
By Dudley’s entropy integral Theorem,

$\mathbb{E}\|f^{(k)}-\hat{f}^{(k)}_{n}\|_{\infty}\leq\int_{0}^{\tilde{\rho}([0,1])}\sqrt{\log
N(\epsilon,[0,1],\tilde{\rho})}d\epsilon.$ (74)

Recall that in (58) the posterior covariance can be expressed as the bias of a
noise-free KRR estimator:

$\sigma^{-2}n\lambda\tilde{V}^{k}_{n}(x,x^{\prime})=\partial^{k}_{x}(K_{0k,x^{\prime}}-\hat{K}_{0k,x^{\prime}})(x).$
(75)

Hence, we have

$\displaystyle\tilde{\rho}(x,x^{\prime})^{2}$ $\displaystyle={\rm
Var}(\tilde{Z}^{k}(x),\tilde{Z}^{k}(x^{\prime}))=\tilde{V}^{k}_{n}(x,x)+\tilde{V}^{k}_{n}(x^{\prime},x^{\prime})-2\tilde{V}^{k}_{n}(x,x^{\prime})$
(76)
$\displaystyle=\sigma^{2}(n\lambda)^{-1}\\{\partial^{k}_{x}K_{0k,x}(x)-\partial^{k}_{x}\hat{K}_{0k,x}(x)+\partial^{k}_{x^{\prime}}K_{0k,x^{\prime}}(x^{\prime})-\partial^{k}_{x^{\prime}}\hat{K}_{0k,x^{\prime}}(x^{\prime})$
(77) $\displaystyle\qquad\qquad\quad\
-2\partial^{k}_{x}K_{0k,x^{\prime}}(x)+2\partial^{k}_{x}\hat{K}_{0k,x^{\prime}}(x)\\}$
(78)
$\displaystyle=\sigma^{2}(n\lambda)^{-1}\\{\partial^{k}_{x}(K_{0k,x}-K_{0j,x^{\prime}})(x)-\partial^{k}_{x^{\prime}}(K_{0k,x}-K_{0k,x^{\prime}})(x^{\prime})$
(79) $\displaystyle\qquad\qquad\quad\
-\partial^{k}_{x}(\hat{K}_{0k,x}-\hat{K}_{0k,x^{\prime}})(x)+\partial^{k}_{x^{\prime}}(\hat{K}_{0k,x}-\hat{K}_{0k,x^{\prime}})(x^{\prime})\\}.$
(80)

Let $g_{0k}=K_{0k,x}-K_{0k,x^{\prime}}$ and
$\hat{g}_{0k}=\hat{K}_{0k,x}-\hat{K}_{0k,x^{\prime}}$. Then the preceding
display implies

$\sigma^{-2}n\lambda\tilde{\rho}(x,x^{\prime})^{2}=\partial^{k}_{x}(g_{0k}-\hat{g}_{0k})(x)-\partial^{k}_{x^{\prime}}(g_{0k}-\hat{g}_{0k})(x^{\prime}).$
(81)

By Lemma 15, with $\mathbb{P}_{0}^{(n)}$-probability at least $1-n^{-10}$ it
holds that

$\begin{split}\sigma^{-2}n\lambda\tilde{\rho}(x,x^{\prime})^{2}&\leq
2\|\partial^{k}(\hat{g}_{0k}-g_{0k})\|_{\infty}\\\ &\leq
2\tilde{\kappa}_{k}\|\hat{g}_{0k}-g_{0k}\|_{\tilde{\mathbb{H}}}\\\ &\leq
4\tilde{\kappa}_{k}\|L_{\tilde{K}}g_{0k}-g_{0k}\|_{\tilde{\mathbb{H}}}.\end{split}$
(82)

Substituting
$g_{0k}=\sum_{i=1}^{\infty}\mu_{i}(\phi^{(k)}_{i}(x)-\phi^{(k)}_{i}(x^{\prime}))\phi_{i}$
yields

$\begin{split}\|L_{\tilde{K}}g_{0k}-g_{0k}\|_{\tilde{\mathbb{H}}}^{2}&=\sum_{i=1}^{\infty}\left(\frac{\lambda\mu_{i}(\phi^{(k)}_{i}(x)-\phi^{(k)}_{i}(x^{\prime}))}{\lambda+\mu_{i}}\right)^{2}\bigg{/}\frac{\mu_{i}}{\lambda+\mu_{i}}\\\
&\leq\lambda^{2}L_{k,\phi}^{2}|x-x^{\prime}|^{2}\sum_{i=1}^{\infty}\frac{i^{2k+2}\mu_{i}}{\lambda+\mu_{i}}\\\
&\leq
L_{k,\phi}^{2}\lambda^{2}\hat{\kappa}_{k+1}^{2}|x-x^{\prime}|^{2}.\end{split}$
(83)

It follows from (82) and (83) that

$\tilde{\rho}(x,x^{\prime})\leq
C\sqrt{\tilde{\kappa}_{k}\hat{\kappa}_{k+1}}n^{-1/2}|x-x^{\prime}|^{1/2}.$
(84)

By the inequality in Ghosal and van der Vaart, (2017, p. 529), we have

$N(\epsilon,[0,1],\tilde{\rho})\leq
N\left(\left(C^{-1}\epsilon\sqrt{n/\tilde{\kappa}_{k}\hat{\kappa}_{k+1}}\right)^{2},[0,1],|\cdot|\right)\lesssim\frac{1}{(\epsilon\sqrt{n/\tilde{\kappa}_{k}\hat{\kappa}_{k+1}})^{2}}.$
(85)

Since
$\tilde{\kappa}_{k}\hat{\kappa}_{k+1}\lesssim\hat{\kappa}_{k}\hat{\kappa}_{k+1}\leq\hat{\kappa}_{k+1}^{2}=O(n)$,
we obtain $\log N(\epsilon,[0,1],\tilde{\rho})\lesssim\log(1/\epsilon)$. On
the other hand, note that

$\displaystyle\|L_{\tilde{K}}g-g\|_{\tilde{\mathbb{H}}}^{2}$
$\displaystyle=\sum_{i=1}^{\infty}\left(\frac{\lambda\mu_{i}(\phi_{i}^{(k)}(x)-\phi_{i}^{(k)}(x^{\prime}))}{\lambda+\mu_{i}}\right)^{2}\bigg{/}\frac{\mu_{i}}{\lambda+\mu_{i}}$
(86)
$\displaystyle\leq\sum_{i=1}^{\infty}\frac{2\lambda^{2}\mu_{i}(\phi_{i}^{(k)}(x)^{2}+\phi_{i}^{(k)}(x^{\prime})^{2})}{\lambda+\mu_{i}}$
(87) $\displaystyle\leq 4\lambda^{2}\tilde{\kappa}_{k}^{2}.$ (88)

The preceding inequality combined with (82) gives that

$\tilde{\rho}(x,x^{\prime})^{2}\leq
4\sigma^{2}(n\lambda)^{-1}\tilde{\kappa}_{k}\|L_{\tilde{K}}g-g\|_{\tilde{\mathbb{H}}}\leq
8\sigma^{2}\tilde{\kappa}_{k}^{2}/n,$ (89)

which implies $\tilde{\rho}([0,1])\lesssim\tilde{\kappa}_{k}/\sqrt{n}$.
Therefore,

$\int_{0}^{\tilde{\rho}([0,1])}\sqrt{\log
N(\epsilon,[0,1],\tilde{\rho})}d\epsilon\lesssim\int_{0}^{\tilde{\kappa}_{k}/\sqrt{n}}\sqrt{\log(1/\epsilon)}d\epsilon\lesssim\tilde{\kappa}_{k}\sqrt{\frac{\log
n}{n}}.$ (90)

∎

######  Lemma 13.

Let the RKHS induced by $K_{\gamma}$ and $\tilde{K}_{\gamma}$ be
$\mathbb{H}_{\gamma}$ and $\tilde{\mathbb{H}}_{\gamma}$, respectively. If
$\gamma>1/2$, then $f\in C^{m}[0,1]$ for any $m\in\mathbb{N}_{0}$ and
$f\in\mathbb{H}_{\gamma}$. Moreover, these exists a constant $C>0$ that does
not depend on $\lambda$ such that $\|f^{(m)}\|_{2}\leq
C\tilde{\kappa}_{\gamma}^{-1}\tilde{\kappa}_{\gamma,m}\|f\|_{\tilde{\mathbb{H}}_{\gamma}}$
for any $0<\lambda<1$ and $f\in\mathbb{H}_{\gamma}$.

###### Proof of Lemma 13.

In view of Corollary 4.36 in Steinwart and Christmann, (2008),
$\tilde{K}_{\gamma}\in C^{2m}([0,1]\times[0,1])$. This implies that $f\in
C^{m}[0,1]$ for any $f\in\tilde{\mathbb{H}}_{\gamma}$, which is also true for
$f\in\mathbb{H}_{\gamma}$ since $\mathbb{H}_{\gamma}$ and
$\tilde{\mathbb{H}}_{\gamma}$ contain the same functions.

Now we prove the norm inequality. Let $f=\sum_{i=1}^{\infty}f_{i}\psi_{i}$
where $\\{\psi_{i}\\}_{i=1}^{\infty}$ is the Fourier basis. Then
$\|f^{(m)}\|_{2}^{2}\asymp\sum_{i=1}^{\infty}(f_{i}i^{m})^{2}$ for any
$m\in\mathbb{N}_{0}$. It is equivalent to show that

$\tilde{\kappa}_{\gamma}^{2}\cdot\sum_{i=1}^{\infty}f_{i}^{2}i^{2m}\leq
C\tilde{\kappa}_{\gamma,m}^{2}\cdot\sum_{i=1}^{\infty}f_{i}^{2}\frac{\lambda+\mu_{i}}{\mu_{i}},$
(91)

for some $C>0$. Hence, it suffices to show that for any $i\in\mathbb{N}$,

$\tilde{\kappa}_{\gamma}^{2}\cdot f_{i}^{2}i^{2m}\leq
C\tilde{\kappa}_{\gamma,m}^{2}\cdot f_{i}^{2}\frac{\lambda+\mu_{i}}{\mu_{i}}.$
(92)

In view of Lemma 9, we have
$\tilde{\kappa}_{\gamma,m}^{2}\asymp(-\log\lambda)^{2m+1}$ and
$\tilde{\kappa}_{\gamma}^{2}\asymp-\log\lambda$. Since $\mu_{i}\asymp
e^{-2\gamma i}$, the preceding inequality becomes

$i^{2m}\leq C(-\log\lambda)^{2m}(1+\lambda e^{2\gamma i}).$ (93)

The above sufficient condition in (93) trivially holds when
$i\leq-\log\lambda$ by taking $C=1$. When $i>-\log\lambda$, it follows that

$\frac{i^{2m}}{(-\log\lambda)^{2m}(1+\lambda e^{2\gamma
i})}\leq\frac{i^{2m}}{(-\log\lambda)^{2m}\lambda e^{2\gamma
i}}\leq\frac{i^{2m}}{(-\log\lambda)^{2m}e^{-i}e^{2\gamma i}}.$ (94)

This is bounded above by some constant $C$ since $\gamma>1/2$, leading to
(93). This completes the proof. ∎

######  Lemma 14.

Suppose $f_{0}\in A^{\gamma_{0}}[0,1]$ for some $\gamma_{0}>\gamma>0$. If
$K_{\gamma}$ is used in the GP prior, then for any $k\in\mathbb{N}_{0}$ it
holds

$\|f_{\lambda}^{(k)}-f_{0}^{(k)}\|_{\infty}\lesssim\sqrt{\lambda}.$ (95)

###### Proof of Lemma 14.

Let $f_{0}=\sum_{i=1}^{\infty}f_{i}\psi_{i}$. Then, for any
$k\in\mathbb{N}_{0}$,

$f_{\lambda}^{(k)}-f_{0}^{(k)}=-\sum_{i=1}^{\infty}\frac{\lambda}{\lambda+\mu_{i}}f_{i}\psi_{i}^{(k)}.$
(96)

Hence,

$\displaystyle\|f_{\lambda}^{(k)}-f_{0}^{(k)}\|_{\infty}$
$\displaystyle\leq\sum_{i=1}^{\infty}\frac{\lambda}{\lambda+\mu_{i}}|f_{i}|\cdot
i^{k}\lesssim\sqrt{\lambda}\sum_{i=1}^{\infty}\frac{\sqrt{\lambda\mu_{i}}}{\lambda+\mu_{i}}\frac{|f_{i}|i^{k}}{\sqrt{\mu_{i}}}$
(97) $\displaystyle\leq\sqrt{\lambda}\sum_{i=1}^{\infty}|f_{i}|i^{k}e^{\gamma
i}\lesssim\sqrt{\lambda}\sum_{i=1}^{\infty}|f_{i}|e^{\gamma_{0}i}\lesssim\sqrt{\lambda}.$
(98)

∎

######  Lemma 15.

Suppose the observations are noiseless. Under Assumption (A), by choosing
$\lambda$ such that $\tilde{\kappa}^{2}=o(\sqrt{n/\log n})$, it holds with
$\mathbb{P}_{0}^{(n)}$-probability at least $1-n^{-10}$ that

$\|\hat{f}_{n}-f_{0}\|_{\tilde{\mathbb{H}}}\leq
2\|f_{\lambda}-f_{0}\|_{\tilde{\mathbb{H}}}.$ (99)

###### Proof of Lemma 15.

A noise-free version of (65) by substituting $\sigma=0$ yields that with
$\mathbb{P}_{0}^{(n)}$-probability at least $1-n^{-10}$,

$\|\hat{f}_{n}-f_{\lambda}\|_{\tilde{\mathbb{H}}}\leq\frac{\tilde{\kappa}^{-1}C(n,\tilde{\kappa})}{1-C(n,\tilde{\kappa})}\|f_{\lambda}-f_{0}\|_{\infty}.$
(100)

Since $\tilde{\kappa}^{2}=o(\sqrt{n/\log n})$, for sufficiently large $n$ we
have $C(n,\tilde{\kappa})\leq 1/2$, and therefore the preceding inequality
becomes

$\|\hat{f}_{n}-f_{\lambda}\|_{\tilde{\mathbb{H}}}\leq\tilde{\kappa}^{-1}\|f_{\lambda}-f_{0}\|_{\infty}.$
(101)

Then we have

$\displaystyle\|\hat{f}_{n}-f_{0}\|_{\tilde{\mathbb{H}}}$
$\displaystyle\leq\|\hat{f}_{n}-f_{\lambda}\|_{\tilde{\mathbb{H}}}+\|f_{\lambda}-f_{0}\|_{\tilde{\mathbb{H}}}$
(102)
$\displaystyle\leq\tilde{\kappa}^{-1}\|f_{\lambda}-f_{0}\|_{\infty}+\|f_{\lambda}-f_{0}\|_{\tilde{\mathbb{H}}}$
(103) $\displaystyle\leq 2\|f_{\lambda}-f_{0}\|_{\tilde{\mathbb{H}}},$ (104)

where the last inequality follows from (57). ∎

##  Appendix D Additional simulation results

In this additional simulation, we compare our plug-in GP estimators with the
inverse method proposed in Holsclaw et al., (2013) using one of their
simulation designs, where the authors perceived the plug-in estimator as
suboptimal. We also compare the empirical Bayes strategy used in the main
paper with a fully Bayesian approach by placing an inverse Gamma and Gamma
priors on $\sigma^{2}$ and $\lambda$, respectively, which do not show
significant differences between the two treatments under the simulation
settings.

We consider the regression function $f_{0}(x)=x\sin(x)/10$ given as one
simulated example in Holsclaw et al., (2013). We generate $n=100$ and 500 data
points on a regular grid in $[0,10]$, and add noise $\varepsilon_{i}\sim
N(0,\sigma^{2})$ with $\sigma=0.1$, 0.2, and 0.3. The goal is to estimate
$f_{0}^{\prime}$. The implementations of the plug-in GP method and the
B-spline prior are the same as in the main paper. For the inverse method, we
follow the priors and general setup in Holsclaw et al., (2013), where a power-
exponential kernel with $\alpha=1.99$ is used.

Table 2 reports the RMSEs of $f_{0}^{\prime}$ calculated at 100 equally spaced
points in $[0,10]$ over 100 repetitions. We notice that in some repetitions
the inverse method produces a zero function as the estimate, and we exclude
these simulations and denote the results as “Inverse*”; this is probably due
to the complicated hyperparameter tuning without clear guides, which hampers
easy implementation and optimal performance. Moreover, the inverse method is
restricted to one particular derivative order, and its generalization to other
derivative orders is nontrivial. From Table 2, we can see that the plug-in
Gaussian process prior with Sobolev kernel achieves best performance among all
scenarios in general. The inverse method and B-splines behave relatively well
for $n=100$ and $n=500$, respectively. In summary, our GP method tends to
produces better estimation for larger $n$ and smaller $\sigma$, i.e., when the
signal-to-noise ratio is large, and we have found no scenario in the
considered settings where the inverse method has significantly smaller RMSE
relative to standard errors.

Table 2: RMSE of estimating $f_{0}^{\prime}$, averaged over 100 repetitions.
The first four rows are the plug-in GP prior with various kernels (Matérn
kernel, squared exponential kernel, second-order Sobolev kernel, and the
selected kernel via cross validation). The fourth and fifth methods are the
random series prior using B-splines and the inverse method proposed in
Holsclaw et al., (2013). The last row “Inverse*” excludes simulations
“Inverse” produces zero estimates. Methods with the smallest RMSE in each
column are boldfaced. Standard errors are provided in parentheses.

| $n=100$ | $n=500$
---|---|---
|  | $\sigma=0.1$
---
| $\sigma=0.2$
---
| $\sigma=0.3$
---
| $\sigma=0.1$
---
| $\sigma=0.2$
---
| $\sigma=0.3$
---
Matérn | 0.09 (0.02) | 0.13 (0.03) | 0.17 (0.03) | 0.06 (0.01) | 0.08 (0.02) | 0.10 (0.02)
SE | 0.10 (0.02) | 0.16 (0.02) | 0.20 (0.03) | 0.07 (0.01) | 0.10 (0.01) | 0.12 (0.02)
Sobolev | 0.06 (0.01) | 0.08 (0.01) | 0.10 (0.02) | 0.03 (0.01) | 0.05 (0.01) | 0.07 (0.01)
CV | 0.06 (0.01) | 0.10 (0.03) | 0.12 (0.05) | 0.03 (0.01) | 0.05 (0.01) | 0.07 (0.01)
Inverse | 0.11 (0.12) | 0.11 (0.10) | 0.13 (0.09) | 0.15 (0.11) | 0.12 (0.12) | 0.11 (0.12)
Inverse* | 0.07 (0.02) | 0.08 (0.03) | 0.10 (0.03) | 0.11 (0.04) | 0.08 (0.03) | 0.07 (0.02)
B-splines | 0.08 (0.01) | 0.11 (0.02) | 0.14 (0.03) | 0.04 (0.01) | 0.06 (0.01) | 0.08 (0.02)

We next compare a fully Bayesian approach with the empirical Bayes strategy
for the proposed method. We use the priors $\sigma^{2}\sim\text{Inverse-
Gamma}(20,1)$ and $\lambda\sim\text{Gamma}(1,1000)$ (parameterized using the
shape and rate parameters with mean $1/1000$) and draw $S=10000$ posterior
samples using the Metropolis-Hastings algorithm; the final estimate is the
average of 10000 posterior samples $\hat{f}_{n,s}^{\prime}$ for
$s=1,\ldots,S$. Similarly, Table 3 reports its RMSEs calculated at 100 equally
spaced points in $[0,10]$ over 100 repetitions, where we also include the
results of empirical Bayes shown in Table 2 to ease comparison. According to
Table 3, we can see that there is no significant difference between the two
treatments.

Table 3: RMSEs for estimating $f_{0}^{\prime}$ over 100 repetitions for
empirical Bayes (EB) and fully Bayesian (FB) approaches. Standard errors are
provided in parentheses.

| $n=100$ | $n=500$
---|---|---
|  | $\sigma=0.1$
---
| $\sigma=0.2$
---
| $\sigma=0.3$
---
| $\sigma=0.1$
---
| $\sigma=0.2$
---
| $\sigma=0.3$
---
Matérn (EB) | 0.09 (0.02) | 0.13 (0.03) | 0.17 (0.03) | 0.06 (0.01) | 0.08 (0.02) | 0.10 (0.02)
Matérn (FB) | 0.08 (0.01) | 0.12 (0.02) | 0.15 (0.03) | 0.05 (0.01) | 0.08 (0.01) | 0.09 (0.02)
SE (EB) | 0.10 (0.02) | 0.16 (0.02) | 0.20 (0.03) | 0.07 (0.01) | 0.10 (0.01) | 0.12 (0.02)
SE (FB) | 0.10 (0.02) | 0.16 (0.02) | 0.19 (0.03) | 0.06 (0.00) | 0.10 (0.01) | 0.13 (0.02)
Sobolev (EB) | 0.06 (0.01) | 0.08 (0.01) | 0.10 (0.02) | 0.03 (0.01) | 0.05 (0.01) | 0.07 (0.01)
Sobolev (FB) | 0.07 (0.01) | 0.08 (0.01) | 0.10 (0.02) | 0.05 (0.01) | 0.06 (0.01) | 0.07 (0.01)

## References

  * Banerjee et al., (2003) Banerjee, S., Gelfand, A. E., and Sirmans, C. (2003).  Directional rates of change under spatial process models.  Journal of the American Statistical Association, 98(464):946–954.
  * Berlinet and Thomas-Agnan, (2011) Berlinet, A. and Thomas-Agnan, C. (2011).  Reproducing Kernel Hilbert Spaces in Probability and Statistics.  Springer Science & Business Media.
  * Bhattacharya et al., (2014) Bhattacharya, A., Pati, D., and Dunson, D. (2014).  Anisotropic function estimation using multi-bandwidth Gaussian processes.  The Annals of Statistics, 42(1):352.
  * Bickel and Ritov, (2003) Bickel, P. J. and Ritov, Y. (2003).  Nonparametric estimators which can be “plugged-in”.  The Annals of Statistics, 31(4):1033–1053.
  * Cahill et al., (2015) Cahill, N., Kemp, A. C., Horton, B. P., and Parnell, A. C. (2015).  Modeling sea-level change using errors-in-variables integrated Gaussian processes.  The Annals of Applied Statistics, 9(2):547–571.
  * Castillo and Nickl, (2013) Castillo, I. and Nickl, R. (2013).  Nonparametric Bernstein–von Mises theorems in Gaussian white noise.  The Annals of Statistics, 41(4):1999–2028.
  * Castillo and Rousseau, (2015) Castillo, I. and Rousseau, J. (2015).  A Bernstein–von Mises theorem for smooth functionals in semiparametric models.  The Annals of Statistics, 43(6):2353–2383.
  * Charnigo et al., (2011) Charnigo, R., Hall, B., and Srinivasan, C. (2011).  A generalized ${C}_{p}$ criterion for derivative estimation.  Technometrics, 53(3):238–253.
  * Church and White, (2006) Church, J. A. and White, N. J. (2006).  A 20th century acceleration in global sea-level rise.  Geophysical Research Letters, 33(1).
  * Church and White, (2011) Church, J. A. and White, N. J. (2011).  Sea-level rise from the late 19th to the early 21st century.  Surveys in Geophysics, 32(4):585–602.
  * Conway, (1994) Conway, J. (1994).  A Course in Functional Analysis.  Graduate Texts in Mathematics. Springer New York.
  * Cucker and Zhou, (2007) Cucker, F. and Zhou, D. X. (2007).  Learning Theory: An Approximation Theory Viewpoint, volume 24.  Cambridge University Press.
  * Dai et al., (2018) Dai, X., Müller, H.-G., and Tao, W. (2018).  Derivative principal component analysis for representing the time dynamics of longitudinal and functional data.  Statistica Sinica, 28(3):1583–1609.
  * de Jonge and van Zanten, (2013) de Jonge, R. and van Zanten, H. (2013).  Semiparametric Bernstein–von Mises for the error standard deviation.  Electronic Journal of Statistics, 7:217–243.
  * Gelfand et al., (2003) Gelfand, A. E., Kim, H.-J., Sirmans, C., and Banerjee, S. (2003).  Spatial modeling with spatially varying coefficient processes.  Journal of the American Statistical Association, 98(462):387–396.
  * Gelfand and Schliep, (2016) Gelfand, A. E. and Schliep, E. M. (2016).  Spatial statistics and Gaussian processes: A beautiful marriage.  Spatial Statistics, 18:86–104.
  * Ghosal and van der Vaart, (2017) Ghosal, S. and van der Vaart, A. W. (2017).  Fundamentals of Nonparametric Bayesian Inference, volume 44.  Cambridge University Press.
  * Holsclaw et al., (2013) Holsclaw, T., Sansó, B., Lee, H. K., Heitmann, K., Habib, S., Higdon, D., and Alam, U. (2013).  Gaussian process modeling of derivative curves.  Technometrics, 55(1):57–67.
  * James et al., (2009) James, G. M., Wang, J., and Zhu, J. (2009).  Functional linear regression that’s interpretable.  The Annals of Statistics, 37(5A):2083–2108.
  * Jiang and Tokdar, (2021) Jiang, S. and Tokdar, S. T. (2021).  Variable selection consistency of Gaussian process regression.  The Annals of Statistics, 49(5):2491–2505.
  * Li and Dunson, (2020) Li, M. and Dunson, D. B. (2020).  Comparing and weighting imperfect models using D-probabilities.  Journal of the American Statistical Association, 115(531):1349–1360.
  * Liu and Li, (2022) Liu, Z. and Li, M. (2022).  On the estimation of derivatives using plug-in KRR estimators.  arXiv preprint arXiv:2006.01350.
  * McDougall and Stoner, (1938) McDougall, J. and Stoner, E. C. (1938).  The computation of fermi-dirac functions.  Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 237(773):67–104.
  * Pati and Bhattacharya, (2015) Pati, D. and Bhattacharya, A. (2015).  Adaptive Bayesian inference in the Gaussian sequence model using exponential-variance priors.  Statistics & Probability Letters, 103:100–104.
  * Rasmussen and Williams, (2006) Rasmussen, C. E. and Williams, C. K. (2006).  Gaussian Process for Machine Learning.  The MIT Press.
  * Riihimäki and Vehtari, (2010) Riihimäki, J. and Vehtari, A. (2010).  Gaussian processes with monotonicity information.  In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, pages 645–652.
  * Rivoirard and Rousseau, (2012) Rivoirard, V. and Rousseau, J. (2012).  Bernstein–von Mises theorem for linear functionals of the density.  The Annals of Statistics, 40(3):1489–1523.
  * Schott, (2016) Schott, J. R. (2016).  Matrix Analysis for Statistics.  John Wiley & Sons.
  * Shi and Choi, (2011) Shi, J. Q. and Choi, T. (2011).  Gaussian Process Regression Analysis for Functional Data.  CRC Press.
  * Smale and Zhou, (2005) Smale, S. and Zhou, D.-X. (2005).  Shannon sampling II: Connections to learning theory.  Applied and Computational Harmonic Analysis, 19(3):285–302.
  * Smale and Zhou, (2007) Smale, S. and Zhou, D.-X. (2007).  Learning theory estimates via integral operators and their approximations.  Constructive Approximation, 26(2):153–172.
  * Stein, (2012) Stein, M. L. (2012).  Interpolation of Spatial Data: Some Theory for Kriging.  Springer Science & Business Media.
  * Steinwart and Christmann, (2008) Steinwart, I. and Christmann, A. (2008).  Support Vector Machines.  Springer Science & Business Media.
  * Stone, (1982) Stone, C. J. (1982).  Optimal global rates of convergence for nonparametric regression.  The Annals of Statistics, 10(4):1040–1053.
  * van der Vaart and Wellner, (1996) van der Vaart, A. and Wellner, J. (1996).  Weak Convergence and Empirical Processes: With Applications to Statistics.  Springer Series in Statistics. Springer.
  * van der Vaart and van Zanten, (2009) van der Vaart, A. W. and van Zanten, J. H. (2009).  Adaptive Bayesian estimation using a Gaussian random field with inverse gamma bandwidth.  The Annals of Statistics, 37(5B):2655–2675.
  * Wahba, (1990) Wahba, G. (1990).  Spline Models for Observational Data, volume 59.  Siam.
  * Wahba and Wang, (1990) Wahba, G. and Wang, Y. (1990).  When is the optimal regularization parameter insensitive to the choice of the loss function?  Communications in Statistics-Theory and Methods, 19(5):1685–1700.
  * Wang and Berger, (2016) Wang, X. and Berger, J. O. (2016).  Estimating shape constrained functions using Gaussian processes.  SIAM/ASA Journal on Uncertainty Quantification, 4(1):1–25.
  * Wang et al., (2022) Wang, Z., Magnotti, J., Beauchamp, M. S., and Li, M. (2022).  Functional group bridge for simultaneous regression and support estimation.  Biometrics.  In press.
  * Yoo and Ghosal, (2016) Yoo, W. W. and Ghosal, S. (2016).  Supremum norm posterior contraction and credible sets for nonparametric multivariate regression.  The Annals of Statistics, 44(3):1069–1102.
  * Yu et al., (2022) Yu, C.-H., Li, M., Noe, C., Fischer-Baum, S., and Vannucci, M. (2022).  Bayesian inference for stationary points in Gaussian process regression models for event-related potentials analysis.  Biometrics.  In press.
  * Zhang, (2005) Zhang, T. (2005).  Learning bounds for kernel regression using effective data dimensionality.  Neural Computation, 17(9):2077–2098.