#  Comparing Scale Parameter Estimators for Gaussian Process Interpolation
with the Brownian Motion Prior: Leave-One-Out Cross Validation and Maximum
Likelihood

Masha Naslidnyk Department of Statistical Science, University College London,
UK. Motonobu Kanagawa Data Science Department, EURECOM, France. Toni Karvonen
Department of Mathematics and Statistics, University of Helsinki, Finland.
Maren Mahsereci University of Tübingen, Germany.

###### Abstract

Gaussian process (GP) regression is a Bayesian nonparametric method for
regression and interpolation, offering a principled way of quantifying the
uncertainties of predicted function values. For the quantified uncertainties
to be well-calibrated, however, the kernel of the GP prior has to be carefully
selected. In this paper, we theoretically compare two methods for choosing the
kernel in GP regression: cross-validation and maximum likelihood estimation.
Focusing on the scale-parameter estimation of a Brownian motion kernel in the
noiseless setting, we prove that cross-validation can yield asymptotically
well-calibrated credible intervals for a broader class of ground-truth
functions than maximum likelihood estimation, suggesting an advantage of the
former over the latter. Finally, motivated by the findings, we propose
_interior cross validation_ , a procedure that adapts to an even broader class
of ground-truth functions.

###### Contents

  1. 1 Introduction
    1. 1.1 Scale Parameter Estimation
    2. 1.2 Contributions
    3. 1.3 Related work
  2. 2 Background
    1. 2.1 Gaussian process regression
    2. 2.2 Kernel parameter estimation
  3. 3 Setting
    1. 3.1 Brownian motion kernel
    2. 3.2 Sequences of partitions
    3. 3.3 Hölder spaces
    4. 3.4 Fractional Brownian motion
    5. 3.5 Functions of finite quadratic variation
  4. 4 Main results
    1. 4.1 Deterministic setting
    2. 4.2 Random setting
    3. 4.3 Interior cross-validation estimators
  5. 5 Consequences for credible intervals
  6. 6 Experiments
    1. 6.1 Test functions
    2. 6.2 Asymptotics of the CV estimator
    3. 6.3 Comparison of CV and ML estimators
  7. 7 Conclusion and future work
  8. 8 Proofs
    1. 8.1 Explicit expressions for the CV and ML estimators
    2. 8.2 Proofs for Section 4.1
    3. 8.3 Proofs for Section 4.2
    4. 8.4 Proofs for Section 4.3
    5. 8.5 Proofs for Section 5
  9. A Connection between the ML and CV estimators
  10. B Further discussion on Theorem 10
  11. C Explicit expression for the leave-$p$-out estimator
  12. D Comparison of CV and ML estimators for Matérn kernels

##  1 Introduction

Gaussian process (GP) regression (or kriging) is a Bayesian nonparametric
method for regression and interpolation that has been extensively studied in
statistics and machine learning (O’Hagan,, 1978; Stein,, 1999; Rasmussen and
Williams,, 2006). Its key property is that it enables uncertainty
quantification of estimated function values in a principled manner, which is
crucial for applications involving decision-making, safety concerns, and
scientific discovery. As such, GP regression has been a core building block of
more applied algorithms, including Bayesian optimisation (Jones et al.,, 1998;
Shahriari et al.,, 2015; Garnett,, 2023), probabilistic numerical computation
(Hennig et al.,, 2015; Cockayne et al.,, 2019; Hennig et al.,, 2022), and
calibration and emulation of computer models  (Sacks et al.,, 1989; Kennedy
and O’Hagan,, 2001; O’Hagan,, 2006; Beck and Guillas,, 2016; Gu et al.,,
2018), to name just a few.

GP regression estimates an unknown function $f$ from its observations as
follows. One first defines a prior distribution for $f$ as a GP by specifying
its kernel (and mean function). Provided $N$ observations about $f$, one then
derives the posterior distribution of $f$, which is another GP with mean
function $m_{N}$ and kernel (or covariance function) $k_{N}$. One can then
predict the function value $f(x)$ at any input $x$ by the posterior mean
$m_{N}(x)$ and quantify its uncertainty using the posterior standard deviation
$\sqrt{\smash[b]{k_{N}(x)}}\coloneqq\sqrt{\smash[b]{k_{N}(x,x)}}$.
Specifically, one can construct a credible interval of $f(x)$ as the interval
$[m_{N}(x)-\alpha\sqrt{\smash[b]{k_{N}(x)}},m_{N}(x)+\alpha\sqrt{\smash[b]{k_{N}(x)}}]$
for a constant $\alpha>0$ (for example, $\alpha\approx 1.96$ leads to the 95%
credible interval). Such uncertainty estimates constitute key ingredients in
the above applications of GP regression.

For GP uncertainty estimates to be reliable, the posterior standard deviation
$\sqrt{\smash[b]{k_{N}(x)}}$ should, ideally, decay at the same rate as the
prediction error $|m_{N}(x)-f(x)|$ decreases, with the increase of sample size
$N$. Otherwise, GP uncertainty estimates are either asymptotically
overconfident or underconfident. For example, if $\sqrt{\smash[b]{k_{N}(x)}}$
goes to $0$ faster than the error $|m_{N}(x)-f(x)|$, then the credible
interval
$[m_{N}(x)-\alpha\sqrt{k_{N}(x)},m_{N}(x)+\alpha\sqrt{\smash[b]{k_{N}(x)}}]$
will not contain the true value $f(x)$ as $N$ increases for any fixed constant
$\alpha>0$ (asymptotically overconfident). If $\sqrt{\smash[b]{k_{N}(x)}}$
goes to $0$ slower than the error $|m_{N}(x)-f(x)|$, then the confidence
interval
$[m_{N}(x)-\alpha\sqrt{\smash[b]{k_{N}(x)}},m_{N}(x)+\alpha\sqrt{\smash[b]{k_{N}(x)}}]$
will get larger than the error $|m_{N}(x)-f(x)|$ as $N$ increases
(asymptotically underconfident). Both of these cases are not desirable in
practice, as GP credible intervals will not be accurate estimates of
prediction errors.

Figure 1: GP interpolation of a fractional Brownian motion with the Hurst
parameter $H=0.2$ (smoothness $l+\alpha=0.2$) using the Brownian motion kernel
(2) with three different scale parameters: $\sigma^{2}=1$ (left),
$\sigma^{2}=\hat{\sigma}_{\rm CV}^{2}=4.752$ given by the LOO-CV estimator
(middle) and $\sigma^{2}=\hat{\sigma}_{\rm ML}^{2}=3.729$ obtained with the ML
estimator (right). In each figure, the red trajectory represents the path of
the fractional Brownian motion, the purple circles the training data, the blue
curve the posterior mean $m_{N}(x)$ and the green shade the 95 % credible
interval
$[m_{N}(x)-1.96\sigma\sqrt{\smash[b]{k_{N}(x)}},m_{N}(x)+1.96\sigma\sqrt{\smash[b]{k_{N}(x)}}]$.
Figure 2: GP interpolation of an integrated fractional Brownian motion with
the Hurst parameter $H=0.5$ (smoothness $l+\alpha=1.5$) using the Brownian
motion kernel (2) with three different scale parameters: $\sigma^{2}=1$
(left), $\sigma^{2}=\hat{\sigma}_{\rm CV}^{2}=0.019$ given by the LOO-CV
estimator (middle) and $\sigma^{2}=\hat{\sigma}_{\rm ML}^{2}=0.067$ obtained
with the ML estimator (right). For the explanation of the figures, see the
caption of Figure 1.

Unfortunately, in general, the posterior standard deviation
$\sqrt{\smash[b]{k_{N}(x)}}$ does not decay at the same rate as the prediction
error $|f(x)-m_{N}(x)|$, because, as is well-known,
$\sqrt{\smash[b]{k_{N}(x)}}$ does not depend on the true function $f$; see
(3b) in Section 2.1. Exceptionally, if the function $f$ is a sample path of
the GP prior (the well-specified case), GP uncertainty estimates can be well-
calibrated. However, in general, the unknown $f$ is not exactly a sample path
of the GP prior (the misspecified case), and the posterior standard deviation
$\sqrt{\smash[b]{k_{N}(x)}}$ does not scale with the prediction error
$|f(x)-m_{N}(x)|$. Figures 1 and 2 (the left panels) show examples where the
true function $f$ is not a sample of the GP prior and where the GP uncertainty
estimates are not well-calibrated.

###  1.1 Scale Parameter Estimation

To obtain sensible uncertainty estimates, one thus needs to adapt the
posterior standard deviation $\sqrt{\smash[b]{k_{N}(x)}}$ to the function $f$.
One simple way to achieve this is to introduce the scale parameter
$\sigma^{2}>0$ and parametrize the kernel as

$k_{\sigma}(x,x^{\prime})\coloneqq\sigma^{2}k(x,x^{\prime}),$ (1)

where $k$ is the original kernel. GP regression with this kernel $k_{\sigma}$
yields the posterior mean function $m_{N}$, which is not influenced by
$\sigma^{2}$, and the posterior covariance function $\sigma^{2}k$, which is
scaled by $\sigma^{2}$. If one estimates $\sigma^{2}$ from observed data of
$f$, the estimate $\hat{\sigma}^{2}$ depends on $f$, and so does the resulting
posterior standard deviation $\hat{\sigma}\sqrt{\smash[b]{k_{N}(x)}}$.

One approach to scale-parameter estimation is the method of maximum likelihood
(ML), which optimizes $\sigma^{2}$ to maximize the marginal likelihood of the
observed data (Rasmussen and Williams,, 2006, Section 5.4). The ML approach is
popular for general hyperparameter optimization in GP regression. Another less
common way in the GP literature is cross-validation (CV), which optimizes
$\sigma^{2}$ to maximize the average predictive likelihood of held-out data
(Sundararajan and Keerthi,, 2001). For either approach, the optimized scale
parameter can be obtained analytically in computational complexity
$\mathcal{O}(N^{3})$. Figures 1 and 2 (middle and right panels) demonstrate
that both approaches yield uncertainty estimates better calibrated than the
original estimates without the scale parameter.

Do these scale parameter estimators lead to asymptotically well-calibrated
uncertainty estimates? To answer this question, one needs to understand their
convergence properties as the sample size $N$ increases. Most existing
theoretical works focus on the well-specified case where there is a “true”
scale parameter $\sigma_{0}^{2}$ such that the unknown $f$ is a GP with the
kernel $\sigma_{0}^{2}k$. In this case, both the ML and CV estimators have
been shown to be consistent in estimating the true $\sigma_{0}^{2}$ (e.g.,
Ying,, 1991; Zhang,, 2004; Bachoc et al.,, 2017, 2020). However, in general,
no “true” scale parameter $\sigma^{2}_{0}$ exists such that the unknown $f$ is
a GP with the covariance $\sigma_{0}^{2}k$. In such misspecified cases, not
much is known about the convergence properties of both estimators. Karvonen et
al., (2020) analyze the ML estimator for the scale parameter, assuming that
$f$ is a deterministic function. They derive upper bounds (and lower bounds in
some cases) for the ML estimator; see Wang, (2021) for closely related work.
To our knowledge, no theoretical work exists for the CV estimator for the
scale parameter in the misspecified case. Bachoc, (2013) and Petit et al.,
(2022) empirically compare the ML and CV estimators under different model
misspecification settings. We will review other related works in Section 1.3.

###  1.2 Contributions

This work studies the convergence properties of the ML and CV estimators,
$\hat{\sigma}_{\rm ML}^{2}$ and $\hat{\sigma}_{\rm CV}^{2}$, of the scale
parameter $\sigma^{2}$ in GP regression, to understand whether they lead to
asymptotically well-calibrated uncertainty estimates. In particular, we
provide the first theoretical analysis of the CV estimator $\hat{\sigma}_{\rm
CV}^{2}$ when the GP prior is misspecified, and also establish novel results
for the ML estimator $\hat{\sigma}_{\rm ML}^{2}$.

To facilitate the analysis, we focus on the following simplified setting. For
a constant $T>0$, let $[0,T]\subset\mathbb{R}$ be the input domain. Let $k$ in
(1) be the Brownian motion kernel

$k(x,x^{\prime})=\min(x,x^{\prime})\quad\text{ for }\quad
x,x^{\prime}\in[0,T].$ (2)

With this choice, a sample path of the GP prior has roughly a smoothness of
1/2 (in terms of the differentiability; we will be more rigorous in later
sections).

We assume that the true unknown function $f$ has the smoothness $l+\alpha$,
where $l\in\\{0\\}\cup\mathbb{N}$ and $0<\alpha\leq 1$. The GP prior has well-
specified smoothness if $l=0$ and $\alpha=1/2$. Other settings of $l$ and
$\alpha$ represent misspecified cases. If $l=0$ and $\alpha<1/2$, the true
function $f$ is rougher than the GP prior (Figure 1); if $l=0$ and
$\alpha>1/2$ or $l\geq 1$, the function $f$ is smoother than the GP prior. We
focus on the noise-free setting where one observes the function values
$f(x_{1}),\dots,f(x_{N})$ at input points $x_{1},\dots,x_{N}\in[0,T]$.

Our main results are new upper and lower bounds for the asymptotic rates of
the CV estimator $\hat{\sigma}_{\rm CV}^{2}$ and the ML estimator
$\hat{\sigma}_{\rm ML}^{2}$ as $N\to\infty$ (Section 4). The results suggest
that the CV estimator can yield asymptotically well-calibrated uncertainty
estimates for a broader class of functions $f$ than the ML estimator; thus,
the former has an advantage over the latter (Section 5). More specifically,
asymptotically well-calibrated uncertainty estimates may be obtained with the
CV estimator for the range $0<l+\alpha\leq 3/2$ of smoothness of the true
function, while this range becomes $0<l+\alpha\leq 1$ with the ML estimator
and is narrower. This finding is consistent with the example in Figure 2,
where the true function has smoothness $l+\alpha=3/2$ and is thus smoother
than the GP prior. The uncertainty estimates of the CV estimator appear to be
well-calibrated, while those of the ML estimator are unnecessarily wide,
failing to adapt to the smoothness. Motivated by these insights, we propose a
method called _interior cross-validation_ , and show it accommodates an even
wider range of smoothness of the true function than the CV estimator.

This paper is structured as follows. After reviewing related works in Section
1.3, we introduce the necessary background on the ML and CV approaches to
scale parameter estimation for GP regression in Section 2. We describe the
setting of the theoretical analysis in Section 3, present our main results in
Section 4, and discuss its consequences on uncertainty quantification in
Section 5. We report simulation experiments in Section 6, conclude in Section
7, and present proofs in Section 8.

###  1.3 Related work

We review here related theoretical works on hyper-parameter selection in GP
regression in the noiseless setting. We categorize them into two groups based
on how the true unknown function $f$ is modelled: random and deterministic.

#### Random setting.

One group of works models the ground truth $f$ as a random function,
specifically as a GP. Most of these works model $f$ as a GP with a Matérn-type
kernel and analyze the ML estimator. Under the assumption that the GP prior is
correctly specified, asymptotic properties of the ML estimator for the scale
parameter and other parameters have been studied (Stein,, 1990; Ying,, 1991,
1993; Loh and Kam,, 2000; Zhang,, 2004; Loh,, 2005; Du et al.,, 2009;
Anderes,, 2010; Wang and Loh,, 2011; Kaufman and Shaby,, 2013; Bevilacqua et
al.,, 2019). Recently Loh et al., (2021) and Loh and Sun, (2023) have
constructed consistent estimators of various parameters for many commonly used
kernels, including Matérns. Chen et al., (2021) and Petit, (2023) consider a
periodic version of Matérn GPs, and show the consistency of the ML estimator
for its smoothness parameter. To our knowledge, the only existing theoretical
result for ML estimation of the scale parameter in the misspecified random
setting considers oversmoothing (Karvonen,, 2021, Theorem 4.2). Oversmoothing
refers to the situation where the chosen kernel is smoother than the true
function. In Section 4.2 (Theorem 12), we provide a result for the
undersmoothing case, which occurs when the chosen kernel is less smooth than
the true function.

In contrast, few theoretical works exist for the CV estimator. Bachoc et al.,
(2017) study the leave-one-out (LOO) CV estimator for the Matérn-1/2 model (or
the Laplace kernel) with one-dimensional inputs, in which case the GP prior is
an Ornstein–Uhlenbeck (OU) process. Assuming the well-specified case where the
true function is also an OU process, they prove the consistency and asymptotic
normality of the CV estimator for the microergodic parameter in the fixed-
domain asymptotic setting. Bachoc, (2018) and Bachoc et al., (2020) discuss
another CV estimator that uses the mean square prediction error as the scoring
criterion of CV (thus different from the one discussed here) in the
increasing-domain asymptotics. Bachoc, (2013) and Petit et al., (2022) perform
empirical comparisons of the ML and CV estimators under different model
misspecification settings. Thus, to our knowledge, no theoretical result
exists for the CV estimator of the scale parameter in the random misspecified
setting, which we provide in Section 4.2 (Theorem 11).

#### Deterministic setting.

Another line of research assumes that the ground truth $f$ is a fixed function
belonging to a specific function space (Stein,, 1993). Xu and Stein, (2017)
assumed that the ground truth $f$ is a monomial on $[0,1]$ and proved some
asymptotic results for the ML estimator when the kernel $k$ is Gaussian. As
mentioned earlier, Karvonen et al., (2020) proved asymptotic upper (and, in
certain cases, also lower) bounds on the ML estimator $\hat{\sigma}_{\rm
ML}^{2}$ of the scale parameter $\sigma^{2}$; see Wang, (2021) for a closely
related work. Karvonen, (2023) has studied the ML and LOO-CV estimators for
the smoothness parameter in the Matérn model; see also Petit, (2023). Ben
Salem et al., (2019) and Karvonen and Oates, (2023) proved non-asymptotic
results on the length-scale parameter in the Matérn and related models. Thus,
there has been no work for the CV estimator of the scale parameter
$\sigma^{2}$ in the deterministic setting, which we provide in Section 4.1
(Theorem 7); we also prove a corresponding result for the ML estimator
(Theorem 8).

##  2 Background

This section briefly reviews GP regression and the ML and LOO-CV estimators of
kernel parameters.

###  2.1 Gaussian process regression

We first explain GP regression (or interpolation). Let $\Omega$ be a set, and
$f\colon\Omega\to\mathbb{R}$ be an unknown function of interest. Suppose one
observes $N$ function values $f(x_{1}),\dots,f(x_{N})$ at pairwise distinct
input points $x_{1},\dots,x_{N}\in\Omega$. The task here is to estimate $f$
based on the data $(\mathbf{x},f(\mathbf{x}))$, where
$f(\mathbf{x})\coloneqq[f(x_{1}),\dots,f(x_{N})]^{\top}\in\mathbb{R}^{N}$ and
$\mathbf{x}\coloneqq[x_{1},\ldots,x_{N}]^{\top}\in\Omega^{N}$.

In GP regression, one first defines a prior distribution of the unknown $f$ as
a GP by specifying its mean function $m\colon\Omega\to\mathbb{R}$ and
covariance function (kernel) $k\colon\Omega\times\Omega\to\mathbb{R}$; we may
write $f\sim\mathcal{GP}(m,k)$ to indicate this. Conditioned on the data
$(\mathbf{x},f(\mathbf{x}))$, the posterior distribution of $f$ is again a GP
whose mean function $m_{N}:\Omega\to\mathbb{R}$ and covariance function
$k_{N}:\Omega\times\Omega\to\mathbb{R}$ are given by

$\displaystyle m_{N}(x)$ $\displaystyle\coloneqq
m(x)+k(x,\mathbf{x})^{\top}k(\mathbf{x},\mathbf{x})^{-1}\left(f(\mathbf{x})-m(\mathbf{x})\right),\quad
x\in\Omega,$ (3a) $\displaystyle k_{N}(x,x^{\prime})$ $\displaystyle\coloneqq
k(x,x^{\prime})-k(x,\mathbf{x})^{\top}k(\mathbf{x},\mathbf{x})^{-1}k(x^{\prime},\mathbf{x}),\quad
x,x^{\prime}\in\Omega,$ (3b)

where
$m(\mathbf{x})\coloneqq[m(x_{1}),\ldots,m(x_{N})]^{\top}\in\mathbb{R}^{N}$ and
$k(x,\mathbf{x})\coloneqq[k(x,x_{1}),\ldots,k(x,x_{N})]^{\top}\in\mathbb{R}^{N}$,
and

$k(\mathbf{x},\mathbf{x})\coloneqq\begin{bmatrix}k(x_{1},x_{1})&\dots&k(x_{1},x_{N})\\\
\vdots&\ddots&\vdots\\\
k(x_{N},x_{1})&\dots&k(x_{N},x_{N})\end{bmatrix}\in\mathbb{R}^{N\times N}$ (4)

is the Gram matrix. Throughout this paper, we assume that the points
$\mathbf{x}$ are such that the Gram matrix is non-singular. For notational
simplicity, we may write the posterior variance as

$k_{N}(x)\coloneqq k_{N}(x,x),\quad x\in\Omega.$

For simplicity and as commonly done, we henceforth assume that the prior mean
function $m$ is the zero function, $m(\cdot)\equiv 0$.

While the GP prior assumes that the unknown function $f$ is a sample path of
the GP with the specified kernel $k$, this assumption does not hold in
general, i.e., model misspecification occurs. In this case, as described in
Figures 1 and 2 (left), the posterior standard deviation
$\sqrt{\smash[b]{k_{N}(x)}}$, which is supposed to quantify the uncertainty of
the unknown function value $f(x)$, may not be well calibrated with the
prediction error $|m_{N}(x)-f(x)|$. One could address this issue by selecting
the kernel $k$ or its parameters from the data $(\mathbf{x},f(\mathbf{x}))$;
we will explain this topic next.

###  2.2 Kernel parameter estimation

The selection of the kernel $k$ is typically performed by defining a
parametric family of kernels $\\{k_{\theta}\\}_{\theta\in\Theta}$ and
selecting the parameter $\theta$ based on an appropriate criterion. Here
$\Theta$ is a parameter set, and $k_{\theta}:\Omega\times\Omega\to\mathbb{R}$
for each $\theta\in\Theta$ is a kernel.

#### Maximum likelihood (ML) estimation.

The ML estimator maximises the log-likelihood of the data
$(\mathbf{x},f(\mathbf{x}))$ given that $f$ is a GP with kernel $k_{\theta}$:

$\log
p(f(\mathbf{x})\,|\,\mathbf{x},\theta)=-\frac{1}{2}\bigg{(}f(\mathbf{x})^{\top}k_{\theta}(\mathbf{x},\mathbf{x})^{-1}f(\mathbf{x})+\log\det
k_{\theta}(\mathbf{x},\mathbf{x})+n\log(2\pi)\bigg{)},$

where $\det k_{\theta}(\mathbf{x},\mathbf{x})$ is the determinant of the Gram
matrix $k_{\theta}(\mathbf{x},\mathbf{x})$ (see, e.g., Rasmussen and Williams,
2006, Section 5.4.1). With the additive terms that do not depend on $\theta$
removed from $\log p(f(\mathbf{x})\,|\,\mathbf{x},\theta)$, this is equivalent
to minimising the loss function

$\mathcal{L}_{\textup{ML}}(\theta):=f(\mathbf{x})^{\top}k_{\theta}(\mathbf{x},\mathbf{x})^{-1}f(\mathbf{x})+\log\det
k_{\theta}(\mathbf{x},\mathbf{x}).$ (5)

In general, $\mathcal{L}_{\textup{ML}}(\theta)$ may not have a unique
minimiser, so that any ML estimator satisfies

$\hat{\theta}_{\textup{ML}}\in\mathop{\text{arg\,min}}_{\theta\in\Theta}\mathcal{L}_{\textup{ML}}(\theta).$

#### Leave-one-out cross-validation (LOO-CV).

The LOO-CV estimator (e.g., Rasmussen and Williams,, 2006, Section 5.4.2),
which we may simply call the CV estimator, is an alternative to the ML
estimator. It maximizes the average log-predictive likelihood

$\sum_{n=1}^{N}\log p(f(x_{n})\,|\,x_{n},\mathbf{x}_{\setminus
n},f(\mathbf{x}_{\setminus n}),\theta)$ (6)

of the held-out data $(x_{n},f(x_{n}))$, where $n=1,\dots,N$, based on the
data $(\mathbf{x}_{\setminus n},f(\mathbf{x}_{\setminus n}))$, where
$\mathbf{x}_{\setminus n}$ denotes the input points with $x_{n}$ removed:

$\mathbf{x}_{\setminus
n}=[x_{1},\dots,x_{n-1},x_{n+1},\dots,x_{N}]^{\top}\in\Omega^{N-1}.$

Let $m_{\theta,\setminus n}$ and $k_{\theta,\setminus n}$ denote the posterior
mean and covariance functions of GP regression with the kernel $k_{\theta}$
and the data $(\mathbf{x}_{\setminus n},f(\mathbf{x}_{\setminus n})$. Because
each $p(f(x_{n})\,|\,x_{n},\mathbf{x}_{\setminus n},f(\mathbf{x}_{\setminus
n}),\theta)$ is the Gaussian density of $f(x_{n})$ with mean
$m_{\theta,\setminus n}(x_{n})$ and variance $k_{\theta,\setminus
n}(x_{n})\coloneqq k_{\theta,\setminus n}(x_{n},x_{n})$, removing additive
terms that do not depend on $\theta$ and reversing the sign in (6) yields the
following CV objective function:

$\mathcal{L}_{\textup{CV}}(\theta)=\sum_{n=1}^{N}\frac{\left[f(x_{n})-m_{\theta,\setminus
n}(x_{n})\right]^{2}}{k_{\theta,\setminus n}(x_{n})}+\log k_{\theta,\setminus
n}(x_{n}){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>(7)

The CV estimator is then defined as its minimizer:

$\hat{\theta}_{\textup{CV}}\in\mathop{\text{arg\,min}}_{\theta\in\Theta}\mathcal{L}_{\textup{CV}}(\theta).$

As for the ML estimator, the CV objective function and its first-order
gradients can be computed in closed form in $\mathcal{O}(N^{3})$ time
(Sundararajan and Keerthi,, 2001).

#### Scale parameter estimation.

As explained in Section 1, we consider the family of kernels
$k_{\sigma}(x,x^{\prime})\coloneqq\sigma^{2}k(x,x^{\prime})$ parametrized with
the scale parameter
$\sigma^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{>0}}$,
where $k$ is a fixed kernel, and study the estimation of $\sigma^{2}$ using
the CV and ML estimators, denoted as $\hat{\sigma}_{\rm CV}^{2}$ and
$\hat{\sigma}_{\rm ML}^{2}$, respectively. In this case, both
$\hat{\sigma}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\rm
ML}}}^{2}$ and $\hat{\sigma}_{\rm CV}^{2}$ can be derived in closed form by
differentiating (5) and (7).

Let $m_{n-1}$ and $k_{n-1}$ be the posterior mean and variance functions of GP
regression using the kernel $k$ and the first $n-1$ training observations
$(x_{1},f(x_{1})),\dots,(x_{n-1},f(x_{n-1}))$. Let $m_{0}(\cdot)\coloneqq 0$
and $k_{0}(x,x)\coloneqq k(x,x)$. Then the ML estimator is given by

$\hat{\sigma}_{\textup{ML}}^{2}=\frac{f(\mathbf{x})^{\top}k(\mathbf{x},\mathbf{x})^{-1}f(\mathbf{x})}{N}=\frac{1}{N}\sum_{n=1}^{N}\frac{[f(x_{n})-m_{n-1}(x_{n})]^{2}}{k_{n-1}(x_{n})}.$
(8)

This expression of the ML estimator is relatively well known; see e.g. Section
4.2.2 in Xu and Stein, (2017) or Proposition 7.5 in Karvonen and Oates,
(2023).

On the other hand, the CV estimator
$\hat{\sigma}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\rm
CV}}}^{2}$ is given by

$\displaystyle\hat{\sigma}_{\textup{CV}}^{2}=\frac{1}{N}\sum_{n=1}^{N}\frac{\left[f(x_{n})-m_{\setminus
n}(x_{n})\right]^{2}}{k_{\setminus n}(x_{n})},$ (9)

where $m_{\backslash n}$ and $k_{\backslash n}$ are the posterior mean and
covariance functions of GP regression using the kernel $k$ and data
$(\mathbf{x}_{\setminus n},f(\mathbf{x}_{\setminus n}))$ with
$(x_{n},f(x_{n}))$ removed:

$\displaystyle m_{\setminus n}(x)$ $\displaystyle=k(\mathbf{x}_{\setminus
n},x)^{\top}k(\mathbf{x}_{\setminus n},\mathbf{x}_{\setminus
n})^{-1}f(\mathbf{x}_{\setminus n}),$ $\displaystyle k_{\setminus
n}(x,x^{\prime})$ $\displaystyle=k(x,x^{\prime})-k(\mathbf{x}_{\setminus
n},x)^{\top}k(\mathbf{x}_{\setminus n},\mathbf{x}_{\setminus
n})^{-1}k(\mathbf{x}_{\setminus n},x^{\prime}).$

Notice the similarity between the two expressions (8) and (9). The difference
is that the ML estimator uses $k_{n-1}$ and $m_{n-1}$, which are based on the
first $n-1$ training observations, while the CV estimator uses $k_{\backslash
n}$ and $m_{\backslash n}$ obtained with $N-1$ observations, for each
$n=1,\dots,N$. Therefore, the CV estimator uses all the data points more
evenly than the ML estimator. This difference may be the source of the
difference in their asymptotic properties established later.

######  Remark 1.

As suggested by the similarity between (8) and (9), there is a deeper
connection between ML and CV estimators in general. For instance, Fong and
Holmes, (2020, Proposition 2) have shown that the Bayesian marginal likelihood
equals the average of leave-$p$-out CV scores. We prove this result for the
special case of scale parameter estimation in GP regression in Appendix A.
Another notable example is the work in Ginsbourger and Schärer, (2021), where
the authors showed that, when corrected for the covariance of residuals, the
CV estimator of the scale parameter reverts to MLE.

##  3 Setting

This section describes the settings and tools for our theoretical analysis:
the Brownian motion kernel in Section 3.1; sequences of partitions in Section
3.2; the Hölder class of functions in Section 3.3; fractional Brownian motion
in Section 3.4; and functions of finite quadratic variation in Section 3.5.

###  3.1 Brownian motion kernel

As explained in Section 1, for the kernel $k$ we focus on the Brownian motion
kernel on the domain $\Omega=[0,T]$ for some $T>0$:

$k(x,x^{\prime})=\min(x,x^{\prime}).$

The resulting kernel $k_{\sigma}(x,x^{\prime})=\sigma^{2}k(x,x^{\prime})$
induces a Brownian motion prior for GP regression.

We assume that the input points $\mathbf{x}=[x_{1},\dots x_{N}]^{\top}$ for GP
regression are positive and ordered:

$0<x_{1}<x_{2}<\dots<x_{N}\leq T.$

The positivity ensures that the Gram matrix (4) is non-singular; the proof is
given in Section 8.1.

As is well known and can be seen in Figures 1 and 2, the posterior mean
function $m_{N}$ in (3) using the Brownian motion kernel becomes the piecewise
linear interpolant of the observations $(\mathbf{x},f(\mathbf{x}))$. See (25)
and (26) in Section 8.1 for the proof and explicit expressions of the
posterior mean and covariance functions.

###  3.2 Sequences of partitions

For our asymptotic analysis, we assume that the input points
$x_{1},\dots,x_{N}\in[0,T]$ cover the domain $[0,T]$ more densely as the
sample size $N$ increases. To make the dependence on the size $N$ explicit, we
write $\mathcal{P}_{N}\coloneqq(x_{N,n})_{n=1}^{N}\subset[0,T]$ as a point set
of size $N$, and assume that they are ordered as

$0\eqqcolon x_{N,0}<x_{N,1}<x_{N,2}<\dots<x_{N,N}=T$

Then $\mathcal{P}_{N}$ defines a partition of $[0,T]$ into $N$ subintervals
$[x_{N,n},x_{N,n+1}]$. When there is no risk of confusion, we may write
$x_{n}$ instead of $x_{N,n}$ for simplicity. Note that we do not require the
nesting $\mathcal{P}_{N}\subset\mathcal{P}_{N+1}$ of partitions.

We define the mesh size of partition $\mathcal{P}_{N}$ as the longest
subinterval in the partition:

$\|\mathcal{P}_{N}\|\coloneqq\max_{n\in\\{0,1,\ldots,N-1\\}}(x_{N,n+1}-x_{N,n})$

The decay rate of the mesh size $\|\mathcal{P}_{N}\|$ quantifies how quickly
the points in $\mathcal{P}_{N}$ cover the interval $[0,T]$. In particular, the
decay rate $\mathcal{P}_{N}=\mathcal{O}(N^{-1})$ implies that the length of
every subinterval is asymptotically upper bounded by $1/N$. At the same time,
if each subinterval is asymptotically lower bounded by $1/N$, we call the
sequence of partitions $(\mathcal{P}_{N})_{N\in\mathbb{N}}$ quasi-uniform,
more formally defined in Wendland, (2005, Definition 4.6) as follows.

######  Definition 2.

For each $N\in\mathbb{N}$, let
$\mathcal{P}_{N}\coloneqq(x_{N,n})_{n=1}^{N}\subset[0,T]$. Define $\Delta
x_{N,n}\coloneqq x_{N,n+1}-x_{N,n}$. Then the sequence of partitions
$(\mathcal{P}_{N})_{N\in\mathbb{N}}$ is called quasi-uniform if there exists a
constant $1\leq C_{\textup{qu}}<\infty$ such that

$\sup_{N\in\mathbb{N}}\frac{\max_{n}\Delta x_{N,n}}{\min_{n}\Delta
x_{N,n}}=C_{\textup{qu}}.$

Quasi-uniformity, as defined here, requires that the ratio of the longest
subinterval, $\max_{n}\Delta x_{N,n}$, to the shortest one, $\min_{n}\Delta
x_{N,n}$, is upper-bounded by $C_{\rm qu}$ for all $N\in\mathbb{N}$. Since
$\min_{n}\Delta x_{N,n}\leq TN^{-1}$ and $\max_{n}\Delta x_{N,n}\geq TN^{-1}$
for any partition of $[0,T]$, quasi-uniformity implies that all subintervals
are asymptotically upper and lower bounded by $1/N$, as we have, for all
$N\in\mathbb{N}$ and
$n_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{0}}}\in\\{0,\dots,N-1\\}$,

$\frac{TN^{-1}}{C_{\textup{qu}}}\leq\min_{n}\Delta x_{N,n}\leq\Delta
x_{N,n_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{0}}}}\leq\max_{n}\Delta
x_{N,n}\leq TC_{\textup{qu}}N^{-1}.$ (10)

Therefore, quasi-uniform sequences of partitions are _space-filling designs_
that cover the space “almost” uniformly. Trivially, equally-spaced points (or
uniform grids) satisfy the quasi-uniformity with $C_{\rm qu}=1$. Wenzel et
al., (2021) showed that points chosen sequentially to minimise GP posterior
variance for a Sobolev kernel are quasi-uniform. We refer to Wynne et al.,
(2021, p. 6) for further examples and a discussion on quasi-uniformity.

###  3.3 Hölder spaces

Section 4.1 studies the deterministic setting where the true unknown function
$f$ is assumed to belong to a Hölder space of functions. To define this space,
we first need the following definition.

######  Definition 3.

For $0<\alpha\leq 1$, a function $f:[0,T]\to\mathbb{R}$ is $\alpha$-Hölder
continuous if there exists a constant $L\geq 0$ such that, for all
$x,x^{\prime}\in[0,T]$,

$|f(x)-f(x^{\prime})|\leq L|x-x^{\prime}|^{\alpha}.$

Any such constant $L$ is called a _Hölder constant_ of $f$.

For $l\in\mathbb{N}\cup\\{0\\}$, denote by $C^{l}([0,T])$ the space of
functions $f\colon[0,T]\to\mathbb{R}$ such that the $l$th derivative $f^{(l)}$
exists and is continuous. For $l=0$, this is the space of continuous
functions. Hölder spaces are now defined as follows.

######  Definition 4.

Let $l\in\mathbb{N}\cup\\{0\\}$ and $0<\alpha\leq 1$. The Hölder space
$C^{l,\alpha}([0,T])$ consists of functions $f\in C^{l}([0,T])$ whose $l$th
derivative $f^{(l)}$ is $\alpha$-Hölder continuous.

Intuitively, $l+\alpha$ represents the smoothness of least-smooth functions in
$C^{l,\alpha}([0,T])$. It is well known that a sample path of Brownian motion
is almost surely $\alpha$-Hölder continuous if and only if $\alpha<1/2$ (e.g.,
Mörters and Peres,, 2010, Corollary 1.20), and thus it belongs to the Hölder
space $C^{l,\alpha}([0,T])$ with $l=0$ and $\alpha=1/2-\varepsilon$ almost
surely for arbitrarily small $\varepsilon>0$; in this sense, the smoothness of
a Brownian motion is $1/2$. As such, as is well known (e.g., Mörters and
Peres,, 2010, Theorem 1.27), a Brownian motion is almost nowhere
differentiable almost surely.

Note that we have the following strict inclusions:111These inclusions follow
from the following facts: By the definition of Hölder continuity, an
$\alpha_{1}$-Hölder continuous function is $\alpha_{2}$-Hölder continuous if
$\alpha_{1}>\alpha_{2}$; continuously differentiable functions are
$\alpha$-Hölder continuous for any $0<\alpha\leq 1$; not all Lipschitz
functions are differentiable.

  * •

$C^{l_{1},\alpha_{1}}([0,T])\subsetneq C^{l_{2},\alpha_{2}}([0,T])$ if (a)
$l_{1}>l_{2}$ or (b) $l_{1}=l_{2}$ and $\alpha_{1}>\alpha_{2}$,

  * •

$C^{l+1}([0,T])\subsetneq C^{l,1}([0,T])$.

###  3.4 Fractional Brownian motion

Section 4.2 considers the random setting where $f$ is a fractional (or
integrated fractional) Brownian motion (see Mandelbrot, (e.g., 1982, Chapter
IX)). Examples of these processes can be seen in Figures 1, 2, 5 and 6.

A fractional Brownian motion on $[0,T]$ with Hurst parameter $0<H<1$ is a
Gaussian process whose kernel is given by

$k_{0,H}(x,x^{\prime})=\big{(}\,\lvert x\rvert^{2H}+\lvert
x^{\prime}\rvert^{2H}-\lvert x-x^{\prime}\rvert^{2H}\big{)}/2.$ (11)

Note that if $H=1/2$, this is the Brownian motion kernel:
$k_{0,1/2}(x,x^{\prime})=\min(x,x^{\prime})$. The Hurst parameter $H$
quantifies the smoothness of the fractional Brownian motion. If $f_{\rm
FBM}\sim\mathcal{GP}(0,k_{0,H})$ for $H\in(0,1)$, then $f_{\rm FBM}\in
C^{0,H-\varepsilon}([0,T])$ almost surely for arbitrarily small
$\varepsilon>0$ (e.g., Nourdin,, 2012, Proposition 1.6).222That $f_{\rm
FBM}\notin C^{0,H}([0,T])$ almost surely for $f_{\rm
FBM}\sim\mathcal{GP}(0,k_{0,H})$ with $H\in(0,1)$ is a straightforward
corollary of, for example, Theorem 3.2 in Wang, (2007).

An integrated Brownian motion with Hurst parameter $H$ is defined via the
integration of a fractional Brownian motion with the same Hurst parameter: if
$f_{\rm FBM}\sim\mathcal{GP}(0,k_{0,H})$, then

$f_{\text{iFBM}}(x)=\int_{0}^{x}f_{\text{FBM}}(z)\,\mathrm{d}z,\quad
x\in[0,T]$

is an integrated Brownian motion with Hurst parameter $H$. It is a zero-mean
GP with the kernel

$\begin{split}k_{1,H}(x,x^{\prime})={}&\int_{0}^{x}\int_{0}^{x^{\prime}}\big{(}\lvert
z\rvert^{2H}+\lvert z^{\prime}\rvert^{2H}-\lvert
z-z^{\prime}\rvert^{2H}\big{)}/2\,\mathrm{d}z\,\mathrm{d}z^{\prime}\\\
={}&\frac{1}{2(2H+1)}\bigg{(}x^{\prime}x^{2H+1}+x(x^{\prime})^{2H+1}\\\
&\hskip
56.9055pt-\frac{1}{2(H+1)}\big{[}x^{2H+2}+(x^{\prime})^{2H+2}-|x-x^{\prime}|^{2H+2}\big{]}\bigg{)}.\end{split}$
(12)

Because differentiating an integrated fractional Brownian motion $f_{\rm
iFBM}\sim\mathcal{GP}(0,k_{1,H})$ yields a fractional Brownian motion $f_{\rm
FBM}\sim\mathcal{GP}(0,k_{0,H})$, a sample path of the former satisfies
$f_{\rm iFBM}\in C^{1,H-\varepsilon}([0,T])$ almost surely for arbitrarily
small $\varepsilon>0$; therefore the smoothness of $f_{\rm iFBM}$ is $1+H$.

###  3.5 Functions of finite quadratic variation

Some of our asymptotic results use the notion of functions of finite quadratic
variation, defined below.

######  Definition 5.

For each $N\in\mathbb{N}$, let
$\mathcal{P}_{N}\coloneqq(x_{N,n})_{n=1}^{N}\subset[0,T]$, and suppose that
$\|\mathcal{P}_{N}\|\to 0$ as $N\to\infty$. Then a function
$f:[0,T]\to\mathbb{R}$ is defined to have finite quadratic variation with
respect to $\mathcal{P}\coloneqq(\mathcal{P}_{N})_{N\in\mathbb{N}}$, if the
limit

$V^{2}(f)\coloneqq\lim_{N\to\infty}\sum_{n={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{0}}}^{N-1}\big{[}f(x_{N,n+1})-f(x_{N,n})\big{]}^{2}$
(13)

exists and is finite. We write $V^{2}(f,\mathcal{P})$ when it is necessary to
indicate the sequence of partitions.

Quadratic variation is defined for a specific sequence of partitions
$(\mathcal{P}_{N})_{N\in\mathbb{N}}$ and may take different values for
different sequences of partitions (Mörters and Peres,, 2010, Remark 1.36). For
conditions that guarantee the invariance of quadratic variation on the
sequence of partitions, see, for instance, Cont and Bas, (2023). Note also
that the notion of quadratic variation differs from that of $p$-variation for
$p=2$, which is defined as the supremum over all possible sequences of
partitions whose mesh sizes tend to zero.

If $f\in C^{0,\alpha}([0,T])$ with $\alpha>1/2$ and
$\|\mathcal{P}_{N}\|=\mathcal{O}(N^{-1})$ as $N\to\infty$, then we have
$V^{2}(f)=0$, because in this case

$\displaystyle\sum_{n={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{0}}}^{N-1}\big{[}f(x_{N,n+1})-f(x_{N,n})\big{]}^{2}\leq
N{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{L}}^{2}\max_{n}(\Delta
x_{N,n})^{2\alpha}$
$\displaystyle=\mathcal{O}(N^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{1-2\alpha}}})\to
0$

as $N\to\infty$. Therefore, given the inclusion properties of Hölder spaces
(see Section 3.3), we arrive at the following standard proposition.

######  Proposition 6.

Suppose that the partitions $(\mathcal{P}_{N})_{N\in\mathbb{N}}$ are such that
$\|\mathcal{P}_{N}\|=\mathcal{O}(N^{-1})$. If $f\in C^{l,\alpha}([0,T])$ for
$l+\alpha>1/2$, then $V^{2}(f)=0$.

If the mesh size tends to zero faster than $1/\log N$, in that
$\|\mathcal{P}_{N}\|=o(1/\log N)$, then the quadratic variation of almost
every sample path of the Brownian motion on the interval $[0,T]$ equals $T$
(Dudley,, 1973). This is of course true for partitions which have the faster
decay $\|\mathcal{P}_{N}\|=\mathcal{O}(N^{-1})$.

##  4 Main results

This section presents our main results on the asymptotic properties of the CV
and ML estimators, $\hat{\sigma}_{\textup{CV}}^{2}$ and
$\hat{\sigma}_{\textup{ML}}^{2}$, for the scale parameter. Section 4.1
considers the deterministic setting where the true function $f$ is fixed and
assumed to belong to a Hölder space. Section 4.2 studies the random setting
where $f$ is an (integrated) fractional Brownian motion. In Section 4.3, we
use the insights obtained in the proofs for the deterministic and random
settings to propose a _interior cross-validation_ (ICV) estimator, and show
its asymptotic properties are an improvement on those of CV and ML estimators.

###  4.1 Deterministic setting

We present our main results for the deterministic case where the true function
$f$ is fixed and assumed to be in a Hölder space $C^{l,\alpha}([0,T])$.
Theorem 7 below provides asymptotic upper bounds on the CV estimator
$\hat{\sigma}_{\textup{CV}}^{2}$ for different values of the smoothness
parameters $l$ and $\alpha$ of the Hölder space.

######  Theorem 7 (Rate of CV decay in Hölder spaces).

Suppose that $f$ is an element of $C^{l,\alpha}([0,T])$, with $l\geq 0$ and
$0<\alpha\leq 1$, such that $f(0)=0$, and the interval partitions
$(\mathcal{P}_{N})_{N\in\mathbb{N}}$ have bounded mesh sizes
$\|\mathcal{P}_{N}\|=\mathcal{O}(N^{-1})$ as $N\to\infty$. Then

$\hat{\sigma}_{\textup{CV}}^{2}=\mathcal{O}\big{(}N^{1-\min\\{2(l+\alpha),3\\}}\big{)}=\begin{cases}\mathcal{O}\left(N^{1-2\alpha}\right)&\text{
if
}\quad{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{l=0,}}\\\
\mathcal{O}\left(N^{-1-2\alpha}\right)&\text{ if }\quad l=1\text{ and
}\alpha<1/2,\\\ \mathcal{O}\left(N^{-2}\right)&\text{ if }\quad l=1\text{ and
}\alpha\geq 1/2,\\\ \mathcal{O}\left(N^{-2}\right)&\text{ if }\quad l\geq
2.\end{cases}$ (14)

###### Proof.

See Section 8.2. ∎

Theorem 8 below is a corresponding result for the ML estimator
$\hat{\sigma}_{\textup{ML}}^{2}$. Note that a similar result has been obtained
by Karvonen et al., (2020, Proposition 4.5), where the function $f$ is assumed
to belong to a Sobolev space and the kernel is a Matérn-type kernel. Theorem 8
is a version of this result where $f$ is in a Hölder space and the kernel is
the Brownian motion kernel; we provide it for completeness and ease of
comparison.

######  Theorem 8 (Rate of ML decay in Hölder spaces).

Suppose that $f$ is a non-zero element of $C^{l,\alpha}([0,T])$, with $l\geq
0$ and $0<\alpha\leq 1$, such that $f(0)=0$, and the interval partitions
$(\mathcal{P}_{N})_{N\in\mathbb{N}}$ have bounded mesh sizes
$\|\mathcal{P}_{N}\|=\mathcal{O}(N^{-1})$ as $N\to\infty$. Then

$\hat{\sigma}_{\textup{ML}}^{2}=\mathcal{O}\big{(}N^{1-\min\\{2(l+\alpha),2\\}}\big{)}=\begin{cases}\mathcal{O}\left(N^{1-2\alpha}\right)&\text{
if
}\quad{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{l=0,}}\\\
\Theta\left(N^{-1}\right)&\text{ if }\quad l\geq 1.\end{cases}$ (15)

###### Proof.

See Section 8.2. The proof is similar to that of Theorem 7. ∎

Figure 3 summarises the rates of Theorems 7 and 8. When $l+\alpha\leq 1$ (or
$l=0$ and $\alpha\leq 1$), the rates of $\hat{\sigma}_{\textup{CV}}^{2}$ and
$\hat{\sigma}_{\textup{ML}}^{2}$ are $\mathcal{O}(N^{1-2\alpha})$, so both of
them may decay (or grow, for $l+\alpha<1/2$) adaptively to the smoothness
$l+\alpha$ of the function $f$. However, when $l+\alpha>1$, the situation is
different: the decay rate of $\hat{\sigma}_{\textup{ML}}^{2}$ is always
$\Theta(N^{-1})$ and thus insensitive to $\alpha$, while that of
$\hat{\sigma}_{\textup{CV}}^{2}$ is $\left(N^{-1-2\alpha}\right)$ for $l=1$
and $\alpha\in(0,1/2]$. Therefore the CV estimator may be adaptive to a
broader range of the smoothness $0<l+\alpha\leq 3/2$ of the function $f$ than
the ML estimator (whose range of adaptation is $0<l+\alpha\leq 1$).

Note that Theorems 7 and 8 provide asymptotic upper bounds (except for the
case $l\geq 1$ of Theorem 8) and may not be tight if the function $f$ is
smoother than “typical” functions in $C^{l,\alpha}([0,T])$.333For example, if
$f(x)=|x-1/2|$ with $T=1$, we have $f\in C^{0,1}([0,T])$, as $f$ is Lipschitz
continuous in this case. However, $f$ is almost everywhere infinitely
differentiable except at one point $x=1/2$, so it is, in this sense, much
smoother than “typical” functions in $C^{0,1}([0,T])$. In Section 4.2, we show
that the bounds are indeed tight in expectation if $f$ is a fractional (or
integrated fractional) Brownian motion with smoothness $l+\alpha$.

In the deterministic setting, a potential approach for obtaining a matching
lower bound could use the rate of decay of the Fourier coefficients as a
notion of smoothness, instead of the Hölder smoothness condition on the
function $f$. Certain self-similarity conditions based on the decay rate and
behaviour of Fourier coefficients are routinely used to study coverage of
Bayesian credible sets (e.g., Szabó et al.,, 2015; Hadji and Szabó,, 2021) as
they define classes of functions that cannot “deceive” parameter estimators.
Motivated by this, we attempted to adapt the argument in Sniekers and van der
Vaart, (2015, Section 4.2) and Sniekers and van der Vaart, (2020, Section 10)
to derive a matching lower bound under a self-similarity assumption on the
Fourier coefficients. However, the bounds obtained through this approach
proved sub-optimal in our setting. A different technique may therefore be
required.

Figure 3: Rates of decay for the ML, CV and ICV estimators from Theorems 7, 8
and 14. Observe that the CV estimator’s range of adaptation to the smoothness
$l+\alpha$ is wider than the ML estimator’s, and the ICV estimator’s range of
adaptation is wider than that for both the CV and ML estimators.

######  Remark 9.

The proof of Theorem 8 shows that for $l=1$ we have
$\hat{\sigma}_{\textup{ML}}^{2}=\Theta(N^{-1})$ whenever
$\|\mathcal{P}_{N}\|\to 0$ as $N\to\infty$. More precisely, it establishes
that

$N\hat{\sigma}_{\textup{ML}}^{2}\to\|f^{\prime}\|_{\mathcal{L}^{2}([0,T])}:=\int_{0}^{T}f^{\prime}(x)^{2}\,\mathrm{d}x\quad\text{
as }\quad N\to\infty.$

Note that the $\mathcal{L}^{2}([0,T])$ norm of $f^{\prime}$ in the right hand
side equals the norm of $f$ in the reproducing kernel Hilbert space of the
Brownian motion kernel (e.g., van der Vaart and van Zanten,, 2008, Section 10)
Therefore, this fact is consistent with a similar more general statement in
Karvonen et al., (2020, Proposition 3.1).

In addition to the above results, Theorem 10 below shows the limit of the CV
estimator $\hat{\sigma}_{\textup{CV}}^{2}$ if the true function $f$ is of
finite quadratic variation.

######  Theorem 10.

For each $N\in\mathbb{N}$, let $\mathcal{P}_{N}\subset[0,T]$ be the equally-
spaced partition of size $N$. Suppose that $f:[0,T]\to\mathbb{R}$ has finite
quadratic variation $V^{2}(f)$ with respect to
$(\mathcal{P}_{N})_{N\in\mathbb{N}}$, $f(0)=0$, and $f$ is continuous on the
boundary, i.e., $\lim_{x\to 0^{+}}f(x)=f(0)$ and $\lim_{x\to T^{-}}f(x)=f(T)$.
Moreover, suppose that the quadratic variation $V^{2}(f)$ remains the same for
all sequences of quasi-uniform partitions with constant
$C_{\mathrm{qu}}=2$.444In Appendix B, we discuss the relaxation of this
requirement. Then

$\lim_{N\to\infty}\hat{\sigma}_{\textup{CV}}^{2}=\frac{V^{2}(f)}{T}.$ (16)

###### Proof.

See Section 8.2. ∎

For the ML estimator $\hat{\sigma}_{\textup{ML}}^{2}$, it is straightforward
to obtain a similar result by using (10) and (28) in Section 8.1: Under the
same conditions as Theorem 10, we have

$\lim_{N\to\infty}\hat{\sigma}_{\textup{ML}}^{2}=\frac{V^{2}(f)}{T}.$ (17)

Theorem 10 and (17) are consistent with Theorems 7 and 8, which assume $f\in
C^{l,\alpha}([0,T])$ with $l+\alpha>1/2$ and imply
$\hat{\sigma}_{\textup{CV}}^{2}\to 0$ and $\hat{\sigma}_{\textup{ML}}^{2}\to
0$ as $N\to\infty$. As summarized in Proposition 6, we have $V(f)=0$ for $f\in
C^{l,\alpha}([0,T])$ with $l+\alpha>1/2$, so Theorem 10 and (17) imply that
$\hat{\sigma}_{\textup{CV}}^{2}\to 0$ and $\hat{\sigma}_{\textup{ML}}^{2}\to
0$ as $N\to\infty$.

When $f$ is a Brownian motion, in which case the Brownian motion prior is
well-specified, the smoothness of $f$ is $l+\alpha=1/2$, and the quadratic
variation $V(f)$ becomes a positive constant (Dudley,, 1973). Proposition 13
in the next subsection shows that this fact, Theorem 10, and (17) lead to the
consistency of the ML and CV estimators in the well-specified setting.

###  4.2 Random setting

In Section 4.1, we obtained asymptotic upper bounds on the CV and ML scale
estimators when the true function $f$ is a fixed function in a Hölder space.
This section shows that these asymptotic bounds are tight in expectation when
$f$ is a fractional (or integrated fractional) Brownian motion.

That is, we consider the asymptotics of the expectations
$\mathbb{E}\hat{\sigma}_{\textup{CV}}^{2}$ and
$\mathbb{E}\hat{\sigma}_{\textup{ML}}^{2}$ under the assumption that
$f\sim\mathcal{GP}(0,k_{l,H})$, where $k_{l,H}$ is the kernel of a fractional
Brownian motion (11) for $l=0$ or that of an integrated fractional Brownian
motion (12) for $l=1$, with $0<H<1$ being the Hurst parameter. Recall that
$f\sim\mathcal{GP}(0,k_{l,H})$ belongs to the Hölder space
$C^{l,H-\varepsilon}([0,T])$ almost surely for arbitrarily small
$\varepsilon>0$, so its smoothness is $l+H$. Figure 4 summarises the obtained
upper and lower rates, corroborating the upper rates in Figure 3.

Theorems 11 and 12 below establish the asymptotic upper and lower bounds for
the CV and ML estimators, respectively.

######  Theorem 11 (Expected CV rate for fractional Brownian motion).

Suppose that $(\mathcal{P}_{N})_{N\in\mathbb{N}}$ are quasi-uniform and
$f\sim\mathcal{GP}(0,k_{l,H})$ with $l\in\\{0,1\\}$ and $0<H<1$. Then

$\mathbb{E}\hat{\sigma}_{\textup{CV}}^{2}=\Theta(N^{1-\min\\{2(l+H),3\\}})=\begin{cases}\Theta\left(N^{1-2H}\right)&\text{
if }\quad l=0\text{ and }H\in(0,1),\\\ \Theta\left(N^{-1-2H}\right)&\text{ if
}\quad l=1\text{ and }H<1/2,\\\ \Theta\left(N^{-2}\right)&\text{ if }\quad
l=1\text{ and }H\geq 1/2.\\\ \end{cases}$

###### Proof.

See Section 8.3. ∎

######  Theorem 12 (Expected ML rate for fractional Brownian motion).

Suppose that $(\mathcal{P}_{N})_{N\in\mathbb{N}}$ are quasi-uniform and
$f\sim\mathcal{GP}(0,k_{l,H})$ with $l\in\\{0,1\\}$ and $0<H<1$. Then

$\mathbb{E}\hat{\sigma}_{\textup{ML}}^{2}=\Theta(N^{1-\min\\{2(l+H),2\\}})=\begin{cases}\Theta\left(N^{1-2H}\right)&\text{
if }\quad l=0\text{ and }H\in(0,1),\\\ \Theta\left(N^{-1}\right)&\text{ if
}\quad l=1\text{ and }H\in(0,1).\end{cases}$

###### Proof.

See Section 8.3. The proof is similar to that of Theorem 11. ∎

Theorems 11 and 12 show that the CV estimator is adaptive to the unknown
smoothness $l+H$ of the function $f$ for a broader range $0<l+H\leq 3/2$ than
the ML estimator, whose range of adaptation is $0<l+H\leq 1$. These results
imply that the CV estimator can be asymptotically well-calibrated for a
broader range of unknown smoothness than the ML estimator, as discussed in
Section 5.

When the smoothness of $f$ is less than $1/2$, i.e., when $l+H<1/2$, the
Brownian motion prior, whose smoothness is $1/2$, is smoother than $f$. In
this case, the expected rates of $\hat{\sigma}_{\textup{CV}}^{2}$ and
$\hat{\sigma}_{\textup{ML}}^{2}$ are $\Theta\left(N^{1-2H}\right)$ and
increase as $N$ increases. The increase of $\hat{\sigma}_{\textup{CV}}^{2}$
and $\hat{\sigma}_{\textup{ML}}^{2}$ can be interpreted as compensating the
overconfidence of the posterior standard deviation
$\sqrt{\smash[b]{k_{N}(x)}}$, which decays too fast to be asymptotically well-
calibrated. This interpretation agrees with the illustration in Figure 1.

On the other hand, when $l+H>1/2$, the function $f$ is smoother than the
Brownian motion prior. In this case, $\hat{\sigma}_{\textup{CV}}^{2}$ and
$\hat{\sigma}_{\textup{ML}}^{2}$ decrease as $N$ increases, compensating the
under-confidence of the posterior standard deviation
$\sqrt{\smash[b]{k_{N}(x)}}$. See Figure 2 for an illustration.

When $l+H=1/2$, this is the well-specified case in that the smoothness of $f$
matches the Brownian motion prior. In this case, Theorems 11 and 12 yield
$\mathbb{E}\hat{\sigma}_{\textup{CV}}^{2}=\Theta(1)$ and
$\mathbb{E}\hat{\sigma}_{\textup{ML}}^{2}=\Theta(1)$, i.e., when the CV and ML
estimators converge, they converge to a positive constant. The following
proposition, which follows from Theorem 10 and (17), shows that this limiting
constant is the true value of the scale parameter $\sigma_{0}^{2}$ in the
well-specified setting $f\sim\mathcal{GP}(0,\sigma_{0}^{2}k)$, recovering
similar results in the literature (e.g., Bachoc et al.,, 2017, Theorem 2).

######  Proposition 13.

Suppose that $f\sim\mathcal{GP}(0,\sigma_{0}^{2}k)$ for $\sigma_{0}>0$ and
that partitions $(\mathcal{P}_{N})_{N\in\mathbb{N}}$ are equally-spaced. Then

$\lim_{N\to\infty}\hat{\sigma}_{\textup{CV}}^{2}=\lim_{N\to\infty}\hat{\sigma}_{\textup{ML}}^{2}=\sigma_{0}^{2}\quad\text{
almost surely}.$

###### Proof.

Since the quadratic variation of almost all sample paths of the unscaled
(i.e., $\sigma_{0}=1$) Brownian motion on $[0,T]$ equals $T$ (Dudley,, 1973),
the claim follows from (16) and (17). ∎

In Section 5, we discuss the implications of the obtained asymptotic rates of
$\hat{\sigma}_{\textup{CV}}^{2}$ and $\hat{\sigma}_{\textup{ML}}^{2}$ on the
reliability of the resulting GP uncertainty estimates. But first, motivated by
the results in Theorem 7 and Theorem 11, we propose a modification to the
cross-validation procedure that may have better asymptotic properties than the
CV estimator.

Figure 4: Expected decay rates for the ML, CV and ICV estimators from Theorems
11, 12 and 15. Observe that the CV estimator’s range of adaptation to the
smoothness $l+H$ is wider than the ML estimator’s, and the ICV estimator’s
range of adaptation is wider than that for both the CV and ML estimators.

###  4.3 Interior cross-validation estimators

The proofs of Theorems 7 and 11 show that when $l=1$ and $\alpha\in(1/2,1]$,
the bound on $\hat{\sigma}_{\textup{CV}}^{2}$ is dominated by the bound on
what we call the boundary terms. These are the terms corresponding to $n=1$
and $n=N$ in (9); see also (30). That the boundary terms dominate is
unsurprising since prediction at boundary points is a more challenging task
than prediction at the interior. Motivated by this observation, we propose an
alternative estimation method called _interior cross validation_ (ICV) that
maximises

$\sum_{n=2}^{N-1}\log p(f(x_{n})\,|\,x_{n},\mathbf{x}_{\setminus
n},f(\mathbf{x}_{\setminus n}),\theta).$

The corresponding scale parameter estimator is

$\hat{\sigma}_{\textup{ICV}}^{2}=\frac{1}{N}\sum_{n=2}^{N-1}\frac{\left[f(x_{n})-m_{\setminus
n}(x_{n})\right]^{2}}{k_{\setminus n}(x_{n})}.$ (18)

With the boundary points removed, the estimator’s range of adaptation to the
smoothness of the true function is greater than that for the CV estimator, as
illustrated in Figure 3 for the deterministic setting and Figure 4 for the
random setting. We present formal results for the deterministic and the random
settings in the following theorems.

######  Theorem 14 (Rate of ICV decay in Hölder spaces).

Suppose that $f$ is an element of $C^{l,\alpha}([0,T])$, with $l\geq 0$ and
$0<\alpha\leq 1$, such that $f(0)=0$, and the interval partitions
$(\mathcal{P}_{N})_{N\in\mathbb{N}}$ have bounded mesh sizes
$\|\mathcal{P}_{N}\|=\mathcal{O}(N^{-1})$ as $N\to\infty$. Then

$\hat{\sigma}_{\textup{ICV}}^{2}=\mathcal{O}\big{(}N^{1-\min\\{2(l+\alpha),4\\}}\big{)}=\begin{cases}\mathcal{O}\left(N^{1-2\alpha}\right)&\text{
if }\quad l=0,\\\ \mathcal{O}\left(N^{-1-2\alpha}\right)&\text{ if }\quad
l=1,\\\ \mathcal{O}\left(N^{-3}\right)&\text{ if }\quad l\geq 2.\end{cases}$

###### Proof.

See Section 8.4. ∎

######  Theorem 15 (Expected ICV rate for fractional Brownian motion).

Suppose that $(\mathcal{P}_{N})_{N\in\mathbb{N}}$ are quasi-uniform and
$f\sim\mathcal{GP}(0,k_{l,H})$ with $l\in\\{0,1\\}$ and $0<H<1$. Then

$\mathbb{E}\hat{\sigma}_{\textup{ICV}}^{2}=\Theta(N^{1-\min\\{2(l+H),4\\}})=\begin{cases}\Theta\left(N^{1-2H}\right)&\text{
if }\quad l=0,\\\ \Theta\left(N^{-1-2H}\right)&\text{ if }\quad
l=1.\end{cases}$

###### Proof.

See Section 8.4. ∎

One may take this idea further. For the Brownian motion kernel, an estimator
that does not attempt to predict on points “close enough” to the boundary,

$\hat{\sigma}_{\textup{ICV}}^{2}[N_{0}]=\frac{1}{N}\sum_{n=N_{0}}^{N-N_{0}}\frac{\left[f(x_{n})-m_{\setminus
n}(x_{n})\right]^{2}}{k_{\setminus n}(x_{n})}$

for some fixed $N_{0}$, has the same range of adaptation as
$\hat{\sigma}_{\textup{ICV}}^{2}=\hat{\sigma}_{\textup{ICV}}^{2}[1]$, the
estimator that only ignores the points on the boundary. However, for smoother
kernels like fractional Brownian motion (iFBM) and the Matérn family,
$\hat{\sigma}_{\textup{ICV}}^{2}[N_{0}]$ may exhibit adaptation beyond the
level $l=2$. The number of boundary points $N_{0}$ to remove would likely
depend on the smoothness of the kernel. Investigating model-dependent cross-
validation estimators that discard a proportion of boundary points would be an
interesting direction for future work.

##  5 Consequences for credible intervals

This section discusses whether the estimated scale parameter, given by the CV
or ML estimator, leads to asymptotically well-calibrated credible intervals.
With the kernel $\hat{\sigma}^{2}k(x,x^{\prime})$, where
$\hat{\sigma}^{2}=\hat{\sigma}_{\textup{CV}}^{2}$ or
$\hat{\sigma}^{2}=\hat{\sigma}_{\textup{ML}}^{2}$, a GP credible interval at
$x\in[0,T]$ is given by

$[m_{N}(x)-\alpha\hat{\sigma}\sqrt{k_{N}(x)},\quad
m_{N}(x)+\alpha\hat{\sigma}\sqrt{k_{N}(x)}]$ (19)

where $\alpha>0$ is a constant (e.g., $\alpha\approx 1.96$ leads to the 95%
credible interval).

As discussed in Section 1, this credible interval (19) is asymptotically well-
calibrated, if it shrinks to $0$ at the same speed as the decay of the error
$|m_{N}(x)-f(x)|$ as $N$ increases, i.e., the ratio

$\frac{|f(x)-m_{N}(x)|}{\hat{\sigma}\sqrt{k_{N}(x)}}$ (20)

should neither diverge to infinity nor converge to $0$. If this ratio diverges
to infinity, the credible interval (19) is asymptotically overconfident, in
that (19) shrinks to $0$ faster than the actual error $|f(x)-m_{N}(x)|$. If
the ratio converges to $0$, the credible interval is asymptotically
underconfident, as it increasingly overestimates the actual error. Therefore,
the ratio (20) should ideally converge to a positive constant for the credible
interval (19) to be reliable.

For ease of analysis, we focus on the random setting in Section 4.2 where $f$
is a fractional (or integrated fractional) Brownian motion and where we
obtained asymptotic upper and lower bounds for
$\mathbb{E}\hat{\sigma}_{\textup{CV}}^{2}$ and
$\mathbb{E}\hat{\sigma}_{\textup{ML}}^{2}$. We study how the expectation of
the posterior variance $\mathbb{E}\hat{\sigma}^{2}k_{N}(x)$ scales with the
expected squared error $\mathbb{E}[f(x)-m_{N}(x)]^{2}$. Specifically, we
analyze their ratio for $\hat{\sigma}^{2}=\hat{\sigma}_{\textup{CV}}^{2}$ and
$\hat{\sigma}^{2}=\hat{\sigma}_{\textup{ML}}^{2}$:

$R_{\textup{CV}}^{\mathbb{E}}(x,N)\coloneqq\frac{\mathbb{E}[f(x)-m_{N}(x)]^{2}}{\mathbb{E}\hat{\sigma}_{\textup{CV}}^{2}k_{N}(x)}\quad\text{
and }\quad
R_{\textup{ML}}^{\mathbb{E}}(x,N)\coloneqq\frac{\mathbb{E}[f(x)-m_{N}(x)]^{2}}{\mathbb{E}\hat{\sigma}_{\textup{ML}}^{2}k_{N}(x)}.$
(21)

The ratio diverging to infinity (resp. converging to $0$) as $N\to\infty$
suggests that the credible interval (19) is asymptotically overconfident
(resp. underconfident) for a non-zero probability of the realisation of $f$.
Thus ideally, the ratio should converge to a positive constant.

Theorem 16 below establishes the asymptotic rates of the ratios in (21).

######  Theorem 16.

Suppose that $(\mathcal{P}_{N})_{N\in\mathbb{N}}$ are quasi-uniform and
$f\sim\mathcal{GP}(0,k_{l,H})$ for $l\in\\{0,1\\}$ and $0<H<1$. Then,

$\sup_{x\in[0,T]}R_{\textup{CV}}^{\mathbb{E}}(x,N)=\begin{cases}\Theta(1)&\text{
if }\quad l=0\text{ and }H\in(0,1),\\\ \Theta(1)&\text{ if }\quad l=1\text{
and }H\in(0,1/2),\\\ \Theta\left(N^{1-2H}\right)&\text{ if }\quad l=1\text{
and }H\in(1/2,1),\\\ \end{cases}$

and

$\sup_{x\in[0,T]}R_{\textup{ML}}^{\mathbb{E}}(x,N)=\begin{cases}\Theta(1)&\text{
if }\quad l=0\text{ and }H\in(0,1),\\\ \Theta\left(N^{-2H}\right)&\text{ if
}\quad l=1\text{ and }H\in(0,1).\end{cases}$

###### Proof.

See Section 8.5. ∎

We have the following observations from Theorem 16, which suggest an advantage
of the CV estimator over the ML estimator for uncertainty quantification:

  * •

The ratio for the CV estimator neither diverges to infinity nor decays to $0$
in the range $0<l+H<3/2$, which is broader than that of the ML estimator,
$0<l+H<1$. This observation suggests that the CV estimator can yield
asymptotically well-calibrated credible intervals for a broader range of the
unknown smoothness $l+H$ of the function $f$ than the ML estimator.

  * •

The ratio decays to $0$ for the CV estimator in the range $3/2<l+H<2$ and for
the ML estimator in the range $1<l+H<2$. Therefore, the ML estimator may yield
asymptotically underconfident credible intervals for a broader range of the
smoothness $l+H$ than the CV estimator.

Moreover, for the interior CV estimator introduced in Equation 18, it
immediately follows from the proof in Section 8.5 that

$\sup_{x\in[0,T]}R_{\textup{ICV}}^{\mathbb{E}}(x,N)=\begin{cases}\Theta(1)&\text{
if }\quad l=0\text{ and }H\in(0,1),\\\ \Theta(1)&\text{ if }\quad l=1\text{
and }H\in(0,1),\end{cases}$

which implies the ICV estimator can yield asymptotically well-calibrated
credible intervals for a broader range of the smoothness than either the CV or
the ML estimator.

##  6 Experiments

This section describes numerical experiments to substantiate the theoretical
results in Section 4. We define test functions in Section 6.1, show empirical
asymptotic results for the CV estimator in Section 6.2, and report comparisons
between the CV and ML estimators in Section 6.3.

To this end, for a continuous function $f$, define
$l[f]\in\mathbb{N}\cup\\{0\\}$ and $\alpha\in(0,1]$ as

$l[f]:=\sup\\{l\in\mathbb{N}\cup\\{0\\}:f\in
C^{l}([0,T])\\},\quad\alpha[f]:=\sup\\{\alpha\in(0,1]:f\in
C^{l[f],\alpha}([0,T])\\}.$ (22)

Then, for arbitrarily small $\varepsilon_{1}\in\mathbb{N}$ and
$\varepsilon_{2}>0$, we have

$f\in
C^{\max(l[f]-\varepsilon_{1},0),\alpha[f]-\varepsilon_{2}}([0,T])\quad\text{
and }\quad f\notin C^{l[f]+\varepsilon_{1},\alpha[f]+\varepsilon_{2}}([0,T]).$

In this sense, $l[f]$ and $\alpha[f]$ characterize the smoothness of $f$.

###  6.1 Test functions

We generate test functions $f:[0,1]\to\mathbb{R}$ as sample paths of
stochastic processes with varying degrees of smoothness, as defined below. The
left columns of Figures 5 and 6 show samples of these functions.

  * •

To generate nowhere differentiable test functions, we use the Brownian motion
(BM), the Ornstein–Uhlenbeck process (OU), and the fractional Brownian motion
(FBM555We use https://github.com/crflynn/fbm to sample from FBM.) which are
zero-mean GPs with kernels

$\displaystyle k_{\textup{BM}}(x,x^{\prime})$
$\displaystyle=\min(x,x^{\prime}),\quad
k_{\textup{OU}}(x,x^{\prime})=\big{(}e^{-\lambda\lvert
x-x^{\prime}\rvert}-e^{-\lambda(x+x^{\prime})}\big{)}/4,$ $\displaystyle
k_{\textup{FBM}}(x,x^{\prime})$ $\displaystyle=\big{(}\,\lvert
x\rvert^{2H}+\lvert x^{\prime}\rvert^{2H}-\lvert
x-x^{\prime}\rvert^{2H}\big{)}/2,$

where $\lambda>0$ and $0<H<1$ is the Hurst parameter (recall that the FBM $=$
BM if $H=1/2$). We set $\lambda=0.2$ in the experiments below. Almost all
samples $f$ from these processes satisfy $l[f]=0$. For BM and OU we have
$\alpha[f]=1/2$ and for FBM $\alpha[f]=H$ (see Section 3.4). It is well known
that the OU process with the kernel $k_{\textup{OU}}$ above satisfies the
stochastic differential equation

$\mathrm{d}f(t)=-\lambda
f(t)\mathrm{d}t+\sqrt{\frac{\lambda}{2}}\,\mathrm{d}B(t),$ (23)

where $B$ is the standard Brownian motion whose kernel is $k_{\rm BM}$.

  * •

To generate differentiable test functions, we use once (iFBM) and twice
(iiFBM) integrated fractional Brownian motions

$f_{\textup{iFBM}}(x)=\int_{0}^{x}f_{\textup{FBM}}(z)\,\mathrm{d}z\quad\text{
and }\quad
f_{\textup{iiFBM}}(x)=\int_{0}^{x}f_{\textup{iFBM}}(z)\,\mathrm{d}z,$

where $f_{\textup{FBM}}\sim\mathcal{GP}(0,k_{\rm FBM})$. See (12) for the iFBM
kernel. With $H$ the Hurst parameter of the original FBM, almost all samples
$f$ from the above processes satisfy $l[f]=1$ and $\alpha[f]=H$ (iFBM) or
$l[f]=2$ and $\alpha[f]=H$ (iiFBM).

  * •

We also consider a piecewise infinitely differentiable function $f(x)=\sin
10x+[x>x_{0}]$, where $x_{0}$ is randomly sampled from the uniform
distribution on $[0,1]$ and $[x>x_{0}]$ is $1$ if $x>x_{0}$ and $0$ otherwise.
This function is of finite quadratic variation with $V^{2}(f)=1$.

Denote $\hat{\sigma}^{2}=\lim_{N\to\infty}\hat{\sigma}_{\textup{CV}}^{2}$. For
the above test functions, with equally-spaced partitions, we expect the
following asymptotic behaviours for the CV estimator from Theorems 7, 10 and
11, Proposition 13, the definition of quadratic variation, and Equation (23):

BM ($l[f]=0$, $\alpha[f]=1/2$):
$\displaystyle\quad\quad\hat{\sigma}_{\textup{CV}}^{2}=\mathcal{O}(1)\quad$
$\displaystyle\text{ and }\quad\hat{\sigma}^{2}=1,$ OU ($l[f]=0$,
$\alpha[f]=1/2$):
$\displaystyle\quad\quad\hat{\sigma}_{\textup{CV}}^{2}=\mathcal{O}(1)\quad$
$\displaystyle\text{ and }\quad\hat{\sigma}^{2}=\lambda/2,$ FBM ($l[f]=0$,
$\alpha[f]=H$):
$\displaystyle\quad\quad\hat{\sigma}_{\textup{CV}}^{2}=\mathcal{O}(N^{1-2H})\quad$
$\displaystyle\text{ and }\quad\hat{\sigma}^{2}=0,$ iFBM ($l[f]=1$,
$\alpha[f]=H$):
$\displaystyle\quad\quad\hat{\sigma}_{\textup{CV}}^{2}=\mathcal{O}(N^{-1-2H})\quad$
$\displaystyle\text{ and }\quad\hat{\sigma}^{2}=0,$ iiFBM ($l[f]=2$,
$\alpha[f]=H$):
$\displaystyle\quad\quad\hat{\sigma}_{\textup{CV}}^{2}=\mathcal{O}(N^{-2})\quad$
$\displaystyle\text{ and }\quad\hat{\sigma}^{2}=0,$ $\displaystyle\sin
10x+[x>x_{0}]:$
$\displaystyle\quad\quad\hat{\sigma}_{\textup{CV}}^{2}=\mathcal{O}(1)\quad$
$\displaystyle\text{ and }\quad\hat{\sigma}^{2}=1.$

Note that the above rate for the iFBM holds for $0<H\leq 1/2$. The chosen
functions allow us to cover a range of $\alpha[f]$ and $l[f]$ relevant to the
varying rate of convergence in Theorems 7 and 11, as well as a range of
$V^{2}(f)$ relevant to the limit in Theorem 10,
$\lim_{N\to\infty}\hat{\sigma}_{\textup{CV}}^{2}=V^{2}(f)/T$.

###  6.2 Asymptotics of the CV estimator

Figure 5 shows the asymptotics of $\hat{\sigma}_{\textup{CV}}^{2}$, where each
row corresponds to one stochastic process generating test functions $f$; the
rows are displayed in the increasing order of smoothness as quantified by
$l[f]+\alpha[f]$. The estimates are obtained for equally-spaced partitions of
sizes $N=10,10^{2},\dots,10^{5}$. In each row, the left panel plots a single
sample of generated test functions $f$. The middle panel shows the mean and
confidence intervals (of two standard deviations) of
$\hat{\sigma}_{\textup{CV}}^{2}$ for 100 sample realisations of $f$ for each
sample size $N$. The right panel describes the convergence rate of
$\hat{\sigma}_{\textup{CV}}^{2}$ to its limit point
$\hat{\sigma}^{2}=\lim_{N\to\infty}\hat{\sigma}_{\textup{CV}}^{2}$ on the log
scale.

We have the following observations:

  * •

The first two rows (the FBM and OU) and the last (the piece-wise infinitely
differentiable function) confirm Theorem 10, which states the convergence
$\hat{\sigma}_{\textup{CV}}^{2}\to V^{2}(f)/T$ as $N\to\infty$. While Theorem
10 does not provide convergence rates, the rates in the first two rows appear
to be $N^{-1/2}$. In the last row the rate is $N^{-2}$.

  * •

The remaining rows show that the observed rates of
$\hat{\sigma}_{\textup{CV}}^{2}$ to $0$ are in complete agreement with the
rates predicted by Theorems 7 and 11. In particular, the rates are adaptive to
the smoothness $l[f]+\alpha[f]$ of the function if $l[f]+\alpha[f]\leq 3/2$,
as predicted.

###  6.3 Comparison of CV and ML estimators

Figure 6 shows the decay rates of $\hat{\sigma}_{\textup{CV}}^{2}$ and
$\hat{\sigma}^{2}_{\textup{ML}}$ to $0$ for test functions $f$ with $l[f]=1$,
under the same setting as for Figure 5. In this case, Theorems 8 and 12
predict that $\hat{\sigma}_{\textup{ML}}^{2}$ decays at the rate
$\Theta(N^{-1})$ regardless of the smoothness; this is confirmed in the right
column. In contrast, the middle column shows again that
$\hat{\sigma}_{\textup{CV}}^{2}$ decays with a rate that adapts to $l[f]$ and
$\alpha[f]$ as long as $l[f]+\alpha[f]\leq 3/2$, as predicted by Theorems 7
and 11. These results empirically support our theoretical finding that the CV
estimator is adaptive to the unknown smoothness $l[f]+\alpha[f]$ of a function
$f$ for a broader range of smoothness than the ML estimator.

Additionally, in Appendix D, we compare the asymptotics of the CV and ML
estimators when the underlying kernel is a Matérn kernel and the Sobolev
smoothness of the true functions differs from that of the kernel. Similarly to
the results presented in this section, we observe that the CV estimator
exhibits a larger range of adaptation than the ML estimator.

Figure 5: Asymptotics of CV estimators for functions of varying smoothness as
quantified by $l[f]$ and $\alpha[l]$ in (22). Runs on individual 100 samples
from $f$ are in gray, means and confidence intervals (of two standard
deviations) are in black. Figure 6: Asymptotics of CV estimator compared to
asymptotics of ML estimators, for once differentiable functions.

##  7 Conclusion and future work

We have analysed the asymptotics of the CV and ML estimators for the scale
parameter in GP interpolation with the Brownian motion kernel. As a novel
contribution, our analysis covers the misspecified case where the smoothness
of the true function $f$ is different from that of the samples from the GP
prior. Our main results in Theorems 7, 8, 11 and 12 indicate that both CV and
ML estimators can adapt to the unknown smoothness of $f$, but the range of
smoothness for which this adaptation happens is broader for the CV estimator.
Accordingly, the CV estimator can make GP uncertainty estimates asymptotically
well-calibrated for a wider range of smoothness than the ML estimator, as
indicated in Theorem 16. In this sense, the CV estimator has an advantage over
the ML estimator. The experiments provide supporting evidence for the
theoretical results.

Natural next steps include the following:

  * •

Supplement the asymptotic upper bounds in Theorems 7 and 8 of the
deterministic setting with matching lower bounds.

  * •

Extend the analyses (of both the deterministic and random settings) to more
generic finitely smooth kernels, higher dimensions, and a noisy setting.

The matching lower bounds, if obtained, would allow one to analyse the ratio
between the prediction error $|f(x)-m_{N}(x)|$ and the posterior standard
deviation $\hat{\sigma}\sqrt{\smash[b]{k_{N}(x)}}$ in the deterministic
setting, corresponding to the one in Section 5 for the random setting. Such an
analysis would need additional assumptions on the true function $f$, such as
the homogeneity of the smoothness of $f$ across the input space. It also
requires a sharp characterisation of the error $|f(x)-m_{N}(x)|$, which could
use super convergence results in Wendland, (2005, Section 11.5) and Schaback,
(2018). Most natural kernel classes for extension are Matérns and other
kernels whose RKHS are norm-equivalent to Sobolev spaces; we conduct initial
empirical analysis in Appendix D and observe results consistent with the main
results in this paper. To this end, it would be possible to adapt the
techniques used in Karvonen et al., (2020) for analyzing the ML estimator to
the CV estimator. In any case, one would need much more advanced techniques
than those used here.

##  8 Proofs

This section provides the proofs of the main results and other lengthy
computations. For $x_{0}=0$ and $x_{1},\dots,x_{N}\in[0,T]$, we will use the
following notation whenever it can improve the readability or highlight a
point:

$\displaystyle\Delta x_{n}$ $\displaystyle\coloneqq x_{n+1}-x_{n},\quad
n=0,1,\dots,N-1,$ $\displaystyle f_{n}$ $\displaystyle\coloneqq f(x_{n}),\quad
n=0,1,\dots,N.$ (24)

###  8.1 Explicit expressions for the CV and ML estimators

Let us define $x_{0}=0$ and use the convention $f(x_{0})=0$.  By a direct
computation it is straightforward to verify that the inverse of the Gram
matrix of the Brownian motion kernel $k(x,x^{\prime})=\min(x,x^{\prime})$ over
the points $0=x_{0}<x_{1}<x_{2}<\dots<x_{N}$ is the band matrix

$\displaystyle
k(\mathbf{x},\mathbf{x})^{-1}=\begin{bmatrix}x_{1}&x_{1}&x_{1}&\dots&x_{1}&x_{1}\\\
x_{1}&x_{2}&x_{2}&\dots&x_{2}&x_{2}\\\ x_{1}&x_{2}&x_{3}&\dots&x_{3}&x_{3}\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
x_{1}&x_{2}&x_{3}&\dots&x_{N-1}&x_{N-1}\\\
x_{1}&x_{2}&x_{3}&\dots&x_{N-1}&x_{N}\end{bmatrix}^{-1}=\begin{bmatrix}b_{1}&c_{1}&0&\dots&0&0\\\
c_{1}&b_{2}&c_{2}&\dots&0&0\\\ 0&c_{2}&b_{3}&\dots&0&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&0&0&\dots&b_{N-1}&c_{N-1}\\\
0&0&0&\dots&c_{N-1}&b_{N}\\\ \end{bmatrix},$

where

$\displaystyle b_{i}$
$\displaystyle=\frac{x_{i+1}-x_{i-1}}{(x_{i-1}-x_{i})(x_{i}-x_{i+1})}\quad\text{for}\quad
i\in\\{2,\dots,N-1\\},\qquad b_{N}$ $\displaystyle=-\frac{1}{x_{N-1}-x_{N}},$
$\displaystyle c_{i}$
$\displaystyle=\frac{1}{(x_{i}-x_{i+1})}\quad\text{for}\quad
i\in\\{1,\dots,N-1\\}.$

It follows that the posterior mean and covariance functions in (3) can be
expressed as

$m_{N}(x)=\begin{dcases}\frac{(x_{n}-x)f(x_{n-1})+(x-x_{n-1})f(x_{n})}{x_{n}-x_{n-1}}&\text{
if }x\in[x_{n-1},x_{n}]\text{ for some }1\leq n\leq N,\\\ f(x_{N})&\text{ if
}x\in[x_{N},T]\end{dcases}$ (25)

and

$k_{N}(x,x^{\prime})=\begin{dcases}\frac{(x_{n}-x^{\prime})(x-x_{n-1})}{x_{n}-x_{n-1}}&\text{
if }x_{n-1}\leq x\leq x^{\prime}\leq x_{n}\text{ for some }1\leq n\leq N,\\\
x-x_{N}&\text{ if }x_{N}\leq x\leq x^{\prime}\leq T,\\\ 0&\text{
otherwise}.\end{dcases}$ (26)

We omit the case $x^{\prime}\leq x$ for $k_{N}(x,x^{\prime})$ as this case is
obtained by the symmetry $k_{N}(x,x^{\prime})=k_{N}(x^{\prime},x)$.

Using these expressions, we have, for each $1\leq n<N$:

$m_{\setminus
n}(x_{n})=\frac{(x_{n}-x_{n+1})f(x_{n-1})+(x_{n-1}-x_{n})f(x_{n+1})}{x_{n-1}-x_{n+1}}$

and

$k_{\setminus n}(x_{n})=k_{\setminus
n}(x_{n},x_{n})=\frac{(x_{n}-x_{n+1})(x_{n}-x_{n-1})}{x_{n-1}-x_{n+1}}$

For $n=N$, we have $m_{\setminus N}(x_{N})=f(x_{N-1})$ and $k_{\setminus
N}(x_{N})=x_{N}-x_{N-1}$. Inserting these expressions in (9) and using the
notation (24), the CV estimator can be written as

$\begin{split}\hat{\sigma}_{\textup{CV}}^{2}&=\frac{1}{N}\Bigg{[}\frac{(x_{2}f_{1}-x_{1}f_{2})^{2}}{x_{1}x_{2}\Delta
x_{1}}+\sum_{n=2}^{N-1}\frac{(\Delta x_{n-1}[f_{n+1}-f_{n}]-\Delta
x_{n}[f_{n}-f_{n-1}])^{2}}{(\Delta x_{n}+\Delta x_{n-1})\Delta x_{n}\Delta
x_{n-1}}\\\ &\hskip 34.14322pt+\frac{(f_{N}-f_{N-1})^{2}}{\Delta
x_{N-1}}\Bigg{]}.\end{split}$ (27)

For the ML estimator (8), we obtain the explicit expression

$\hat{\sigma}_{\textup{ML}}^{2}=\frac{1}{N}\sum_{n=1}^{N}\frac{[f(x_{n})-f(x_{n-1})]^{2}}{\Delta
x_{n-1}}$ (28)

by observing that $m_{n-1}(x_{n})=f(x_{n})$ and
$k_{n-1}(x_{n})=x_{n}-x_{n-1}$.

######  Remark 17.

The leave-$p$-out estimator $\hat{\sigma}_{\textup{CV}(p)}^{2}$ can be
expressed in a form similar (albeit more complicated) to Equation 27. We
derive this expression in Appendix C. This suggests that the analysis in
Section 4 could potentially be generalised to apply to the leave-$p$-out
estimators, a possibility that we leave open for future research to explore.

###  8.2 Proofs for Section 4.1

###### Proof of Theorem 7.

The estimator $\hat{\sigma}_{\textup{CV}}^{2}$ in (27) may be written as

$\hat{\sigma}_{\textup{CV}}^{2}=B_{1,N}+I_{N}+B_{2,N}$ (29)

in terms of the boundary terms

$B_{1,N}=\frac{1}{N}\cdot\frac{(x_{2}f_{1}-x_{1}f_{2})^{2}}{x_{1}x_{2}\Delta
x_{1}}\quad\text{ and }\quad
B_{2,N}=\frac{1}{N}\cdot\frac{(f_{N}-f_{N-1})^{2}}{\Delta x_{N-1}}$ (30)

and the interior term

$I_{N}=\frac{1}{N}\sum_{n=2}^{N-1}\frac{(\Delta x_{n-1}[f_{n+1}-f_{n}]-\Delta
x_{n}[f_{n}-f_{n-1}])^{2}}{(\Delta x_{n}+\Delta x_{n-1})\Delta x_{n}\Delta
x_{n-1}}.$ (31)

The claimed rate in (14) is $\mathcal{O}(N^{-2})$ if $l\geq 2$ or $l=1$ and
$\alpha\geq 1/2$. By the inclusion properties of Hölder spaces in Section 3.3,
it is therefore sufficient to consider the cases (a) $l=0$ and (b) $l=1$ and
$\alpha\in(0,1/2]$.

Suppose first that $l=0$. Let $L$ be a Hölder constant of a function $f\in
C^{0,\alpha}([0,T])$. Using the Hölder condition, the bounding assumption on
$\Delta x_{n}$, and $f_{0}=f(0)=0$, the boundary terms can be bounded as

$\displaystyle B_{1,N}=\frac{1}{N}\cdot\frac{(x_{1}(f_{1}-f_{2})+\Delta
x_{1}(f_{1}-f_{0}))^{2}}{x_{1}x_{2}\Delta x_{1}}$
$\displaystyle\leq\frac{1}{N}\cdot\frac{2(x_{1}^{2}(f_{1}-f_{2})^{2}+\Delta
x_{1}^{2}(f_{1}-f_{0})^{2})}{x_{1}x_{2}\Delta x_{1}}$
$\displaystyle\leq\frac{1}{N}\cdot\frac{2L^{2}(x_{1}^{2}\Delta
x_{1}^{2\alpha}+x_{1}^{2\alpha}\Delta x_{1}^{2})}{x_{1}x_{2}\Delta x_{1}}$
$\displaystyle=\mathcal{O}(N^{-1}\Delta x_{1}^{2\alpha-1})$
$\displaystyle=\mathcal{O}(N^{-2\alpha})$ (32)

and

$B_{2,N}=\frac{1}{N}\cdot\frac{(f_{N}-f_{N-1})^{2}}{\Delta
x_{N-1}}\leq\frac{1}{N}L^{2}\Delta
x_{N-1}^{2\alpha-1}=\mathcal{O}(N^{-2\alpha}).$ (33)

Similarly, the interior term is bounded as

$\displaystyle I_{N}$
$\displaystyle\leq\frac{2}{N}\sum_{n=2}^{N-1}\frac{\Delta
x_{n-1}^{2}(f_{n+1}-f_{n})^{2}+\Delta x_{n}^{2}(f_{n}-f_{n-1})^{2}}{(\Delta
x_{n}+\Delta x_{n-1})\Delta x_{n}\Delta x_{n-1}}$
$\displaystyle\leq\frac{2L^{2}}{N}\sum_{n=2}^{N-1}\frac{\Delta
x_{n-1}^{2}\Delta x_{n}^{2\alpha}+\Delta x_{n}^{2}\Delta
x_{n-1}^{2\alpha}}{(\Delta x_{n}+\Delta x_{n-1})\Delta x_{n}\Delta x_{n-1}}$
$\displaystyle=\frac{2L^{2}}{N}\sum_{n=2}^{N-1}\frac{\Delta x_{n-1}\Delta
x_{n}^{2\alpha-1}+\Delta x_{n}\Delta x_{n-1}^{2\alpha-1}}{\Delta x_{n}+\Delta
x_{n-1}}$ $\displaystyle=\frac{2L^{2}}{N}\sum_{n=2}^{N-1}\bigg{(}\frac{\Delta
x_{n-1}}{\Delta x_{n}+\Delta x_{n-1}}\Delta x_{n}^{2\alpha-1}+\frac{\Delta
x_{n}}{\Delta x_{n}+\Delta x_{n-1}}\Delta x_{n-1}^{2\alpha-1}\bigg{)}$
$\displaystyle\leq\frac{2L^{2}}{N}\sum_{n=2}^{N-1}\big{(}\Delta
x_{n}^{2\alpha-1}+\Delta x_{n-1}^{2\alpha-1}\big{)}$
$\displaystyle=\mathcal{O}(N^{1-2\alpha}).$

Inserting the above bounds in (29) yields
$\hat{\sigma}_{\textup{CV}}^{2}=\mathcal{O}(N^{-2\alpha}+N^{1-2\alpha})=\mathcal{O}(N^{1-2\alpha})$,
which is the claimed rate when $l=0$.

Suppose then that $l=1$ and $\alpha\in(0,1/2]$, so that the first derivative
$f^{\prime}$ of $f\in C^{1,\alpha}([0,T])$ is $\alpha$-Hölder and hence
continuous. Because a continuously differentiable function is Lipschitz, we
may set $\alpha=1$ in the estimates (32) and (33) for the boundary terms
$B_{1,N}$ and $B_{2,N}$ in the preceding case. This shows these terms are
$\mathcal{O}(N^{-2})$. Because $f$ is differentiable, we may use the mean
value theorem to write the interior term as

$\displaystyle I_{N}$ $\displaystyle=\frac{1}{N}\sum_{n=2}^{N-1}\frac{\Delta
x_{n-1}\Delta x_{n}}{\Delta x_{n-1}+\Delta
x_{n}}\bigg{(}\frac{f_{n+1}-f_{n}}{\Delta x_{n}}-\frac{f_{n}-f_{n-1}}{\Delta
x_{n-1}}\bigg{)}^{2}$ $\displaystyle=\frac{1}{N}\sum_{n=2}^{N-1}\frac{\Delta
x_{n-1}\Delta x_{n}}{\Delta x_{n-1}+\Delta
x_{n}}\big{[}f^{\prime}(\tilde{x}_{n})-f^{\prime}(\tilde{x}_{n-1})\big{]}^{2},$

where $\tilde{x}_{n}\in(x_{n},x_{n+1})$. Let $L$ be a Hölder constant of
$f^{\prime}$. Then the Hölder continuity of $f^{\prime}$ and the assumption
that $\|\mathcal{P}_{N}\|=\mathcal{O}(N^{-1})$ yield

$\displaystyle I_{N}\leq\frac{L^{2}}{N}\sum_{n=2}^{N-1}\frac{\Delta
x_{n-1}\Delta x_{n}}{\Delta x_{n-1}+\Delta
x_{n}}\lvert\tilde{x}_{n}-\tilde{x}_{n-1}\rvert^{2\alpha}$
$\displaystyle\leq\frac{L^{2}}{N}\sum_{n=2}^{N-1}\frac{\Delta x_{n-1}\Delta
x_{n}}{\Delta x_{n-1}+\Delta x_{n}}(\Delta x_{n-1}+\Delta x_{n})^{2\alpha}$
$\displaystyle\leq\frac{L^{2}}{N}\sum_{n=2}^{N-1}\Delta x_{n}(\Delta
x_{n-1}+\Delta x_{n})^{2\alpha}$ $\displaystyle=\mathcal{O}(N^{-2\alpha-1}).$

Using the above bounds in (29) yields
$\hat{\sigma}_{\textup{CV}}^{2}=\mathcal{O}(N^{-2}+N^{-2\alpha-1})=\mathcal{O}(N^{-2\alpha-1})$,
which is the claimed rate when $l=1$. ∎

###### Proof of Theorem 8.

From (28) we have

$\hat{\sigma}_{\textup{ML}}^{2}=\frac{1}{N}\sum_{n=1}^{N}\frac{(f_{n}-f_{n-1})^{2}}{\Delta
x_{n-1}}.$

Suppose first that $l=0$. As in the proof of Theorem 7, we get

$\hat{\sigma}_{\textup{ML}}^{2}=\frac{1}{N}\sum_{n=1}^{N}\frac{(f_{n}-f_{n-1})^{2}}{\Delta
x_{n-1}}\leq\frac{L^{2}}{N}\sum_{n=1}^{N}\Delta
x_{n-1}^{2\alpha-1}=\mathcal{O}\big{(}N^{1-2\alpha}\big{)}$ (34)

when $\|\mathcal{P}_{N}\|=\mathcal{O}(N^{-1})$. Suppose then that $l=1$. By
the mean value theorem there are $\xi_{n}\in(x_{n-1},x_{n})$ such that

$\hat{\sigma}_{\textup{ML}}^{2}=\frac{1}{N}\sum_{n=1}^{N}\frac{(f_{n}-f_{n-1})^{2}}{\Delta
x_{n-1}}=\frac{1}{N}\sum_{n=1}^{N}\Delta
x_{n-1}\bigg{(}\frac{f_{n}-f_{n-1}}{\Delta
x_{n-1}}\bigg{)}^{2}=\frac{1}{N}\sum_{n=1}^{N}\Delta
x_{n-1}f^{\prime}(\xi_{n})^{2}.$

Since $f^{\prime}$ is continuous on $[0,T]$ and hence Riemann integrable, we
obtain the asymptotic equivalence

$N\hat{\sigma}_{\textup{ML}}^{2}\to\int_{0}^{T}f^{\prime}(x)^{2}\,\mathrm{d}x\quad\text{
as }\quad N\to\infty$

when $\|\mathcal{P}_{N}\|\to 0$ as $N\to\infty$. The integral is positive
because $f$ has been assumed non-constant. ∎

###### Proof of Theorem 10.

For equally-spaced partitions, $\Delta x_{n}=x_{1}=T/N$ for all
$n\in\\{0,\dots,N-1\\}$, the estimator $\hat{\sigma}_{\textup{CV}}^{2}$ in
(27) takes the form

$\hat{\sigma}_{\textup{CV}}^{2}=\frac{1}{T}\Bigg{[}\frac{(x_{2}f_{1}-x_{1}f_{2})^{2}}{x_{1}x_{2}}+\frac{1}{2}\sum_{n=2}^{N-1}((f_{n+1}-f_{n})-(f_{n}-f_{n-1}))^{2}+(f_{N}-f_{N-1})^{2}\Bigg{]}.$

Recall from the proof of Theorem 7 the decomposition

$\hat{\sigma}_{\textup{CV}}^{2}=B_{1,N}+I_{N}+B_{2,N}$

in terms of the boundary terms $B_{1,N}$ and $B_{2,N}$ in (30) and the
interior term $I_{N}$ in (31). Because $f$ is assumed continuous on the
boundary and equispaced partitions are quasi-uniform, both $B_{1,N}$ and
$B_{2,N}$ tend to zero as $N\to\infty$. We may therefore focus on the interior
term, which decomposes as

$\begin{split}I_{N}={}&\frac{1}{2}\sum_{n=2}^{N-1}\left((f_{n+1}-f_{n})-(f_{n}-f_{n-1})\right)^{2}\\\
={}&\sum_{n=2}^{N-1}(f_{n+1}-f_{n})^{2}+(f_{n}-f_{n-1})^{2}-\frac{1}{2}(f_{n+1}-f_{n-1})^{2}\end{split}$

The sums $\sum_{n=2}^{N-1}(f_{n+1}-f_{n})^{2}$ and
$\sum_{n=2}^{N-1}(f_{n}-f_{n-1})^{2}$ tend to $V^{2}(f)$ by definition. To
establish the claimed bound we are therefore left to prove that

$\sum_{n=2}^{N-1}(f_{n+1}-f_{n-1})^{2}\to 2V^{2}(f)\qquad\text{as}\qquad
N\to\infty.$ (35)

We may write the sum as

$\displaystyle\sum_{n=2}^{N-1}(f_{n+1}-f_{n-1})^{2}$
$\displaystyle=\sum_{n=1}^{\lfloor\frac{N-1}{2}\rfloor}(f_{2n+1}-f_{2n-1})^{2}+\sum_{n=1}^{\lfloor\frac{N-2}{2}\rfloor}(f_{2n+2}-f_{2n})^{2}.$

Consider a sub-partition of $\mathcal{P}_{N}$ that consists of odd-index
points $x_{1},x_{3},\dots
x_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{2\lfloor\frac{N-1}{2}\rfloor+1}}}$
of $\mathcal{P}_{N}$. The sequence of these sub-partitions is quasi-uniform
with constant $2$. The assumption that the quadratic variation is $V^{2}(f)$
for all partitions with quasi-uniformity constant $2$ implies that

$\displaystyle\lim_{N\to\infty}\sum_{n=1}^{\lfloor\frac{N-1}{2}\rfloor}(f_{2n+1}-f_{2n-1})^{2}=V^{2}(f).$

The same will hold for sub-partitions formed of even-index points of
$\mathcal{P}_{N}$, giving

$\displaystyle\lim_{N\to\infty}\sum_{n=1}^{\lfloor\frac{N-2}{2}\rfloor}(f_{2n+2}-f_{2n})^{2}=V^{2}(f).$

Thus, (35) holds. This completes the proof. ∎

###  8.3 Proofs for Section 4.2

###### Proof of Theorem 11.

Recall the explicit expression of $\hat{\sigma}_{\textup{CV}}^{2}$ in (27):

$\begin{split}\hat{\sigma}_{\textup{CV}}^{2}&=\frac{1}{N}\Bigg{[}\frac{(x_{2}f_{1}-x_{1}f_{2})^{2}}{x_{1}x_{2}\Delta
x_{1}}+\sum_{n=2}^{N-1}\frac{(\Delta x_{n-1}[f_{n+1}-f_{n}]-\Delta
x_{n}[f_{n}-f_{n-1}])^{2}}{(\Delta x_{n}+\Delta x_{n-1})\Delta x_{n}\Delta
x_{n-1}}\\\ &\hskip 34.14322pt+\frac{(f_{N}-f_{N-1})^{2}}{\Delta
x_{N-1}}\Bigg{]}.\end{split}$ (36)

We consider the cases $l=0$ and $l=1$ separately. Recall that
$f\sim\mathcal{GP}(0,k_{l,H})$ implies that
$\mathbb{E}[f(x)f(x^{\prime})]=k_{l,H}(x,x^{\prime})$.

Suppose first that $l=0$, in which case $f\sim\mathcal{GP}(0,k_{0,H})$ for the
fractional Brownian motion kernel $k_{0,H}$ in (11). In this case the expected
values of squared terms in the expression for $\hat{\sigma}_{\textup{CV}}^{2}$
are $\mathbb{E}[x_{2}f_{1}-x_{1}f_{2}]^{2}=x_{1}x_{2}\Delta
x_{1}(x_{1}^{2H-1}+\Delta x_{1}^{2H-1}-(x_{1}+\Delta x_{1})^{2H-1})$,

$\displaystyle\mathbb{E}\big{[}$ $\displaystyle\Delta
x_{n-1}(f_{n+1}-f_{n})-\Delta x_{n}(f_{n}-f_{n-1})\big{]}^{2}$
$\displaystyle=\big{(}\Delta x_{n}^{2H-1}+\Delta x_{n-1}^{2H-1}-(\Delta
x_{n-1}+\Delta x_{n})^{2H-1}\big{)}\Delta x_{n-1}\Delta x_{n}(\Delta
x_{n}+\Delta x_{n-1}),$

and $\mathbb{E}[f_{N}-f_{N-1}]^{2}=\Delta x_{N-1}^{2H}$. Substituting these in
the expectation of $\hat{\sigma}_{\textup{CV}}^{2}$ and using the fact that
$\Delta x_{n}=\Theta(N^{-1})$ for all $n$ by quasi-uniformity we get

$\begin{split}\mathbb{E}\hat{\sigma}_{\textup{CV}}^{2}&=\frac{1}{N}\Bigg{[}(x_{1}^{2H-1}+\Delta
x_{1}^{2H-1}-(x_{1}+\Delta x_{1})^{2H-1})\\\ &\hskip
34.14322pt+\sum_{n=2}^{N-1}\left(\Delta x_{n-1}^{2H-1}+\Delta
x_{n}^{2H-1}-(\Delta x_{n-1}+\Delta x_{n})^{2H-1}\right)+\Delta
x_{N-1}^{2H-1}\Bigg{]}\\\ &=\frac{1}{N}\Bigg{[}\Delta
x_{1}^{2H-1}\left(\left(\frac{x_{1}}{\Delta
x_{1}}\right)^{2H-1}+1-\left(\frac{x_{1}}{\Delta
x_{1}^{2H-1}}+1\right)^{2H-1}\right)\\\ &\hskip 34.14322pt+\Delta
x_{n}^{2H-1}\sum_{n=2}^{N-1}\left(\left(\frac{\Delta x_{n-1}}{\Delta
x_{n}}\right)^{2H-1}+1-\left(\frac{\Delta x_{n-1}}{\Delta
x_{n}}+1\right)^{2H-1}\right)+\Delta x_{N-1}^{2H-1}\Bigg{]}\\\
&\eqqcolon{}\frac{1}{N}\Bigg{[}\Delta x_{1}^{2H-1}c_{1}+\Delta
x_{n}^{2H-1}\sum_{n=2}^{N-1}c_{n}+\Delta x_{N-1}^{2H-1}\Bigg{]}.\end{split}$

Notice that the function $x\mapsto x^{2H-1}+1-(x+1)^{2H-1}$ is positive for
$x>0$ and $H\in(0,1)$, and increasing for $H\in(1/2,1)$ and non-increasing for
$H\in(0,1/2]$. By quasi-uniformity we have $C_{\textup{qu}}^{-1}\leq\Delta
x_{n-1}/\Delta x_{n}\leq C_{\textup{qu}}$, and can bound $c_{n}$ for any $n$
and $N$ as

$\begin{split}0<C_{\textup{qu}}^{2H-1}+1-(C_{\textup{qu}}+1)^{2H-1}\leq
c_{n}\leq C_{\textup{qu}}^{1-2H}+1-(C_{\textup{qu}}^{-1}+1)^{2H-1}\>&\text{ if
}\>H\in(0,1/2],\\\
0<C_{\textup{qu}}^{1-2H}+1-(C_{\textup{qu}}^{-1}+1)^{2H-1}\leq c_{n}\leq
C_{\textup{qu}}^{2H-1}+1-(C_{\textup{qu}}+1)^{2H-1}\>&\text{ if
}\>H\in(1/2,1).\end{split}$

Finally, by quasi-uniformity $\Delta x_{n}=\Theta(N^{-1})$, and
$\mathbb{E}\hat{\sigma}_{\textup{CV}}^{2}=\Theta(N^{-2H})+\Theta(N^{1-2H})+\Theta(N^{-2H})=\Theta(N^{1-2H})$.

Suppose then that $l=1$, in which case $f\sim\mathcal{GP}(0,k_{1,H})$ for the
integrated fractional Brownian motion kernel $k_{1,H}$ in (12). It is
straightforward (though, in the case of the second expectation, somewhat
tedious) to compute that the expected values of squared terms in the
expression (36) for $\hat{\sigma}_{\textup{CV}}^{2}$ are

$\mathbb{E}[x_{2}f_{1}-x_{1}f_{2}]^{2}=\frac{x_{1}x_{2}\Delta
x_{1}}{2(H+1)(2H+1)}\big{(}x_{2}^{2H+1}-x_{1}^{2H+1}-\Delta
x_{1}^{2H+1}\big{)}$

and

$\begin{split}\mathbb{E}\big{[}\Delta&x_{n-1}(f_{n+1}-f_{n})-\Delta
x_{n}(f_{n}-f_{n-1})\big{]}^{2}\\\ ={}&\frac{\Delta x_{n}\Delta x_{n-1}(\Delta
x_{n}+\Delta x_{n-1})}{2(H+1)(2H+1)}\big{[}(\Delta x_{n}+\Delta
x_{n-1})^{2H+1}-\Delta x_{n}^{2H+1}-\Delta x_{n-1}^{2H+1}\big{]}\end{split}$
(37)

and

$\begin{split}\mathbb{E}[f_{N}-f_{N-1}]^{2}=\frac{\Delta
x_{N-1}}{2H+1}\bigg{(}x_{N}^{2H+1}-x_{N-1}^{2H+1}-\frac{1}{2(H+1)}\Delta
x_{N-1}^{2H+1}\bigg{)}.\end{split}$

Therefore, by (36),

$\begin{split}\mathbb{E}\hat{\sigma}_{\textup{CV}}^{2}={}&\frac{\big{(}x_{2}^{2H+1}-x_{1}^{2H+1}-\Delta
x_{1}^{2H+1}\big{)}}{2(H+1)(2H+1)N}\\\
&+\frac{1}{2(H+1)(2H+1)N}\sum_{n=2}^{N-1}\big{[}(\Delta x_{n}+\Delta
x_{n-1})^{2H+1}-\Delta x_{n}^{2H+1}-\Delta x_{n-1}^{2H+1}\big{]}\\\
&+\frac{1}{(2H+1)N}\bigg{(}x_{N}^{2H+1}-x_{N-1}^{2H+1}-\frac{1}{2(H+1)}\Delta
x_{N-1}^{2H+1}\bigg{)}\\\
\eqqcolon{}&\frac{1}{2(H+1)(2H+1)}B_{1,N}+\frac{1}{2(H+1)(2H+1)}I_{N}+\frac{1}{(2H+1)}B_{2,N}.\end{split}$

By quasi-uniformity, $B_{1,N}\leq N^{-1}x_{2}^{2H+1}=\mathcal{O}(N^{-2-2H})$.
Consider then the interior term

$\begin{split}I_{N}&=\frac{1}{N}\sum_{n=2}^{N-1}\Delta
x_{n}^{2H+1}\Bigg{[}\bigg{(}1+\frac{\Delta x_{n-1}}{\Delta
x_{n}}\bigg{)}^{2H+1}-\bigg{(}1+\bigg{(}\frac{\Delta x_{n-1}}{\Delta
x_{n}}\bigg{)}^{2H+1}\bigg{)}\Bigg{]}\\\
&\eqqcolon\frac{1}{N}\sum_{n=2}^{N-1}\Delta
x_{n}^{2H+1}c{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>(38)

Because the function $x\mapsto(1+x)^{2H+1}-(1+x^{2H+1})$ is positive and
increasing for $x>0$ if $H\in(0,1)$ and $C_{\textup{qu}}^{-1}\leq\Delta
x_{n-1}/\Delta x_{n}\leq C_{\textup{qu}}$ by quasi-uniformity, we have

$0<(1+C_{\textup{qu}}^{-1})^{2H+1}-(1+C_{\textup{qu}}^{-(2H+1)})\leq
c{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{{}^{\prime}}}_{n}\leq\bigg{(}1+\frac{\Delta
x_{n-1}}{\Delta x_{n}}\bigg{)}^{2H+1}\leq(1+C_{\textup{qu}})^{2H+1}$

for every $n$. Because $N^{-1}\sum_{n=2}^{N-1}\Delta
x_{n}^{2H+1}=\Theta(N^{-1-2H})$ by quasi-uniformity, we conclude from (38)
that $I_{N}=\Theta(N^{-1-2H})$. For the last term $B_{2,N}$, recall that we
have set $x_{N}=T$. Thus

$\begin{split}B_{2,N}&=\frac{1}{N}\bigg{(}T^{2H+1}-(T-\Delta
x_{N-1})^{2H+1}-\frac{1}{2(H+1)}\Delta x_{N-1}^{2H+1}\bigg{)}.\end{split}$

By the generalised binomial theorem,

$T^{2H+1}-(T-\Delta x_{N-1})^{2H+1}=(2H+1)T^{2H}\Delta
x_{N-1}+\mathcal{O}(\Delta x_{N-1}^{2})$

as $\Delta x_{N-1}\to 0$. It follows that under quasi-uniformity we have
$B_{2,N}=\Theta(N^{-2})$ for every $H\in(0,1)$. Putting these bounds for
$B_{1,N}$, $I_{N}$ and $B_{2,N}$ together we conclude that

$\begin{split}\mathbb{E}\hat{\sigma}_{\textup{CV}}^{2}&=\frac{1}{2(H+1)(2H+1)}B_{1,N}+\frac{1}{2(H+1)(2H+1)}I_{N}+\frac{1}{(2H+1)}B_{2,N}\\\
&=\mathcal{O}(N^{-2-2H})+\Theta(N^{-1-2H})+\Theta(N^{-2}),\end{split}$

which gives $\mathbb{E}\hat{\sigma}_{\textup{CV}}^{2}=\Theta(N^{-1-2H})$ if
$H\in(0,1/2]$ and $\mathbb{E}\hat{\sigma}_{\textup{CV}}^{2}=\Theta(N^{-2})$ if
$H\in[1/2,1)$. ∎

Observe that in the proof of Theorem 11 it is the boundary term $B_{2,N}$ that
determines the rate when there is sufficient smoothness, in that $l=1$ and
$H\in[1/2,1)$. Similar phenomenon occurs in the proof of Theorem 7. The
smoother a process is, the more correlation there is between its values at
far-away points. Because the Brownian motion (as well as fractional and
integrated Brownian motions) has a zero boundary condition at $x=0$ but no
boundary condition at $x=T$ and no information is available at points beyond
$T$, the importance of $B_{2,N}$ is caused by the fact that around $T$ one has
the least information about the process.

###### Proof of Theorem 12.

From (28) we get

$\mathbb{E}\hat{\sigma}_{\textup{ML}}^{2}=\frac{1}{N}\sum_{n=1}^{N}\frac{\mathbb{E}[f_{n}-f_{n-1}]^{2}}{\Delta
x_{n-1}}.$

We may then proceed as in the proof of Theorem 11 and use quasi-uniformity to
show that

$\mathbb{E}\hat{\sigma}_{\textup{ML}}^{2}=\frac{1}{N}\sum_{n=1}^{N}\frac{\mathbb{E}[f_{n}-f_{n-1}]^{2}}{\Delta
x_{n-1}}=\frac{1}{N}\sum_{n=1}^{N}\frac{\Delta x_{n-1}^{2H}}{\Delta
x_{n-1}}=\frac{1}{N}\sum_{n=1}^{N}\Delta x_{n-1}^{2H-1}=\Theta(N^{1-2H})$

when $l=0$ and

$\begin{split}\mathbb{E}\hat{\sigma}_{\textup{ML}}^{2}&=\sum_{n=1}^{N}\frac{\mathbb{E}[f_{n}-f_{n-1}]^{2}}{\Delta
x_{n-1}}\\\
&=\frac{1}{(2H+1)N}\sum_{n=1}^{N}\bigg{(}x_{n}^{2H+1}-x_{n-1}^{2H+1}-\frac{1}{2(H+1)}\Delta
x_{n-1}^{2H+1}\bigg{)}\\\
&=\frac{1}{(2H+1)N}\sum_{n=1}^{N}\bigg{(}(2H+1)x_{n}^{2H}\Delta
x_{n-1}+\mathcal{O}(\Delta x_{n-1}^{2})-\frac{1}{2(H+1)}\Delta
x_{n-1}^{2H+1}\bigg{)}\\\ &=\Theta(N^{-1})\end{split}$

when $l=1$. ∎

###  8.4 Proofs for Section 4.3

For the Brownian motion kernel, the ICV estimator defined in (18) takes the
explicit form

$\hat{\sigma}_{\textup{ICV}}^{2}=\frac{1}{N}\sum_{n=2}^{N-1}\frac{(\Delta
x_{n-1}[f_{n+1}-f_{n}]-\Delta x_{n}[f_{n}-f_{n-1}])^{2}}{(\Delta x_{n}+\Delta
x_{n-1})\Delta x_{n}\Delta x_{n-1}}.$

We analyse this estimator below.

###### Proof of Theorem 14.

The proof of Theorem 7 shows that when $l=1$ and $\alpha\in(1/2,1]$, the bound
is dominated by the bound on the boundary terms, $B_{1,N}=\mathcal{O}(N^{-2})$
and $B_{2,N}=\mathcal{O}(N^{-2})$, since

$\hat{\sigma}_{\textup{CV}}^{2}=B_{1,N}+I_{N}+B_{2,N}=\mathcal{O}(N^{-2})+\mathcal{O}(N^{-1-2\alpha})+\mathcal{O}(N^{-2})=\mathcal{O}(N^{-2}).$

As $\hat{\sigma}_{\textup{ICV}}^{2}=I_{N}$, it follows that
$\hat{\sigma}_{\textup{ICV}}^{2}=\mathcal{O}(N^{-1-2\alpha})$ when $l=1$. ∎

###### Proof of Theorem 15.

The proof of Theorem 11 shows that when $l=1$ and $H\in[1/2,1)$, the bound is
dominated by the bound on the right boundary terms, $B_{2,N}=\Theta(N^{-2})$,
since

$\begin{split}\mathbb{E}\hat{\sigma}_{\textup{CV}}^{2}&=\frac{1}{2(H+1)(2H+1)}B_{1,N}+\frac{1}{2(H+1)(2H+1)}I_{N}+\frac{1}{(2H+1)}B_{2,N}\\\
&=\mathcal{O}(N^{-2-2H})+\Theta(N^{-1-2H})+\Theta(N^{-2})\end{split}$

As $\mathbb{E}\hat{\sigma}_{\textup{ICV}}^{2}=I_{N}/(2(H+1)(2H+1))$, it
follows that $\hat{\sigma}_{\textup{ICV}}^{2}=\Theta(N^{-1-2H})$ when $l=1$. ∎

###  8.5 Proofs for Section 5

###### Proof of Theorem 16.

We only provide the proof for the case $l=1$ and leave the simpler case $l=0$
to the reader. Let $x\in(x_{n-1},x_{n})$. From the expression for $m_{N}$ in
Section 8.1, we get

$\begin{split}\mathbb{E}[f(x)-m_{N}(x)]^{2}&=\mathbb{E}\Bigg{[}f(x)-\frac{(x_{n}-x)f(x_{n-1})+(x-x_{n-1})f(x_{n})}{\Delta
x_{n-1}}\Bigg{]}^{2}\\\ &=\frac{1}{\Delta
x_{n-1}^{2}}\mathbb{E}\big{[}(x-x_{n-1})(f(x_{n})-f(x))-(x_{n}-x)(f(x)-f(x_{n-1}))\big{]}^{2}.\end{split}$

Then, we can use (37) with $x_{n}$ instead of $x_{n+1}$ and $x$ instead of
$x_{n}$ to get

$\mathbb{E}[f(x)-m_{N}(x)]^{2}=\frac{(x_{n}-x)(x-x_{n-1})}{C_{H}\Delta
x_{n-1}}\big{[}\Delta
x_{n-1}^{2H+1}-(x_{n}-x)^{2H+1}-(x-x_{n-1})^{2H+1}\big{]},$

where $C_{H}=2(H+1)(2H+1)$. The expression for $k_{N}$ in Section 8.1 gives

$\frac{\mathbb{E}[f(x)-m_{N}(x)]^{2}}{k_{N}(x)}=\frac{1}{C_{H}}\big{[}\Delta
x_{n-1}^{2H+1}-(x_{n}-x)^{2H+1}-(x-x_{n-1})^{2H+1}\big{]}.$

By removing the negative terms and using the quasi-uniformity (10), we obtain

$\sup_{x\in[0,T]}\frac{\mathbb{E}[f(x)-m_{N}(x)]^{2}}{k_{N}(x)}\leq\frac{(TC_{\textup{qu}})^{2H+1}}{C_{H}}N^{-1-2H},$

To see that this bound is tight, observe that for the midpoint
$x=(x_{n}+x_{n-1})/2$ we have $x_{n}-x=x-x_{n-1}=\Delta x_{n-1}/2$ and

$\frac{\mathbb{E}[f(x)-m_{N}(x)]^{2}}{k_{N}(x)}=\frac{1}{C_{H}}\bigg{(}1-\frac{1}{2^{2H}}\bigg{)}\Delta
x_{n-1}^{2H+1}\geq\frac{T^{2H+1}}{C_{H}C_{\textup{qu}}^{2H+1}}\bigg{(}1-\frac{1}{2^{2H}}\bigg{)}N^{-1-2H}$

by the quasi-uniformity. Therefore

$\sup_{x\in[0,T]}\frac{\mathbb{E}[f(x)-m_{N}(x)]^{2}}{k_{N}(x)}=\Theta(N^{-1-2H})$

when $l=1$. One can similarly show that

$\sup_{x\in[0,T]}\frac{\mathbb{E}[f(x)-m_{N}(x)]^{2}}{k_{N}(x)}=\Theta(N^{1-2H})$

when $l=0$. The claims then follow from the rates for
$\mathbb{E}\hat{\sigma}_{\textup{CV}}^{2}$ and
$\mathbb{E}\hat{\sigma}_{\textup{ML}}^{2}$ in Theorems 11 and 12. ∎

## Acknowledgements

The authors are grateful to the anonymous reviewers and the Associate Editor
for their time and effort in reviewing the manuscript, and the comments and
suggestions that improved the paper. In particular, one of the reviewers
suggested removal of the rightmost boundary point in the CV estimator, which
led to the introduction on the interior cross validation estimator. MN
acknowledges support from the U.K. Research and Innovation under grant number
EP/S021566/1. MK has been supported by the French government, through the 3IA
Cote d’Azur Investment in the Future Project managed by the National Research
Agency (ANR) with the reference number ANR-19-P3IA-0002. TK was supported by
the Academy of Finland postdoctoral researcher grant #338567 “Scalable,
adaptive and reliable probabilistic integration”. Part of this research was
carried out during a visit by TK to EURECOM in May 2023 that was funded by the
Institut français de Finlande, the Embassy of France to Finland, and the
Finnish Society of Sciences and Letters. MM gratefully acknowledges financial
support by the European Research Council through ERC StG Action 757275 /
PANAMA; the DFG Cluster of Excellence “Machine Learning - New Perspectives for
Science”, EXC 2064/1, project number 390727645; the German Federal Ministry of
Education and Research (BMBF) through the Tübingen AI Center (FKZ:
01IS18039A); and funds from the Ministry of Science, Research and Arts of the
State of Baden-Württemberg.

##  Appendix A Connection between the ML and CV estimators

Here we prove a connection between the ML and CV estimators; see Remark 1. Let

$C(N,p)=\binom{N}{p}=\frac{N!}{p!(N-p)!}$

denote the binomial coefficient. The leave-$p$-out cross-validation (LPO-CV)
estimator of $\sigma^{2}$ is

$\hat{\sigma}_{\textup{CV}(p)}^{2}=\frac{1}{C(N,p)}\sum_{i=1}^{C(N,p)}\frac{1}{p}\sum_{n=1}^{p}\frac{[f(x_{p,i,n})-m_{\setminus\\{p,i\\}}(x_{p,i,n})]^{2}}{k_{\setminus\\{p,i\\}}(x_{p,i,n})},$
(39)

where $i$ indexes the $N$-choose-$p$ possible sets of held-out datapoints,
$\mathbf{x}_{\setminus\\{p,i\\}}$, among $\mathbf{x}$ and $n\leq p$ the data
points left out of each of these sets. That is, for each $p$ and $i$ we have

$\mathbf{x}=\mathbf{x}_{\setminus\\{p,i\\}}\cup\\{x_{p,i,1},\ldots,x_{p,i,p}\\}.$

The functions $m_{\setminus\\{p,i\\}}$ and $k_{\setminus\\{p,i\\}}$ are the GP
conditional mean and variance based on the set
$\mathbf{x}_{\setminus\\{p,i\\}}$, which contains $N-p$ points. The purpose of
this section is to prove that

$\hat{\sigma}_{\textup{ML}}^{2}=\frac{1}{N}\sum_{p=1}^{N}\hat{\sigma}_{\textup{CV}(p)}^{2}.$
(40)

Denote
$\nu(\mathbf{x})=f(\mathbf{x})^{\top}k(\mathbf{x},\mathbf{x})^{-1}f(\mathbf{x})$.
The block matrix inversion formula applied to
$g(\mathbf{x}_{\setminus\\{p,i\\}}\cup\\{x\\})$ and the equations in Section 2
for the conditional mean and variance yield

$\frac{[f(x)-m_{\setminus\\{p,i\\}}(x)]^{2}}{k_{\setminus\\{p,i\\}}(x)}=\nu(\mathbf{x}_{\setminus\\{p,i\\}}\cup\\{x\\})-\nu(\mathbf{x}_{\setminus\\{p,i\\}})$
(41)

for any $1\leq p\leq N$ and $x\notin\mathbf{x}_{\\{p,i\\}}$, where we use the
convention $\nu(\mathbf{x}_{\setminus\\{N,i\\}})=\nu(\emptyset)=0$. For each
$1\leq p\leq N$, $i\leq C(N,p)$ and $n\leq p$ there is a unique index
$j(p,i,n)\leq C(N,p-1)$ such that

$\mathbf{x}_{\setminus\\{p,i\\}}\cup\\{x_{p,i,n}\\}=\mathbf{x}_{\setminus\\{p-1,j(p,i,n)\\}}.$
(42)

Setting $x=x_{p,i,n}$ in (41) gives

$\frac{[f(x_{p,i,n})-m_{\setminus\\{p,i\\}}(x_{p,i,n})]^{2}}{k_{\setminus\\{p,i\\}}(x_{p,i,n})}=\nu(\mathbf{x}_{\setminus\\{p,i\\}}\cup\\{x_{p,i,n}\\})-\nu(\mathbf{x}_{\setminus\\{p,i\\}}).$

Therefore

$\begin{split}\sum_{p=1}^{N}\hat{\sigma}_{\textup{CV}(p)}^{2}&=\frac{1}{N}\sum_{p=1}^{N}\frac{1}{C(N,p)}\sum_{i=1}^{C(N,p)}\frac{1}{p}\sum_{n=1}^{p}\frac{[f(x_{p,i,n})-m_{\setminus\\{p,i\\}}(x_{p,i,n})]^{2}}{k_{\setminus\\{p,i\\}}(x_{p,i,n})}\\\
&=\sum_{p=1}^{N}\frac{1}{C(N,p)}\sum_{i=1}^{C(N,p)}\frac{1}{p}\sum_{n=1}^{p}\big{[}\nu(\mathbf{x}_{\setminus\\{p,i\\}}\cup\\{x_{p,i,n}\\})-\nu(\mathbf{x}_{\setminus\\{p,i\\}})\big{]}.\end{split}$
(43)

By (42) from each set $\mathbf{x}_{\setminus\\{p,i\\}}$ on level $p$ (i.e.,
sets from which $p$ points have been left out) one can obtain $p$ sets on
level $p-1$ by adding one of the left-out datapoints. However, there are
$C(N,p)$ sets on level $p$ and $C(N,p-1)$ sets on level $p-1$. Hence for each
set $\mathbf{x}_{\setminus\\{p-1,j\\}}$ on level $p-1$ there are

$p\cdot\frac{C(N,p)}{C(N,p-1)}=p\cdot\frac{N!(p-1)!(N-p+1)!}{N!p!(N-p)!}=N-p+1$

combinations of sets $\mathbf{x}_{\setminus\\{p,i\\}}$ on level $p$ and points
$x_{p,i,n}$ left out of these sets such that
$\mathbf{x}_{\setminus\\{p,i\\}}\cup\\{x_{p,i,n}\\}=\mathbf{x}_{\setminus\\{p-1,j\\}}$.
Therefore

$\begin{split}\sum_{i=1}^{C(N,p)}&\frac{1}{p}\sum_{n=1}^{p}\big{[}\nu(\mathbf{x}_{\setminus\\{p,i\\}}\cup\\{x_{p,i,n}\\})-\nu(\mathbf{x}_{\setminus\\{p,i\\}})\big{]}\\\
&=\sum_{i=1}^{C(N,p)}\frac{1}{p}\sum_{n=1}^{p}\nu(\mathbf{x}_{\setminus\\{p,i\\}}\cup\\{x_{p,i,n}\\})-\sum_{i=1}^{C(N,p)}\frac{1}{p}\sum_{n=1}^{p}\nu(\mathbf{x}_{\setminus\\{p,i\\}})\\\
&=\frac{N-p+1}{p}\sum_{j=1}^{C(N,p-1)}\nu(\mathbf{x}_{\setminus\\{p-1,j\\}})-\sum_{i=1}^{C(N,p)}\nu(\mathbf{x}_{\setminus\\{p,i\\}})\end{split}$

and consequently (43) writes

$\begin{split}\sum_{p=1}^{N}\hat{\sigma}_{\textup{CV}(p)}^{2}&=\sum_{p=1}^{N}\frac{1}{C(N,p)}\Bigg{[}\frac{N-p+1}{p}\sum_{j=1}^{C(N,p-1)}\nu(\mathbf{x}_{\setminus\\{p-1,j\\}})-\sum_{i=1}^{C(N,p)}\nu(\mathbf{x}_{\setminus\\{p,i\\}})\Bigg{]}\\\
&=\sum_{p=1}^{N}\Bigg{[}\frac{1}{C(N,p-1)}\sum_{j=1}^{C(N,p-1)}\nu(\mathbf{x}_{\setminus\\{p-1,j\\}})-\frac{1}{C(N,p)}\sum_{i=1}^{C(N,p)}\nu(\mathbf{x}_{\setminus\\{p,i\\}})\Bigg{]},\end{split}$

which is a telescoping sum. We are left with

$\sum_{p=1}^{N}\hat{\sigma}_{\textup{CV}(p)}^{2}=\frac{1}{C(N,0)}\sum_{j=1}^{C(N,0)}\nu(\mathbf{x}_{\setminus\\{0,j\\}})-\frac{1}{C(N,N)}\sum_{i=1}^{C(N,N)}\nu(\mathbf{x}_{\setminus\\{N,i\\}}),$

where
$\nu(\mathbf{x}_{\setminus\\{0,j\\}})=f(\mathbf{x})^{\top}k(\mathbf{x},\mathbf{x})^{-1}f(\mathbf{x})$
and $\nu(\mathbf{x}_{\setminus\\{N,i\\}})=\nu(\emptyset)=0$. Thus

$\frac{1}{N}\sum_{p=1}^{N}\hat{\sigma}_{\textup{CV}(p)}^{2}=\frac{f(\mathbf{x})^{\top}k(\mathbf{x},\mathbf{x})^{-1}f(\mathbf{x})}{N}=\hat{\sigma}_{\textup{ML}}^{2},$

which establishes (40).

##  Appendix B Further discussion on Theorem 10

The requirement of having the same $V^{2}(f)$ for all sequences of partitions
quasi-uniform with constant $2$ can be relaxed somewhat: trivially, it is
sufficient that the quadratic variation is $V^{2}(f)$ specifically with
respect to even-points and odd-points sequences of sub-partitions used in the
proof in Section 8.2. Furthermore, we may even have different quadratic
variations with respect to said sequences. Then the results becomes

$\lim_{N\to\infty}\hat{\sigma}^{2}_{\mathrm{CV}}=\frac{\nu}{T}\qquad\text{for}\qquad\nu=\frac{V_{0}^{2}(f)+V_{1}^{2}(f)}{2},$

where $V_{0}^{2}(f)$ and $V_{1}^{2}(f)$ are quadratic variations with respect
to the even- and odd-points sub-partitions respectively, meaning that

$\displaystyle V^{2}(f)$
$\displaystyle=\lim_{N\to\infty}\sum_{n=1}^{N-1}(f_{n+1}-f_{n})^{2},$
$\displaystyle V_{0}^{2}(f)$
$\displaystyle=\lim_{N\to\infty}\sum_{n=1}^{\lfloor\frac{N-2}{2}\rfloor}(f_{2n+2}-f_{2n})^{2},$
$\displaystyle V_{1}^{2}(f)$
$\displaystyle=\lim_{N\to\infty}\sum_{n=1}^{\lfloor\frac{N-1}{2}\rfloor}(f_{2n+1}-f_{2n-1})^{2}.$

##  Appendix C Explicit expression for the leave-$p$-out estimator

Using the expressions for posterior mean and covariance functions in (25) and
(26), we may derive an explicit expression for the leave-$p$-out cross-
validation (LPO-CV) estimator of the scale parameter, given in (39) by

$\hat{\sigma}_{\textup{CV}(p)}^{2}=\frac{1}{C(N,p)}\sum_{i=1}^{C(N,p)}\frac{1}{p}\sum_{n=1}^{p}\frac{[f(x_{p,i,n})-m_{\setminus\\{p,i\\}}(x_{p,i,n})]^{2}}{k_{\setminus\\{p,i\\}}(x_{p,i,n})}.$

The expression is less straightforward than that for $p=1$. Denote by
$x_{\lfloor p,i,n\rfloor}$ the largest point in the set
$\mathbf{x}_{\setminus\\{p,i\\}}=\mathbf{x}\setminus\\{x_{p,i,1},\ldots,x_{p,i,p}\\}$
that does not exceed $x_{p,i,n}$, and by $x_{\lceil p,i,n\rceil}$ the smallest
point in the set $\mathbf{x}_{\setminus\\{p,i\\}}$ that exceeds $x_{p,i,n}$.
Through somewhat cumbersome arithmetic derivations one can show that the
estimator takes the form

$\displaystyle\hat{\sigma}_{\textup{CV}(p)}^{2}=\frac{1}{C(N,p)}\sum_{i=1}^{C(N,p)}\bigg{[}B_{p,i,1}+\sum_{n=2}^{p-1}I_{p,i,n}+B_{p,i,p}\bigg{]}$

where, for $\Delta x^{-}_{p,i,n}=(x_{p,i,n}-x_{\lfloor p,i,n\rfloor})$ and
$\Delta x^{+}_{p,i,n}=(x_{\lceil p,i,n\rceil}-x_{p,i,n})$, the inner term is

$I_{p,i,n}=\frac{\Delta x^{-}_{p,i,n}(f_{\lceil p,i,n\rceil}-f_{p,i,n})-\Delta
x^{+}_{p,i,n}(f_{p,i,n}-f_{\lfloor p,i,n\rfloor})}{(\Delta
x^{+}_{p,i,n}+\Delta x^{-}_{p,i,n})\Delta x^{+}_{p,i,n}\Delta x^{-}_{p,i,n}},$

and the boundary terms $B_{p,i,1}$ and $B_{p,i,p}$ depend on whether the $i$th
set contains $x_{1}$ or $x_{N}$, respectively. Specifically,

$\displaystyle B_{p,i,1}$ $\displaystyle=\begin{cases}\frac{(x_{\lceil
p,i,1\rceil}f_{p,i,1}-x_{p,i,1}f_{\lceil p,i,1\rceil})^{2}}{x_{p,i,1}x_{\lceil
p,i,1\rceil}\Delta x^{+}_{p,i,1}}&\text{if the $i$\textsuperscript{th} set
contains $x_{1}$,}\\\ I_{p,i,1}&\text{otherwise,}\end{cases}$ $\displaystyle
B_{p,i,p}$ $\displaystyle=\begin{cases}\frac{(f_{p,i,p}-f_{\lfloor
p,i,p\rfloor})^{2}}{\Delta x^{-}_{p,i,p}}&\text{if the $i$'th set contains
$x_{N}$,}\\\ I_{p,i,p}&\text{otherwise.}\end{cases}$

Though more cumbersome, it may be feasible to conduct convergence analysis
similar to that in Section 4 for $\hat{\sigma}_{\textup{CV}(p)}^{2}$. We leave
this up to future work.

##  Appendix D Comparison of CV and ML estimators for Matérn kernels

A natural next step is extending the analysis to kernels whose reproducing
kernel Hilbert spaces (RKHSs) are norm-equivalent to Sobolev spaces, such as
the commonly used Matérn kernels. The ML estimator for Matérn kernels was
analysed in Karvonen et al., (2020). Their experiments in Section 5.1 suggest
that, for $x_{+}\coloneqq\max(x,0)$,

$\hat{\sigma}_{\textup{ML}}^{2}=\Theta(N^{2(\nu_{\mathrm{model}}-2\nu_{\mathrm{true}})_{+}-1})$
(44)

when $k_{\nu_{\mathrm{model}}}$ is a Matérn kernel of order
$\nu_{\mathrm{model}}$ and $f$ is a finite linear combination of the form
$f=\sum_{i=1}^{m}\alpha_{i}k_{\nu_{\mathrm{true}}}(\cdot,x_{i})$ for some
$m\in\mathbb{N}$, $\alpha_{i}\in\mathbb{R}$, $z_{i}\in[0,1]$, and the Matérn
kernel $k_{\nu_{\mathrm{true}}}$ of order $\nu_{\mathrm{true}}$. Empirically,
we compare this to the rate of the CV estimator in Figure 8. The test
functions $f$ are posterior means of a GP with the $k_{\nu_{\mathrm{true}}}$
kernel conditioned on points $\\{(x_{1},y_{1}),\dots,(x_{10},y_{10})\\}$,
where each $x_{i}$ and $y_{i}$ is sampled i.i.d from the uniform distribution
on $[0,1]$. Since such $f$ are of the form
$f=\sum_{i=1}^{10}\alpha_{i}k_{\nu_{\mathrm{true}}}(\cdot,x_{i})$, we expect
the MLE rate in (44) to apply; we use experimental data and the results in
Theorems 7 and 12 to hypothesise what the rate in each individual example is.
Similarly to the observations for the Brownian motion kernel, we see that the
CV estimator adapts to the smoothness of the true function over a larger
ranger of smoothness compared to the ML estimator. For instance, for
$\nu_{\mathrm{model}}=1$, the experimental results suggest that the dependence
of rate on $\nu_{\mathrm{true}}$ is as illustrated in Figure 7. While the CV
and the ML estimators adapt to the function smoothness when
$\nu_{\mathrm{true}}\leq 1/2$, for $\nu_{\mathrm{true}}\in[1/2,3/4]$ only the
CV estimator continues adapting to the smoothness. This implies the CV
estimator is less likely to become asymptotically overconfident in the event
of undersmoothing.

Figure 7: Rates of decay for the ML and CV estimators for the Matérn kernel of
order $1$, and a true function that is a linear combination of Matérn kernels
of order $\nu_{\mathrm{true}}$. The ML rate is given in Karvonen et al.,
(2020, Equation 5.2). The CV rate is empirically observed in Figure 8. Observe
that the CV estimator’s range of adaptation to the smoothness
$\nu_{\mathrm{true}}$ is wider than the ML estimator’s.  Figure 8: Asymptotics
of CV estimator compared to asymptotics of the ML estimator, for the Matérn
kernel $\nu_{\mathrm{model}}$, and a true function that is a finite linear
combination of Matérn kernels $\nu_{\mathrm{true}}$.

## References

  * Anderes, (2010) Anderes, E. (2010).  On the consistent separation of scale and variance for Gaussian random fields.  The Annals of Statistics, 38(2):870–893.
  * Bachoc, (2013) Bachoc, F. (2013).  Cross validation and maximum likelihood estimations of hyper-parameters of Gaussian processes with model misspecification.  Computational Statistics & Data Analysis, 66:55–69.
  * Bachoc, (2018) Bachoc, F. (2018).  Asymptotic analysis of covariance parameter estimation for Gaussian processes in the misspecified case.  Bernoulli, 24(2):1531–1575.
  * Bachoc et al., (2020) Bachoc, F., Betancourt, J., Furrer, R., and Klein, T. (2020).  Asymptotic properties of the maximum likelihood and cross validation estimators for transformed Gaussian processes.  Electronic Journal of Statistics, 14(1):1962–2008.
  * Bachoc et al., (2017) Bachoc, F., Lagnoux, A., and Nguyen, T. M. N. (2017).  Cross-validation estimation of covariance parameters under fixed-domain asymptotics.  Journal of Multivariate Analysis, 160:42–67.
  * Beck and Guillas, (2016) Beck, J. and Guillas, S. (2016).  Sequential design with mutual information for computer experiments (MICE): Emulation of a tsunami model.  SIAM/ASA Journal on Uncertainty Quantification, 4(1):739–766.
  * Ben Salem et al., (2019) Ben Salem, M., Bachoc, F., Roustant, O., Gamboa, F., and Tomaso, L. (2019).  Gaussian process-based dimension reduction for goal-oriented sequential design.  SIAM/ASA Journal on Uncertainty Quantification, 7(4):1369–1397.
  * Bevilacqua et al., (2019) Bevilacqua, M., Faouzi, T., Furrer, R., and Porcu, E. (2019).  Estimation and prediction using generalized Wendland covariance functions under fixed domain asymptotics.  The Annals of Statistics, 47(2):828–856.
  * Chen et al., (2021) Chen, Y., Owhadi, H., and Stuart, A. M. (2021).  Consistency of empirical Bayes and kernel flow for hierarchical parameter estimation.  Mathematics of Computation, 90(332):2527–2578.
  * Cockayne et al., (2019) Cockayne, J., Oates, C. J., Sullivan, T. J., and Girolami, M. (2019).  Bayesian probabilistic numerical methods.  SIAM Review, 61(4):756–789.
  * Cont and Bas, (2023) Cont, R. and Bas, P. (2023).  Quadratic variation and quadratic roughness.  Bernoulli, 29(1):496–522.
  * Du et al., (2009) Du, J., Zhang, H., and Mandrekar, V. S. (2009).  Fixed-domain asymptotic properties of tapered maximum likelihood estimators.  The Annals of Statistics, 37(6A):3330–3361.
  * Dudley, (1973) Dudley, R. M. (1973).  Sample functions of the Gaussian process.  The Annals of Probability, 1(1):66–103.
  * Fong and Holmes, (2020) Fong, E. and Holmes, C. C. (2020).  On the marginal likelihood and cross-validation.  Biometrika, 107(2):489–496.
  * Garnett, (2023) Garnett, R. (2023).  Bayesian optimization.  Cambridge University Press.
  * Ginsbourger and Schärer, (2021) Ginsbourger, D. and Schärer, C. (2021).  Fast calculation of gaussian process multiple-fold cross-validation residuals and their covariances.  arXiv preprint arXiv:2101.03108.
  * Gu et al., (2018) Gu, M., Wang, X., and Berger, J. O. (2018).  Robust gaussian stochastic process emulation.  The Annals of Statistics, 46(6A):3038–3066.
  * Hadji and Szabó, (2021) Hadji, A. and Szabó, B. (2021).  Can we trust Bayesian uncertainty quantification from Gaussian process priors with squared exponential covariance kernel?  SIAM/ASA Journal on Uncertainty Quantification, 9(1):185–230.
  * Hennig et al., (2015) Hennig, P., Osborne, M. A., and Girolami, M. (2015).  Probabilistic numerics and uncertainty in computations.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471(2179):20150142.
  * Hennig et al., (2022) Hennig, P., Osborne, M. A., and Kersting, H. P. (2022).  Probabilistic Numerics: Computation as Machine Learning.  Cambridge University Press.
  * Jones et al., (1998) Jones, D. R., Schonlau, M., and Welch, W. J. (1998).  Efficient global optimization of expensive black-box functions.  Journal of Global Optimization, 13(4):455.
  * Karvonen, (2021) Karvonen, T. (2021).  Estimation of the scale parameter for a misspecified gaussian process model.  arXiv preprint arXiv:2110.02810.
  * Karvonen, (2023) Karvonen, T. (2023).  Asymptotic bounds for smoothness parameter estimates in Gaussian process interpolation.  SIAM/ASA Journal on Uncertainty Quantification, 11(4):1225–1257.
  * Karvonen and Oates, (2023) Karvonen, T. and Oates, C. J. (2023).  Maximum likelihood estimation in Gaussian process regression is ill-posed.  Journal of Machine Learning Research, 24(120):1–47.
  * Karvonen et al., (2020) Karvonen, T., Wynne, G., Tronarp, F., Oates, C. J., and Särkkä, S. (2020).  Maximum likelihood estimation and uncertainty quantification for Gaussian process approximation of deterministic functions.  SIAM/ASA Journal on Uncertainty Quantification, 8(3):926–958.
  * Kaufman and Shaby, (2013) Kaufman, C. G. and Shaby, B. A. (2013).  The role of the range parameter for estimation and prediction in geostatistics.  Biometrika, 100(2):473–484.
  * Kennedy and O’Hagan, (2001) Kennedy, M. C. and O’Hagan, A. (2001).  Bayesian calibration of computer models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63(3):425–464.
  * Loh, (2005) Loh, W.-L. (2005).  Fixed-domain asymptotics for a subclass of Matérn-type Gaussian random fields.  The Annals of Statistics, 33(5):2344–2394.
  * Loh and Kam, (2000) Loh, W.-L. and Kam, T.-K. (2000).  Estimating structured correlation matrices in smooth Gaussian random field models.  Annals of Statistics, 28(3):880–904.
  * Loh and Sun, (2023) Loh, W.-L. and Sun, S. (2023).  Estimating the parameters of some common Gaussian random fields with nugget under fixed-domain asymptotics.  Bernoulli, 29(3):2519–2534.
  * Loh et al., (2021) Loh, W.-L., Sun, S., and Wen, J. (2021).  On fixed-domain asymptotics, parameter estimation and isotropic Gaussian random fields with Matérn covariance functions.  The Annals of Statistics, 49(6):3127–3152.
  * Mandelbrot, (1982) Mandelbrot, B. B. (1982).  The Fractal Geometry of Nature.  WH Freeman New York.
  * Mörters and Peres, (2010) Mörters, P. and Peres, Y. (2010).  Brownian motion, volume 30.  Cambridge University Press.
  * Nourdin, (2012) Nourdin, I. (2012).  Selected aspects of fractional Brownian motion.  Number 4 in Bocconi & Springer Series. Springer.
  * O’Hagan, (1978) O’Hagan, A. (1978).  Curve fitting and optimal design for prediction.  Journal of the Royal Statistical Society. Series B (Methodological), 40(1):1–42.
  * O’Hagan, (2006) O’Hagan, A. (2006).  Bayesian analysis of computer code outputs: A tutorial.  Reliability Engineering & System Safety, 91(10-11):1290–1300.
  * Petit, (2023) Petit, S. (2023).  Maximum likelihood estimation and prediction error for a Matérn model on the circle.  arXiv:2209.07791v5.
  * Petit et al., (2022) Petit, S., Bect, J., Feliot, P., and Vazquez, E. (2022).  Parameter selection in Gaussian process interpolation: An empirical study of selection criteria.  arXiv:2107.06006v4.
  * Rasmussen and Williams, (2006) Rasmussen, C. E. and Williams, C. K. (2006).  Gaussian processes for machine learning, volume 2.  MIT press Cambridge, MA.
  * Sacks et al., (1989) Sacks, J., Welch, W. J., Mitchell, T. J., and Wynn, H. P. (1989).  Design and analysis of computer experiments.  Statistical science, 4(4):409–423.
  * Schaback, (2018) Schaback, R. (2018).  Superconvergence of kernel-based interpolation.  Journal of Approximation Theory, 235:1–19.
  * Shahriari et al., (2015) Shahriari, B., Swersky, K., Wang, Z., Adams, R. P., and De Freitas, N. (2015).  Taking the human out of the loop: A review of bayesian optimization.  Proceedings of the IEEE, 104(1):148–175.
  * Sniekers and van der Vaart, (2015) Sniekers, S. and van der Vaart, A. (2015).  Adaptive Bayesian credible sets in regression with a Gaussian process prior.  Electronic Journal of Statistics, 9(2):2475 – 2527.
  * Sniekers and van der Vaart, (2020) Sniekers, S. and van der Vaart, A. (2020).  Adaptive bayesian credible bands in regression with a gaussian process prior.  Sankhya A, 82(2):386–425.
  * Stein, (1990) Stein, M. L. (1990).  A comparison of generalized cross validation and modified maximum likelihood for estimating the parameters of a stochastic process.  The Annals of Statistics, 18(3):1139–1157.
  * Stein, (1993) Stein, M. L. (1993).  Spline smoothing with an estimated order parameter.  The Annals of Statistics, 21(3):1522–1544.
  * Stein, (1999) Stein, M. L. (1999).  Interpolation of Spatial Data: Some Theory for Kriging.  Springer Series in Statistics. Springer.
  * Sundararajan and Keerthi, (2001) Sundararajan, S. and Keerthi, S. S. (2001).  Predictive approaches for choosing hyperparameters in Gaussian processes.  Neural Computation, 13(5):1103–1118.
  * Szabó et al., (2015) Szabó, B., van der Vaart, A. W., and van Zanten, J. H. (2015).  Frequentist coverage of adaptive nonparametric Bayesian credible sets.  The Annals of Statistics, 43(4):1391–1428.
  * van der Vaart and van Zanten, (2008) van der Vaart, A. W. and van Zanten, J. H. (2008).  Reproducing Kernel Hilbert spaces of Gaussian Priors, volume 3 of IMS Collections, pages 200–222.  Institute of Mathematical Statistics.
  * Wang and Loh, (2011) Wang, D. and Loh, W.-L. (2011).  On fixed-domain asymptotics and covariance tapering in Gaussian random field models.  Electronic Journal of Statistics, 5:238–269.
  * Wang, (2007) Wang, W. (2007).  Almost-sure path properties of fractional Brownian sheet.  Annales de l’Institut Henri Poincare (B) Probability and Statistics, 43(5):619–631.
  * Wang, (2021) Wang, W. (2021).  On the inference of applying Gaussian process modeling to a deterministic function.  Electronic Journal of Statistics, 15(2):5014–5066.
  * Wendland, (2005) Wendland, H. (2005).  Scattered Data Approximation.  Number 17 in Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press.
  * Wenzel et al., (2021) Wenzel, T., Santin, G., and Haasdonk, B. (2021).  A novel class of stabilized greedy kernel approximation algorithms: Convergence, stability and uniform point distribution.  Journal of Approximation Theory, 262:105508.
  * Wynne et al., (2021) Wynne, G., Briol, F.-X., and Girolami, M. (2021).  Convergence guarantees for Gaussian process means with misspecified likelihoods and smoothness.  Journal of Machine Learning Research, 22(123):1–40.
  * Xu and Stein, (2017) Xu, W. and Stein, M. L. (2017).  Maximum likelihood estimation for a smooth Gaussian random field model.  SIAM/ASA Journal on Uncertainty Quantification, 5(1):138–175.
  * Ying, (1991) Ying, Z. (1991).  Asymptotic properties of a maximum likelihood estimator with data from a Gaussian process.  Journal of Multivariate Analysis, 36(2):280–296.
  * Ying, (1993) Ying, Z. (1993).  Maximum likelihood estimation of parameters under a spatial sampling scheme.  The Annals of Statistics, 21(3):1567–1590.
  * Zhang, (2004) Zhang, H. (2004).  Inconsistent estimation and asymptotically equal interpolations in model-based geostatistics.  Journal of the American Statistical Association, 99(465):250–261.