scale=1, angle=0, opacity=1, contents=

# Convergence rates of non-stationary and deep Gaussian process regression

Conor Moriarty-Osborne<EMAIL_ADDRESS>School of Mathematics and
Maxwell Institute for Mathematical Sciences, University of Edinburgh, King’s
Buildings, Edinburgh EH9 3FD, UK  Aretha L. Teckentrup<EMAIL_ADDRESS>School of Mathematics and Maxwell Institute for Mathematical Sciences,
University of Edinburgh, King’s Buildings, Edinburgh EH9 3FD, UK

###### Abstract

The focus of this work is the convergence of non-stationary and deep Gaussian
process regression. More precisely, we follow a Bayesian approach to
regression or interpolation, where the prior placed on the unknown function
$f$ is a non-stationary or deep Gaussian process, and we derive convergence
rates of the posterior mean to the true function $f$ in terms of the number of
observed training points. In some cases, we also show convergence of the
posterior variance to zero. The only assumption imposed on the function $f$ is
that it is an element of a certain reproducing kernel Hilbert space, which we
in particular cases show to be norm-equivalent to a Sobolev space. Our
analysis includes the case of estimated hyper-parameters in the covariance
kernels employed, both in an empirical Bayes’ setting and the particular
hierarchical setting constructed through deep Gaussian processes. We consider
the settings of noise-free or noisy observations on deterministic or random
training points. We establish general assumptions sufficient for the
convergence of deep Gaussian process regression, along with explicit examples
demonstrating the fulfilment of these assumptions. Specifically, our examples
require that the Hölder or Sobolev norms of the penultimate layer are bounded
almost surely.

AMS 2020 subject classifications: 62G08, 62G20, 60G17, 65C20, 68T07

##  1 Introduction

Given a model, such as a simulator of a real world physical process, we are
interested in forming an approximation, or emulator, of this model. Following
the Bayesian framework [27, 46], we construct a _prior distribution_ over the
space of potential models. This prior distribution encompasses our initial
knowledge of the model. Given some observations, which typically come in the
form of input-output pairs, we then condition the prior distribution on the
observations to construct a more informative probability distribution, known
as the _posterior distribution_. This posterior distribution then encompasses
our initial knowledge, as well as the observations. An advantage of the
Bayesian framework is that the posterior distribution allows for uncertainty
quantification in our approximation.

Mathematically, the model of interest corresponds to a function, mapping
inputs to corresponding outputs. A popular choice of prior distributions on
functions are Gaussian processes (GPs), see e.g. [43, 27, 20] for an
introduction to the topic. A GP is completely described by its mean and
covariance functions, and these functions can then be designed to describe
features we expect to see in our model, such as typical length scales and
amplitudes, or a prescribed regularity. In this work, we assume a zero mean in
the prior distribution, and focus on sophisticated covariance kernels to
capture the features of the model. The posterior mean, which is a linear
combination of kernel evaluations, then inherits the features from the
covariance kernel.

When the model being recovered is non-stationary, here meaning that it behaves
differently in different parts of its domain, standard choices of covariance
kernels such as the Matérn family [43, 28] fail to capture essential
properties of the model, and as a result, the accuracy of the emulator is
often poor. Many modern applications involve non-stationary models (see e.g.
[23, 32, 34]), and a considerable amount of work has been done over the last
decades to design non-stationary covariance kernels (see e.g. the recent
survey [48]). However, although the superior performance of these methods has
been observed in practice, convergence properties of these methods are often
still less understood than for classical, stationary constructions.

In the last decade, the use of deep GPs (DGPs) has further gained traction as
a method to create flexible priors, including capturing non-stationary
behaviour [11, 13, 49, 15]. DGPs are no longer GPs themselves, but are defined
through a hierarchy of stochastic processes that are conditionally Gaussian.
In particular, the stochastic process at layer $n$ of the hierarchy is used to
define the covariance structure of stochastic process at layer $n+1$, and
conditioned on the $n$-th layer, the distribution of the $n+1$-th layer is
Gaussian. DGPs can be seen as a particular example of a hierarchical
statistical model. Although DGPs have been shown to work much better than
standard stationary GPs in a wide variety of applications (see e.g. [31, 12,
52]), the theoretical understanding of DGPs, including convergence properties,
is still limited.

The focus of this work is to study convergence properties of regression (in
the case of noisy training data) or interpolation (in the case of noise-free
training data) with non-stationary and deep GP priors. More precisely, we
follow a Bayesian approach where the prior placed on the unknown function $f$
is a non-stationary or deep Gaussian process and the data comes in the form of
noisy or noise-free point evaluations of $f$, and we derive convergence rates
of the posterior mean to the true function $f$ in terms of the number of
observed training points. In some cases, we also show convergence of the
posterior variance to zero. Whether or not the training data is assumed noisy
depends on the context. If the training data consists of real data, the
assumption that the observed function values contain noise is common. If on
the other hand the training data is simulated data, obtained from simulating a
mathematical model on a computer, assuming no noise in the function values may
be more appropriate.

A large body of work exists exploring the convergence analysis of GP
regression (or interpolation), under various assumptions on the training data
and the true function $f$ generating the data, see e.g. the seminal works [8,
53, 59, 58, 51] as well as the more recent works [22, 26, 64, 36]. However,
the results known in the literature for GP regression are often either proven
for specific, stationary kernels, or they are proven for general kernels under
assumptions that have so far only been verified for stationary kernels. Hence,
there are still many open questions in the non-stationary case. For DGP
regression, the only work on convergence properties with the number of
training points that the authors are aware of is [15].

In this work, our primary interest is in the case where the training data is
noise-free, and we hence use proof techniques that give optimal convergence
rates in this case. The same approach was already used in for example [51, 26,
57], and uses tools mostly from the scattered data approximation literature
[62]. Although we do obtain convergence rates also in the case of noisy data,
the rates we obtain are not necessarily optimal. We consider the case where
the locations of the training data are fixed, given e.g. by a rule with known
theoretical properties, or sampled randomly from a specified distribution.

We further want to impose minimal assumptions on the true function $f$ being
approximated to guarantee convergence of the methodology. In particular, in
the vast majority of our results we only assume that $f$ belongs to the
Sobolev space $H^{\beta}(\Omega)$ for a bounded domain
$\Omega\subseteq\mathbb{R}^{d}$ and $\beta\geq d/2$. This is in contrast to
[15], where the authors construct a function space relating to the DGP prior
distribution and show convergence of the methodology when recovering a
function from this particular function space. This imposes structural
assumptions on $f$, which may not be straight-forward to verify in practice.

Finally, we note here that as in [15], our final results on the convergence of
DGP regression for specific constructions require a truncation of the
(conditionally) Gaussian process on the penultimate layer, such that
appropriate Hölder or Sobolev norms of the penultimate layer are uniformly
bounded almost surely. Although we do present general convergence results for
DGP regression that do not require this truncation, we were unable to verify
the assumptions of these results without enforcing the truncation.

###  1.1 Our contributions

In this paper we make the following novel contributions to the analysis of
non-stationary and deep GP regression, of which the first two are also of
independent interest:

  1. 1.

We show that for two typical non-stationary kernels, namely the warping and
mixture kernels introduced in section 2.2.2, the reproducing kernel Hilbert
space (RKHS) can be found explicitly, and we prove that under given conditions
the RKHS is norm-equivalent to a Sobolev space. These results are given in
Theorem 3.7 and 3.8.

  2. 2.

We present results on the regularity properties of sample paths of a non-
stationary Gaussian process for general kernels, as given in Lemmas 4.3 and
4.8. Additionally, specific applications to the warping and mixture kernels
are discussed in Corollaries 4.4 and 4.9.

  3. 3.

We prove convergence rates for non-stationary GP regression in terms of the
number of training points, under various assumptions on the true function $f$
being approximated and the training data used.

For noise-free data, these results are largely based on the known results for
general kernels given in Propositions 5.1 and 5.6, with the required
assumptions proven here for the warping and mixture kernels. The main final
results are given in Corollaries 5.5, 5.7, 5.11, 5.13, 5.18 and 5.19.
Corollaries 5.5 and 5.11 also apply to the convolution kernel from section
2.2.2.

For noisy data, results for general kernels are given in Theorem 5.8 (adapted
from [64, Theorem 2]) and [58, Theorem 2], and results for specific kernels
are given in Corollaries 5.9, 5.14 and 5.20.

  4. 4.

We prove convergence rates for deep GP regression in terms of the number of
training points used, under general assumptions on the true function $f$ being
approximated and the training data used. The results for general kernels and
constructions is given in Theorems 6.1 and 6.7.

Explicit constructions of (constrained) deep GPs based on the warping and
mixture kernels, and satisfying the assumptions of Theorem 6.1 or 6.7, are
provided in Lemmas 6.2, 6.3 and 6.8.

To the best of our knowledge, the error bounds proven in items 3 and 4 are the
first to apply in the settings considered in this work. In particular, the
analysis in section 6 is the first to show convergence of deep GP regression
under the general assumption that the true function $f$ belongs to a Sobolev
space, rather than imposing further structural assumptions (as in e.g. [15]).

###  1.2 Paper structure

This paper is structured as follows. In section 2, we introduce the general
set-up of GP regression, together with standard choices of stationary and non-
stationary covariance kernels and the estimation of related hyper-parameters.
In section 3, we introduce and construct native spaces for the non-stationary
covariance kernels, and show that these are norm-equivalent to Sobolev spaces
in special cases. Section 4 gives results on the sample path regulairty of
non-stationary Gaussian processes. Sections 5 and 6 are devoted to the error
analysis on non-stationary and deep GP regression, respectively. Finally, in
section 7, we present simple numerical simulations that illustrate the theory.
Section 8 offers some conclusions and discussion.

###  1.3 Preliminaries and notation

Throughout this work, $\Omega\subset\mathbb{R}^{d}$, for $d\in\mathbb{N}$,
will be a bounded Lipschitz domain that satisfies an interior cone condition
[62]. For any multi-index $q=(q_{1},\ldots,q_{n})\in\mathbb{N}^{d}_{0}$ with
length $|q|=q_{1}+\ldots+q_{n}$ we denote the $q$-th derivative of a function
$g:\Omega\rightarrow\mathbb{R}$ by

$D^{q}g=\frac{\partial^{|q|}}{\partial^{q_{1}}\ldots\partial^{q_{d}}}g.$

In the case $d=1$, we use the simpler notation $g^{\prime}:=D^{1}g$ and
$g^{(q)}:=D^{q}g$. When $g$ is defined on $\Omega\times\Omega$, the operator
$D^{q_{1},q_{2}}$ computes the $q_{1}$-th partial derivative with respect to
the first input and the $q_{2}$-th partial derivative with respect to the
second input, where $(q_{1},q_{2})\in\mathbb{N}_{0}^{2d}$. For $p\geq 0$ we
use the notation $g\in C^{p}(\Omega)$ if for all multi-indices $|q|\leq\lfloor
p\rfloor$ the derivative $D^{q}g$ exists and is $p-\lfloor p\rfloor$-Hölder
continuous. The norm on $C^{p}(\Omega)$ is given by

$\|g\|_{C^{p}(\Omega)}=\sum_{|\beta|\leq\lfloor
p\rfloor}\sup_{u\in\Omega}|D^{\beta}g(u)|+\sum_{|\beta|=\lfloor
p\rfloor}\sup_{u,u^{\prime}\in\Omega,u\neq
u^{\prime}}\frac{|D^{\beta}g(u)-D^{\beta}g(u^{\prime})|}{\|u-u^{\prime}\|_{2}^{p-\lfloor
p\rfloor}}.$

When $p$ is an integer this norm is the standard norm on continuously
differentiable functions. The space of smooth functions is defined as
$C^{\infty}(\Omega)=\cap_{p\geq 0}C^{p}(\Omega)$. Additionally, we use the
notation $k\in C^{p,p}(\Omega\times\Omega)$ for $p\geq 0$ when
$k(u_{1},\cdot)\in C^{p}(\Omega)$ and $k(\cdot,u_{2})\in C^{p}(\Omega)$ for
all $u_{1},u_{2}\in\Omega$. Note that $C^{2p}(\Omega\times\Omega)\subset
C^{p,p}(\Omega\times\Omega)$ for all $p\geq 0$.

For $\beta\in\mathbb{N}$, the Sobolev space of functions with square
integrable weak derivatives up to the $\beta$-th order is denoted
$H^{\beta}(\Omega)$. For $\beta\in\mathbb{R}_{+}\setminus\mathbb{N}$, with
$1\leq\beta<\infty$, we define $H^{\beta}(\Omega)$ as a fractional order
Sobolev space known as a Sobolev–Slobodeckij space (see, e.g. [5]), equipped
with the norm

$\|g\|_{H^{\beta}(\Omega)}:=\sqrt{\sum_{|q|\leq\lfloor\beta\rfloor}\int_{\Omega}|D^{q}g(u)|^{2}\mathrm{d}u}+\sup_{|q|=\lfloor\beta\rfloor}\sqrt{\int_{\Omega}\int_{\Omega}\frac{|D^{q}g(u)-D^{q}g(v)|^{2}}{|u-v|^{2(\beta-\lfloor\beta\rfloor)+1}}\mathrm{d}u\mathrm{d}v}.$

The special case $\beta=0$ is denoted by $L^{2}(\Omega):=H^{0}(\Omega)$. For
$\beta\in\mathbb{N}$, $W^{\beta,\infty}(\Omega)$ denotes the space of
functions with essentially bounded weak derivatives up to the $\beta$-th
order, equipped with the norm

$\|g\|_{W^{\beta,\infty}(\Omega)}=\sum_{|\beta|\leq\lfloor
p\rfloor}\operatorname*{ess\,sup}_{u\in\Omega}|D^{\beta}g(u)|.$

For vector spaces $V_{1},V_{2}$ the notation $V_{1}\hookrightarrow V_{2}$
denotes the continuous embedding of $V_{1}$ into $V_{2}$ and $V_{1}\cong
V_{2}$ denotes that $V_{1}$ and $V_{2}$ are equal as vector spaces and their
respective norms are equivalent. For any function $g:V_{1}\to V_{2}$ we denote
by $\operatorname{inv}(g):V_{2}\to V_{1}$ the inverse function of $g$, whereas
$1/g$ refers to its reciprocal function.

We have further summarised the notation used in this work in Table 1 .

Table 1: Notation $\Omega$ | Bounded Lipschitz domain satisfying an interior cone condition
---|---
$d$ | Dimension of domain $\Omega\subset\mathbb{R}^{d}$
$N$ | Number training points $U_{N}\subset\Omega$
$L$ | Number of components in mixture kernel $k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}}$
$D$ | Depth of a deep GP $\\{f^{n}\\}_{n=0}^{D-1}$
$\mathcal{GP}$ | A Gaussian process
$k_{\text{Mat}(\nu)}$ | Matérn kernel
$\sigma^{2}$ | Marginal variance of Matérn kernel
$\lambda$ | Correlation length of Matérn kernel
$\nu$ | Smoothness parameter of Matérn kernel
$k_{\operatorname{warp}}^{w,k_{s}}$ | Warping kernel, see (5)
$k_{\operatorname{mix}}^{{\\{\sigma_{\ell}\\}}_{\ell=1}^{L}}$ | Mixture kernel, see (6)
$k_{\operatorname{conv}}^{\lambda,k_{i}}$ | Convolution kernel, see (7)
$\mathcal{H}_{k}(\Omega)$ | The native space relating to kernel $k$ on $\Omega$
$V_{1}\cong V_{2}$ |  $V_{1}$ and $V_{2}$ are equal as vector spaces and their respective norms are equivalent
$V_{1}\hookrightarrow V_{2}$ |  $V_{1}$ is continuously embedded into $V_{2}$
$C_{\text{low}}$, $C_{\text{up}}$ | Constants for norm equivalences
$\operatorname{inv}(g)$ | The inverse function of a function $g$
$1/g$ | The reciprocal function of a function $g$
$g{{}^{\prime}}$ | the $1$-st derivative of a function univariate function $g$
$g^{(r)}$ | the $r$-th derivative of a function $g$
$\lceil x\rceil$ | the ceiling of $x\in\mathbb{R}$, the smallest integer greater than or equal to $x$
$\lfloor x\rfloor$ | the floor of $x\in\mathbb{R}$, the greatest integer less than or equal to $x$
$\mathbb{N}_{0}$ | Set of natural numbers $\mathbb{N}$ and $0$
$\mathbb{R}_{+}$ | Set of positive real numbers
$\mathbbm{1}_{A}$ | Indicator function on a set $A$
$\mathbf{I}_{N}$ | Identity matrix of size $N$
$H^{\beta}(\Omega)$ | Sobolev space $W^{\beta,2}(\Omega)$

##  2 GP regression

We want to use GP regression to derive approximations of functions
$f:\Omega\subset\mathbb{R}^{d}\to\mathbb{R}$. We will focus on the convergence
of these approximations as the number of training points $N$ tends to
infinity, but first give an introduction to the general methodology.

###  2.1 Set up

Let $f:\Omega\to\mathbb{R}$ be an arbitrary function, with
$\Omega\subset\mathbb{R}^{d}$ a bounded Lipschitz domain that satisfies an
interior cone condition [62]. Denote by
$U_{N}:=\\{u_{n}\\}_{n=1}^{N}\subset\Omega$ a set of $N\in\mathbb{N}$ distinct
training points where $f(u_{n})$ is observed. Collectively, we denote this
data as $y_{N}:=\\{u_{n},f(u_{n})+\varepsilon_{n}\\}_{n=1}^{N}$, where
$\varepsilon_{n}\sim N(0,\delta^{2})$ i.i.d. for $\delta^{2}\geq 0$. In our
analysis, we will be interested in the setting of noisy data, where
$\delta^{2}>0$, as well as the setting of noise-free data, where
$\varepsilon_{n}\equiv 0$ and formally $\delta^{2}=0$.

To recover $f$ from $y_{N}$ in the Bayesian framework, we first assign a GP
_prior_ to $f$, denoted by

(1) $f_{0}\sim\mathcal{GP}(0,k(u,u^{\prime})).$

Here, $k:\Omega\times\Omega\to\mathbb{R}$ is a positive semi-definite
covariance function (or covariance kernel), that is, for any $n\in\mathbb{N}$,
$c\in\mathbb{R}^{n}\setminus\\{0\\}$ and $u_{1}\ldots,u_{n}\in\Omega$, we have
that111We use the convention of [62] whereby a non-strict inequality is used
to define positive semi-definiteness. The kernel is referred to as positive
definite if strict inequality holds for pairwise distinct
$u_{1}\ldots,u_{n}\in\Omega$.

$\sum_{i=1}^{n}\sum_{j=1}^{n}c_{i}c_{j}k(u_{i},u_{j})\geq 0.$

For ease of presentation, we have chosen the mean function in (1) to be zero,
but all results in this paper extend to the case where a non-zero mean
function is used, under suitable assumptions (see, e.g. [57]).

The prior distribution encapsulates our prior knowledge of the function $f$,
and in particular is independent of the data $y_{N}$. Intuitively, it should
give a higher probability to the types of functions we expect to see, and this
is typically reflected in the choice of mean and covariance function. Typical
choices for the covariance function $k$ are discussed in section 2.2.

We condition the prior on the training data $y_{N}$ to obtain the _posterior_
distribution

(2)
$f_{N}:=f_{0}|y_{N}\sim\mathcal{GP}(m_{N,\delta^{2}}^{f}(u),k_{N,\delta^{2}}(u,u^{\prime})).$

with mean and the covariance function given by (see e.g. [43])

(3) $\displaystyle m_{N,\delta^{2}}^{f}(u)$
$\displaystyle=k(u,U_{N})^{T}(K(U_{N},U_{N})+\delta^{2}\mathrm{I})^{-1}f^{\varepsilon}(U_{N}),$
(4) $\displaystyle k_{N,\delta^{2}}(u,u^{\prime})$
$\displaystyle=k(u,u^{\prime})-k(u,U_{N})^{T}(K(U_{N},U_{N})+\delta^{2}\mathrm{I})^{-1}k(u^{\prime},U_{N}),$

where
$f^{\varepsilon}(U_{N})=\left[f(u_{1})+\varepsilon_{1},\ldots,f(u_{N})+\varepsilon_{N}\right]\in\mathbb{R}^{N}$,
$k(u,U_{N})=\left[k(u,u_{1}),\ldots,k(u,u_{N})\right]\in\mathbb{R}^{N}$ and
$K(U_{N},U_{N})\in\mathbb{R}^{N\times N}$ is the matrix with $ij$th entry
$k(u_{i},u_{j})$. Note that $K(U_{N},U_{N})+\delta^{2}\mathrm{I}$ is
invertible for $\delta^{2}>0$ since we have assumed that $k$ is positive semi-
definite. In the noise-free case, in which the above equations hold with
$\delta^{2}=0$, we require $k$ to be positive definite.

One of the major advantages of using the Bayesian framework outlined above is
that it allows us to perform uncertainty quantification. We are not given a
point estimate for the unknown function $f$ given the data $y_{N}$, but rather
a distribution over a suitable function space. Calculating variances, and
hence error estimates, is therefore possible, and this can be crucial in
applications of GP regression in computational pipelines (see, e.g. [55]).

###  2.2 Covariance kernels

We begin our discussion on covariance kernels by introducing the widely used
family of stationary kernels known as Matérn kernels. Following this, we
discuss various methods to construct non-stationary kernels from stationary
ones.

####  2.2.1 Matérn kernel

The family of Matérn kernels (see e.g. [43]) is given by

$k_{\text{Mat}(\nu)}(u,u^{\prime})=\sigma^{2}\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\sqrt{2\nu}\frac{\|u-u^{\prime}\|_{2}}{\lambda}\right)K_{\nu}\left(\sqrt{2\nu}\frac{\|u-u^{\prime}\|_{2}}{\lambda}\right),$

where $\Gamma$ is the gamma function, $B_{\nu}$ is the modified Bessel
function of the second kind and $\sigma^{2},\lambda,\nu>0$ are positive
parameters. The parameter $\sigma^{2}$ is the marginal variance, $\lambda$ is
the correlation length scale and $\nu$ is the smoothness parameter. Figures 6
and 6 illustrate the effect of changing the correlation length $\lambda$ on
sample paths from a Matérn kernel with $\nu=5/2$. Notice that with a smaller
length scale of $0.5$ the samples fluctuate much more rapidly than with a
larger length scale of $3$.

Figure 1: Sample paths from $\mathcal{GP}(0,k_{\text{Mat}(5/2)})$, with
$\\{\sigma^{2},\lambda\\}=\\{1,0.5\\}$.

Figure 2: Sample paths from $\mathcal{GP}(0,k_{\text{Mat}(5/2)})$, with
$\\{\sigma^{2},\lambda\\}=\\{1,3\\}$.

Figure 3: Sample paths from
$\mathcal{GP}(0,k_{\operatorname{warp}}^{w,k_{s}})$, with stationary kernel
$k_{\text{Mat}(5/2)}$ and $w(u)=(u+1.5)^{4}$.

Figure 4: Sample paths from
$\mathcal{GP}(0,k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{2}})$,
with stationary kernels
$k_{1}=k_{\text{Mat}(\infty)},\,k_{2}=k_{\text{Mat}(0.5)}$ and
$\sigma_{1}(u)=\mathbbm{1}_{\\{u<0\\}}$, $\sigma_{2}(u)=\mathbbm{1}_{\\{u\geq
0\\}}$.

Figure 5: Samples from
$\mathcal{GP}(0,k_{\operatorname{conv}}^{{\lambda_{a}},k_{i}})$ with
stationary Gaussian kernel and
${\lambda_{a}}(u)=\mathbbm{1}_{\\{u<0\\}}+\mathbbm{1}_{\\{u\geq 0\\}}/100$.

Figure 6: Posterior mean with prior kernel a Matérn kernel with $\nu=1/2$,
with ${\lambda}=0.8,0.08,0.05$.

When $\nu=p+1/2$ with $p\in\mathbb{N}_{0}$, the Matérn kernel can be written
as a product of an exponential and a polynomial of order $p$ [1]. For
instance,

  * •

$p=0$ gives
$k_{\text{Mat}(1/2)}(u,u^{\prime})=\sigma^{2}\exp\left(-\frac{\|u-u^{\prime}\|_{2}}{\lambda}\right)$,

  * •

$p=1$ gives
$k_{\text{Mat}(3/2)}(u,u^{\prime})=\sigma^{2}\left(1+\frac{\sqrt{3}\|u-u^{\prime}\|_{2}}{\lambda}\right)\exp\left(-\frac{\sqrt{3}\|u-u^{\prime}\|_{2}}{\lambda}\right)$,

  * •

$p=2$ gives
$k_{\text{Mat}(5/2)}(u,u^{\prime})=\sigma^{2}\left(1+\frac{\sqrt{5}\|u-u^{\prime}\|_{2}}{\lambda}+\frac{5\|u-u^{\prime}\|_{2}^{2}}{3\lambda^{2}}\right)\exp\left(-\frac{\sqrt{5}\|u-u^{\prime}\|_{2}}{\lambda}\right)$.

Additionally, as $\nu\to\infty$ the Matérn covariance kernel converges to the
Gaussian kernel

$k_{\text{Mat}(\infty)}(u,u^{\prime})=\sigma^{2}\exp\left(-\frac{\|u-u^{\prime}\|_{2}^{2}}{2\lambda^{2}}\right).$

Matérn kernels are often the standard choice for GP regression, but can behave
poorly when the observations exhibit non-stationary trends. Figure 6 gives a
simple example of when non-stationary kernels are much better suited. The
function to be approximated is continuous but not continuously differentiable,
and it varies smoothly in the left hand side of the domain, whereas it has
very steep gradients in the right hand side. Figure 6 shows the posterior mean
$m_{N,0}^{f}$, varying the length scale $\lambda$ in the prior with a Matérn
kernel with $\nu=1/2$. We can see that no single length scale works for the
entire function. We desire a large length scale when the function is very
flat, and a small length scale at the three jump points where the gradient is
very large.

####  2.2.2 Non-stationary covariance kernels

As illustrated by the simple example in Figure 6, many applications require a
GP prior that can model different length scales in the observations. We now
introduce three specific constructions for non-stationary kernels that allow
for this. There are of course many more methods known to create non-stationary
kernels, see for instance [44, 37, 40] and the references therein, but we
believe the constructions covered in this paper are prototypical of approaches
that efficiently model separate length scales in the data. It may be possible
to analyse alternative methods using techniques similar to those presented in
this paper.

Before we go on, we introduce here precise definitions of stationarity, non-
stationarity and anisotropy. A kernel $k:\Omega\times\Omega\to\mathbb{R}$ is
stationary if there exists a function $k_{s}:\Omega\to\mathbb{R}$ such that
that $k(u,u)=k_{s}(u-u^{\prime})$. In other words, the covariance between the
GP at $u$ and $u^{\prime}$ does not depend on the location of $u$ and
$u^{\prime}$. A kernel is thus non-stationary if no such $k_{s}$ exists. We
refer to a GP as (non-)stationary if it has a (non-)stationary covariance
kernel. A stationary kernel is isotropic if there exists a function
$k_{i}:[0,\infty)\to\mathbb{R}$ such that that
$k(u,u)=k_{i}(\|u-u^{\prime}\|_{2})$. In an abuse of notation, we shall for
simplicity often write
$k(u,u^{\prime})=k(u-u^{\prime})=k(\|u-u^{\prime}\|_{2})$ when referring to a
stationary and/or isotopic kernel.

#### Warping kernel

The warping kernel [11, 47] uses a homeomorphism to warp the input space of a
GP. For a warping function $w:\Omega\to\Omega$ and a stationary kernel
$k_{s}$, the warping kernel is defined by

(5)
$k_{\operatorname{warp}}^{w,k_{s}}(u,u^{\prime})=k_{s}\left({w(u)-w(u^{\prime})}\right).$

The positive (semi)-definiteness of this kernel follows from a simple
argument.

###### Proposition 2.1

If $k_{s}$ is positive semi-definite, then $k_{\operatorname{warp}}^{w,k_{s}}$
is positive semi-definite. If $w$ is injective and $k_{s}$ is positive
definite, then $k_{\operatorname{warp}}^{w,k_{s}}$ is positive definite.

##### Proof:

We have

$\sum_{i=1}^{n}\sum_{j=1}^{n}c_{i}c_{j}k_{\operatorname{warp}}^{w,k_{s}}(u^{i},u^{j})=\sum_{i=1}^{n}\sum_{j=1}^{n}c_{i}c_{j}k_{s}(w(u^{i})-w(u^{j}))=\sum_{j=1}^{n}c_{i}c_{j}k_{s}(\tilde{u}^{i}-\tilde{u}^{j}),$

where $\tilde{u}^{i}=w(u^{i})$, for $i=1,\dots,n$. The claim then follows.
$\square$

See Figure 6 for an example of sample paths of a GP with a warping kernel.
Notice that the samples fluctuate much more in the positive than the negative
part of the domain. This happens because the negative part is “pushed
together” by the warping function $w(u)=(u+1.5)^{4}$, while the positive part
of the domain is “pulled apart”. For example,
$|w(0.75)-w(0.5)|>|0.25|>|w(-0.5)-w(-0.75)|$, and so for this non-stationary
kernel the correlation between $(u,u^{\prime})=(0.75,0.5)$ will be lower than
for $(u,u^{\prime})=(-0.5,-0.75)$.

#### Kernel mixture

The kernel mixture formulation [17, 60] in general allows the regularity,
variance and length scales to change over the domain. For $L\in\mathbb{N}$,
the non-stationary covariance kernel mixture is defined as

(6)
$k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}}(u,u^{\prime}):=\sum_{\ell=1}^{L}\sigma_{\ell}(u)\sigma_{\ell}(u^{\prime})k_{\ell}\left(u,u^{\prime}\right),$

where the functions $\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}$ represent
mixture component coefficients and stationary kernels, respectively. The
kernel mixture is positive (semi-)definite under mild assumptions.

###### Proposition 2.2

Suppose $\sigma_{\ell}:\Omega\to\mathbb{R}$ satisfies $\sigma_{\ell}(u)\geq 0$
for all $u\in\Omega$ and $\ell=1,\dots,L$. If $k_{\ell}$ is positive semi-
definite for all $\ell=1,\dots,L$, then
$k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}}$ is
positive semi-definite. If there exists $j\in\\{1,\ldots,L\\}$ such that
$\sigma_{j}(u)>0$ for all $u\in\Omega$ and $k_{j}$ is positive definite, then
$k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}}$ is
positive definite.

##### Proof:

We have

$\sum_{i,j=1}^{n}c_{i}c_{j}k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}}(u^{i},u^{j})=\sum_{i,j=1}^{n}c_{i}c_{j}\sum_{\ell=1}^{L}\sigma_{\ell}(u^{i})\sigma_{\ell}(u^{j})k_{\ell}(u^{i},u^{j})=\sum_{\ell=1}^{L}\sum_{i,j=1}^{n}\tilde{c}_{i,\ell}\tilde{c}_{j,\ell}k_{\ell}(u^{i},u^{j})$

where $\tilde{c}_{i,\ell}=c_{i}\sigma_{\ell}(u^{i})$, for $i=1,\dots,n$ and
$\ell=1,\dots,L$. The claim then follows. $\square$

To model separate length scales, the kernel mixture could be constructed such
that each $k_{\ell}$ is a Matérn kernel with a different length scale
$\lambda_{\ell}$. See Figure 6 for an example of a GP with a mixture kernel.
In this example $L=2$, with two stationary kernels
$k_{1}=k_{\text{Mat}(\infty)}$ and $k_{2}=k_{\text{Mat}(0.5)}$, and non-
stationary functions $\sigma_{1},\sigma_{2}:\Omega\to\mathbb{R}$ defined
$\sigma_{1}(u)=\mathbbm{1}_{\\{u<0\\}}$, $\sigma_{2}(u)=\mathbbm{1}_{\\{u\geq
0\\}}$. Notice how the regularity of the samples changes drastically at $u=0$.
This is exactly the behaviour we would expect, since when $u,u^{\prime}>0$ the
only kernel not taking non zero values is $k_{\text{Mat}(0.5)}$ and similarly
for $u,u^{\prime}<0$ and the Gaussian kernel $k_{\text{Mat}(\infty)}$.

In the case $L=1$, the mixture kernel is no longer a mixture, but can instead
be interpreted as a version of the stationary kernel $k_{1}$ with non-
stationary marginal variance $\sigma_{1}(u)^{2}$.

#### Convolution kernel

Paciorek [41] gives a general method to construct anisotropic versions of
isotropic covariance kernels, based on the convolution of two Gaussian
kernels. The resulting convolution kernel allows non-stationarity in the form
of different length scales in the domain. For $k_{i}$ an isotropic covariance
function and ${\lambda_{a}}:\Omega\to\mathbb{R}_{+}$ representing the length
scale, we define the non-stationary convolution covariance kernel by

(7)
$k_{\operatorname{conv}}^{{\lambda_{a}},k_{i}}(u,u^{\prime})=\frac{2^{d/2}{\lambda_{a}}(u)^{d/4}{\lambda_{a}}(u^{\prime})^{d/4}}{({\lambda_{a}}(u)+{\lambda_{a}}(u^{\prime}))^{d/2}}k_{i}\left(\frac{\|u-u^{\prime}\|_{2}}{\sqrt{({\lambda_{a}}(u)+{\lambda_{a}}(u^{\prime}))/2)}}\right).$

Crucially, ${\lambda_{a}}$ is a function rather than a scalar parameter as in
the Matérn covariance kernel.

###### Proposition 2.3

Suppose ${\lambda_{a}}:\Omega\to\mathbb{R}$ satisfies ${\lambda_{a}}(u)>0$ for
all $u\in\Omega$. If $k_{i}$ is positive (semi-) definite, then
$k_{\operatorname{conv}}^{{\lambda_{a}},k_{i}}$ is positive (semi-)definite.

##### Proof:

This follows from [41, Theorem 1] and [13, Propostion 1], since the matrix
${\lambda_{a}}(u)\mathrm{I}$ is positive definite for all $u\in\Omega$.
$\square$

Figure 6 gives an example for the convolution kernel being used as the kernel
for a GP prior. The isotropic kernel in this example is the Gaussian kernel,
and the non-stationary function ${\lambda_{a}}:[-1,1]\to\mathbb{R}$ is given
by ${\lambda_{a}}(u)=\mathbbm{1}_{\\{u<0\\}}+\mathbbm{1}_{\\{u\geq 0\\}}/100$.
For negative $u$ we see that ${\lambda_{a}}(u)=1$, meaning that there is a
high correlation between any two points $u,u^{\prime}<0$, whereas for positive
$u$ we see that ${\lambda_{a}}(u)=1/100$ and so there is a low correlation
between any two points $u,u^{\prime}>0$.

###  2.3 Estimating hyper-parameters

To obtain good practical performance, it is crucial to use suitable choices of
the functions $w$, $\\{\sigma_{\ell}\\}_{\ell=1}^{L}$ and $\lambda_{a}$
defining the non-stationary structure of the warping, mixture and convolution
kernels, respectively. Good choices are typically not known a-priori, and
these functions are usually learnt from the observations $y_{N}$. Including
the hyper-parameters appearing in the stationary and/or isotropic kernels used
in the constructions (5)-(7), we then need to estimate a set of hyper-
parameters $\theta$, containing both scalars and functions.

In the following, we discuss two options for dealing with the hyper-
parameters. Firstly, we present an empirical Bayesian approach, whereby a
point-estimate of the hyper-parameters is computed, and then “plugged-in” to
the GP prior before we perform the standard prior to posterior update.
Secondly, we discuss a particular hierarchical Bayesian approach, where a GP
prior is placed on the functional hyper-parameters resulting in a DGP model.
We are then interested in the marginal posterior on the unknown function $f$
given $y_{N}$.

For a fixed value of the hyper-parameters $\theta$, recall the prior (1) on
$f$ given by

$f_{0}|\theta\sim\mathcal{GP}(0,k(\cdot,\cdot;\theta)),$

where we have now explicitly indicated the dependence on $\theta$. In the
empirical Bayesian approach, we use the available data $y_{N}$ to estimate
suitable values for the hyper-parameters, denoted by $\widehat{\theta}_{N}$,
and then use the conditional posterior

(8)
$f_{N}(\widehat{\theta}_{N}):=f_{0}|\widehat{\theta}_{N},y_{N}\sim\mathcal{GP}(m_{N,0}^{f}(\cdot;\widehat{\theta}_{N}),k_{N}(\cdot,\cdot;\widehat{\theta}_{N}))$

to recover $f$. Here, $m_{N,0}^{f}$ and $k_{N}$ are as in (3)-(4), with the
dependence on $\theta$ made explicit. Note that this distribution is given in
closed form, and in particular, its mean and covariance functions are known
analytically. This approach is discussed in more detail in section 5.4.

In a hierarchical approach, we instead place a prior distribution on $\theta$.
A particular example is a DGP, which is constructed using a sequence of random
fields $\\{f^{n}\\}_{n\in\mathbb{N}_{0}}$ that are conditionally Gaussian [11,
13]:

$\displaystyle f^{0}$ $\displaystyle\sim\mathcal{GP}(0,k_{0}(u,u^{\prime})),$
$\displaystyle f^{n+1}|f^{n}$
$\displaystyle\sim\mathcal{GP}(0,k_{n}(u,u^{\prime};f^{n})).$

We denote by $f^{D}$ a DGP of depth $D$. In the context of constructing priors
for regression, we can then use the DGP $f^{D}$ as a prior on $f$. The
penultimate layer $f^{D-1}$ will be used to construct the non-stationarity in
$f^{D}$, for example by defining the warping function $w=F(f^{D-1})$ or the
length scale function ${\lambda_{a}}=F(f^{D-1})$. The remaining layers define
the prior on $f^{D-1}$. In contrast to standard GP regression, the posterior
$f^{D}|y_{N}$ can no longer be found explicitly, and one typically uses Markov
chain Monte Carlo (MCMC) methods for sampling [13, 29]. DGPs are discussed in
more detail in section 6.

##  3 Reproducing kernel Hilbert spaces

A crucial ingredient in the analysis of GP and DGP regression are reproducing
kernel Hilbert spaces (RKHSs) [4, 42], also referred to as native spaces [62]
or Cameron-Martin spaces [10] (in the context of Gaussian measures). This
section is devoted to studying the RKHSs corresponding to the non-stationary
kernels $k^{w,k_{s}}_{\mathrm{warp}}$ and
$k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}}$
introduced in section 2.2.2. In particular, we show that these are norm-
equivalent to Sobolev spaces in special cases. We were unable to find an
explicit expression for the RKHS of the convolution kernel
$k^{{\lambda_{a}},k_{i}}_{\mathrm{conv}}$.

Let us start with the definition of an RKHS.

###### Definition 3.1

Let $\mathcal{H}(\Omega)$ be a real Hilbert space of functions
$f:\Omega\to\mathbb{R}$. $\mathcal{H}(\Omega)$ is the RKHS of a kernel
$k:\Omega\times\Omega\to\mathbb{R}$ if

  * (i)

$k(\cdot,u)\in\mathcal{H}(\Omega)$ for all $u\in\Omega$,

  * (ii)

$f(u)=\left(f,k(\cdot,u)\right)_{\mathcal{H}}(\Omega)$ for all
$f\in\mathcal{H}(\Omega)$ and all $u\in\Omega$.

In the remainder of this work, we shall use $\mathcal{H}_{k}(\Omega)$ to
denote the RKHS corresponding to a kernel $k$. It is not always easy to find
the RKHS of a specific kernel in closed form. A notable exception is the
family of Matérn covariance kernels defined in section 2.2.1, for which the
RKHS can be characterised by its norm-equivalence with a Sobolev space [62,
57].

######  Proposition 3.2 ([57, Lemma 3.4])

We have $\mathcal{H}_{k_{\text{Mat}(\nu)}}(\Omega)\cong H^{\nu+d/2}(\Omega)$.
In particular,

$C_{\text{low}}(k_{\text{Mat}(\nu)})\|g\|_{\mathcal{H}_{k_{\text{Mat}(\nu)}}(\Omega)}\leq\|g\|_{H^{\nu+d/2}(\Omega)}\leq
C_{\text{up}}(k_{\text{Mat}(\nu)})\|g\|_{\mathcal{H}_{k_{\text{Mat}(\nu)}}(\Omega)},$

for all $g\in H^{\nu+d/2}(\Omega)$, where

$\displaystyle C_{\text{low}}(k_{\text{Mat}(\nu)})$
$\displaystyle=\frac{\pi^{d/2}\Gamma(\nu)}{\sigma^{2}\Gamma(\nu+d/2)\lambda^{d}}\max\\{1,\lambda^{2}\\},$
$\displaystyle C_{\text{up}}(k_{\text{Mat}(\nu)})$
$\displaystyle=\frac{\sigma^{2}\Gamma(\nu+d/2)\lambda^{d}}{\pi^{d/2}\Gamma(\nu)}\max\\{1,\lambda^{-2}\\}.$

Note that the only hyper-parameter that affects the order of the Sobolev space
is $\nu$, whereas $\lambda$ and $\sigma^{2}$ also affect the constants in the
norm equivalence.

###  3.1 Operations on kernels

For the warping kernel $k^{w,k_{s}}_{\mathrm{warp}}$ and the kernel mixture
$k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}}$, we can
use algebraic operations on kernels to construct the corresponding RKHSs
explicitly. We recall the following results from [4, 42].

###### Proposition 3.3

[42, Theorem 5.7] Let $k:\Omega\times\Omega\to\mathbb{R}$ and
$w:\Omega\to\Omega$ be given. Then the RKHS corresponding to the kernel
$k^{w}(u,u^{\prime})=k\left({w(u),w(u^{\prime})}\right)$ is given by

$\mathcal{H}_{k^{w}}(\Omega)=\left\\{h\circ
w:h\in\mathcal{H}_{k}(\Omega)\right\\}$

equipped with the norm
$\|u\|_{\mathcal{H}_{k^{w}}}=\min\\{\|h\|_{\mathcal{H}_{k}}:u=h\circ w\\}$.

###### Proposition 3.4

[42, Proposition 5.20] Let $k:\Omega\times\Omega\to\mathbb{R}$ and
$\sigma:\Omega\to\mathbb{R}$ be given, and define

$\mathcal{H}_{k}^{0}(\Omega):=\\{h\in\mathcal{H}_{k}(\Omega):\sigma h\equiv
0\\}.$

Then the RKHS corresponding to the kernel
$k^{\sigma}(u,u^{\prime})=\sigma(u)\sigma(u^{\prime})k(u,u^{\prime})$ is given
by

$\mathcal{H}_{k^{\sigma}}(\Omega)=\left\\{\sigma
h:h\in{\left(\mathcal{H}_{k}^{0}(\Omega)\right)}^{\perp}\right\\},$

equipped with the norm $\|\sigma
h\|_{\mathcal{H}_{k^{\sigma}}}=\|h\|_{\mathcal{H}_{k}}$.

###### Remark 3.5

Note that if we in Proposition 3.4 additionally assume $\sigma(u)>0$ for all
$u\in\Omega$, then $\sigma h\equiv 0$ if and only if $h\equiv 0$ and hence
$\mathcal{H}_{k^{\sigma}}(\Omega)=\\{\sigma h:h\in\mathcal{H}_{k}(\Omega)\\}$.

###### Proposition 3.6

[42, Theorem 5.4] Let $k_{1}:\Omega\times\Omega\to\mathbb{R}$ and
$k_{2}:\Omega\times\Omega\to\mathbb{R}$ be given. Then the RKHS corresponding
to the kernel $k=k_{1}+k_{2}$ is given by

$\mathcal{H}_{k}(\Omega):=\\{u_{1}+u_{2}:u_{i}\in\mathcal{H}_{k_{i}}(\Omega)\,i=1,2\\}$

equipped with the norm

$\|u\|^{2}_{\mathcal{H}_{k}(\Omega)}=\min\\{\|u_{1}\|^{2}_{\mathcal{H}_{k_{1}}(\Omega)}+\|u_{2}\|^{2}_{\mathcal{H}_{k_{2}}(\Omega)}:u=u_{1}+u_{2},u_{i}\in\mathcal{H}_{k_{i}}(\Omega),i=1,2\\}.$

###  3.2 RKHSs of non-stationary kernels as Sobolev spaces

The results presented in the previous section allow us to derive the RKHSs of
the non-stationary kernels $k^{w,k_{s}}_{\mathrm{warp}}$ and
$k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}}$
explicitly in terms of the RKHS of the stationary/isotropic kernel used to
construct them and the hyper-parameters $w$ or
$\\{\sigma_{\ell}\\}_{\ell=1}^{L}$, respectively. By Proposition 3.2, we know
that the native space associated with the Matérn kernel can be characterised
by its equivalence with a Sobolev space, and we use this fact to prove the
following equivalences.

######  Theorem 3.7 (RKHS for warping kernel)

Let $\Omega\subset\mathbb{R}$. Suppose $\\{w,k_{s}\\}$ are such that $(i)$
$k_{s}=k_{\text{Mat}(\nu)}$ with $\nu=p+1/2$ for some $p\in\mathbb{N}$, $(ii)$
$w\in C^{\beta}(\Omega)$ for $\beta=\nu+1/2$, and $(iii)$
$|w{{}^{\prime}}|\geq c>0$ and $|\operatorname{inv}(w){{}^{\prime}}|\geq
c^{\prime}>0$ a.e. in $\Omega$. Then
$\mathcal{H}_{k_{\operatorname{warp}}^{w,k_{s}}}(\Omega)\cong
H^{\beta}(\Omega)$. In particular,

$C_{\text{low}}(k_{\operatorname{warp}}^{w,k_{s}})\|g\|_{\mathcal{H}_{k_{\operatorname{warp}}^{w,k_{s}}(\Omega)}}\leq\|g\|_{H^{\beta}(\Omega)}\leq
C_{\text{up}}(k_{\operatorname{warp}}^{w,k_{s}})\|g\|_{\mathcal{H}_{k_{\operatorname{warp}}^{w,k_{s}}(\Omega)}}$

for all $g\in H^{\beta}(\Omega)$, where for a constant $C(\beta)$ depending
only $\beta$ we have

$\displaystyle C_{\text{low}}({k_{\operatorname{warp}}^{w,k_{s}}})$
$\displaystyle=C_{\text{low}}(k_{s})C(\beta)^{-1}\sqrt{c^{\prime}}\max\left\\{1,\|\operatorname{inv}(w)\|_{C^{\beta}(\Omega)}^{\beta}\right\\}^{-1}$
$\displaystyle C_{\text{up}}({k_{\operatorname{warp}}^{w,k_{s}}})$
$\displaystyle=C_{\text{up}}(k_{s})C(\beta)\frac{1}{\sqrt{c}}\max\left\\{1,\|w\|_{C^{\beta}(\Omega)}^{\beta}\right\\}.$

######  Theorem 3.8 (RKHS for mixture kernel)

Let $j\in\\{1,\dots,L\\}$ be such that $\nu_{j}=\min_{1\leq\ell\leq
L}\nu_{\ell}$. Suppose $\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}$ are such
that for $\ell=1,\dots,L$ we have $(i)$ $k_{\ell}=k_{\text{Mat}(\nu_{\ell})}$,
$(ii)$ $\sigma_{\ell}\in H^{\beta}(\Omega)$ for $\beta=\nu_{j}+d/2$ and
$(iii)$ $1/\sigma_{j}\in H^{\beta}(\Omega)$. Then
$\mathcal{H}_{k_{\operatorname{mix}}^{{\\{\sigma_{\ell},k_{\ell}\\}}_{\ell=1}^{L}}}(\Omega)\cong
H^{\beta}(\Omega)$. In particular,

$C_{\text{low}}(k_{\operatorname{mix}}^{{\\{\sigma_{\ell},k_{\ell}\\}}_{\ell=1}^{L}})\|g\|_{\mathcal{H}_{k_{\operatorname{mix}}^{{\\{\sigma_{\ell},k_{\ell}\\}}_{\ell=1}^{L}}}(\Omega)}\leq\|g\|_{H^{\beta}(\Omega)}\leq
C_{\text{up}}(k_{\operatorname{mix}}^{{\\{\sigma_{\ell},k_{\ell}\\}}_{\ell=1}^{L}})\|g\|_{\mathcal{H}_{k_{\operatorname{mix}}^{{\\{\sigma_{\ell},k_{\ell}\\}}_{\ell=1}^{L}}}(\Omega)}$

for all $g\in H^{\beta}(\Omega)$, where

$\displaystyle
C_{\text{low}}(k_{\operatorname{mix}}^{{\\{\sigma_{\ell},k_{\ell}\\}}_{\ell=1}^{L}})$
$\displaystyle=C_{\ref{prop:banachalg}}^{-1}C_{\text{low}}(k_{j})\|1/\sigma_{j}\|_{H^{\beta}}^{-1},$
$\displaystyle
C_{\text{up}}(k_{\operatorname{mix}}^{{\\{\sigma_{\ell},k_{\ell}\\}}_{\ell=1}^{L}})$
$\displaystyle=\sqrt{2}C_{\ref{prop:banachalg}}\max_{{1\leq\ell\leq
L}}\\{C_{\text{up}}(k_{\ell})\|\sigma_{\ell}\|_{H^{\beta}(\Omega)}\\}.$

###### Remark 3.9

In Theorem 3.7, it is possible to replace the assumption that $w\in
C^{\beta}(\Omega)$ with $w\in W^{\beta,\infty}(\Omega)$. If additionally
$\operatorname{inv}(w)\in W^{\beta,\infty}(\Omega)$, the conclusions of the
Theorem still hold, with the obvious replacement of norms.

###### Remark 3.10

In Theorem 3.8, it might be possible to change the assumption that
$\sigma_{j}(u)\geq\sigma_{\min}$ (and hence $1/\sigma_{j}\in
H^{\beta}(\Omega)$), where $\nu_{j}=\min_{1\leq\ell\leq L}\nu_{\ell}$, to
assuming that for each $u\in\Omega$ there exists $\sigma_{\ell}$ such that
$\sigma_{\ell}(u)\geq\sigma_{\min}$. This has been not been investigated for
brevity. Note however, that we would still require that each $\sigma_{\ell}\in
H^{\beta}(\Omega)$, and thus, that each $\sigma_{\ell}$ is continuous. This
would preclude examples where $\sigma_{\ell}$ are piecewise constant.

Theorems 3.7 and 3.8 provide assumptions on the non-stationary hyper-
parameters $w$ and $\\{\sigma_{\ell}\\}_{\ell=1}^{L}$ such that the RKHS of
the non-stationary kernel is identical (up to norm-equivalence) to the RKHS of
the Matérn kernel used to construct it. Choosing $w$ and
$\\{\sigma_{\ell}\\}_{\ell=1}^{L}$ more regular will not change the RKHS,
since we are limited by the regularity of $k_{s}$ and
$\\{k_{\ell}\\}_{\ell=1}^{L}$. Choosing $w$ and
$\\{\sigma_{\ell}\\}_{\ell=1}^{L}$ less regular on the other hand is
unfortunately not within the framework of this theory, and the RKHS is no
longer a Sobolev space. It is a subspace of a Sobolev space, dependent on $w$
or $\\{\sigma_{\ell}\\}_{\ell=1}^{L}$, respectively. To see this suppose, for
instance, that $L=1$, $k_{1}=k_{\text{Mat}(\nu_{1})}$ for $\beta=\nu_{1}+d/2$
such that $\mathcal{H}_{k_{1}}(\Omega)\cong H^{\beta}(\Omega)$, and
$\sigma_{1}\in H^{\alpha}(\Omega)$ with $\alpha<\beta$. By Proposition 3.4, we
know that every element of
$\mathcal{H}_{k_{\operatorname{mix}}^{\sigma_{1},k_{1}}}$ can be written in
the form $\sigma_{1}h$ for $h\in H^{\beta}(\Omega)$, which is not true for all
elements of $H^{\alpha}(\Omega)$.

It is possible to extend Theorem 3.7 to the case of general $\nu>0$ and
$\Omega\subset\mathbb{R}^{d}$ following the same proof technique, but this
becomes very technical and has been omitted for brevity. In the statement of
Theorem 3.7, we have further implicitly used the following result, which
allows us to conclude that $\operatorname{inv}(w)\in C^{\beta}(\Omega)$ under
the assumptions of Theorem 3.7.

###### Lemma 3.11

Suppose $\Omega\subseteq\mathbb{R}$, $w\in C^{\beta}(\Omega)$ for
$\beta\in\mathbb{N}$, and $w^{\prime}(u)\geq w^{\prime}_{\min}$ for all
$u\in\Omega$, for some $1\geq w^{\prime}_{\min}>0$. Then
$\operatorname{inv}(w)\in C^{\beta}(\Omega)$ and

$\|\operatorname{inv}(w)\|_{C^{\beta}(\Omega)}\leq\left(1+\frac{1}{(w^{\prime}_{\min})^{2\beta}}\right)(\beta+1)(2\beta-1)!B_{2\beta-1}\left(1+\left\|w\right\|^{\beta-1}_{C^{\beta}(\Omega)}\right),$

where $B_{2\beta-1}$ is the Bell number, the count of the possible partitions
of the set $\\{1,\ldots,2\beta-1\\}$.

##### Proof:

By the inverse function rule (see, e.g. [24]), we see that

(9)
$\operatorname{inv}(w)^{\prime}(u)=\frac{1}{w^{\prime}(\operatorname{inv}(w)(u))}.$

Then $w^{\prime}\in C^{0}(\Omega)$, $\operatorname{inv}(w)\in C^{0}(\Omega)$
and $w^{\prime}(u)\geq w^{\prime}_{\min}$ together imply that
$\operatorname{inv}(w)\in C^{1}(\Omega)$. Using this argument inductively we
see that $\operatorname{inv}(w)\in C^{\beta}(\Omega)$.

For the second claim, note that an explicit expression the $\alpha$-th
derivative of the inverse function is given by

$\operatorname{inv}(w)^{(\alpha)}(u)=\sum_{i=0}^{\alpha-1}\frac{1}{(w^{\prime}(\operatorname{inv}(w)(u)))^{\alpha+i}}(-1)^{i}\frac{1}{i!}\sum_{\begin{subarray}{c}b_{1}+\dots+b_{i}=\alpha+i-1\\\
b_{j}\geq
2\end{subarray}}\binom{\alpha+i-1}{b_{1},\ldots,b_{\alpha-1}}\prod_{1\leq
i\leq\alpha-1}w^{(b_{i})}(u),$

for $u\in\Omega$, where $b_{1},\ldots,b_{i}$ are the sizes of the blocks of
all partitions of size $i$ of the set $\\{1,\ldots,\alpha+i-1\\}$ [24]. First
notice that maximum size of any individual block over all
$i=1,\ldots,\alpha-1$ and over all $i$-tuples $(b_{1},b_{2},\ldots,b_{i})$ is
$\alpha$. This occurs when $i=1$ so that $b_{1}=\alpha-i+1=\alpha$. Hence we
can see that for any $b_{i}$ in any such $i$-tuple we have that

(10)
$\displaystyle\|w^{(b_{i})}\|_{C^{0}(\Omega)}\leq\sup_{\tilde{\alpha}\leq\alpha}\|w^{(\tilde{\alpha})}\|_{C^{0}(\Omega)}=\left\|w\right\|_{C^{\alpha}(\Omega)}.$

Using the fact that $w^{\prime}$ is bounded away from zero, the triangle
inequality and (10) we see that

$\displaystyle\|\operatorname{inv}(w)^{(\alpha)}\|_{C^{0}(\Omega)}$
$\displaystyle\leq\left(1+\frac{1}{(w^{\prime}_{\min})^{2\alpha}}\right)(2\alpha-1)!\sum_{\begin{subarray}{c}b_{1}+\dots+b_{i}=\alpha+i-1\\\
b_{j}\geq 2\end{subarray}}\prod_{1\leq
i\leq\alpha-1}\left\|w^{(b_{i})}\right\|_{C^{0}(\Omega)}$
$\displaystyle\leq\left(1+\frac{1}{(w^{\prime}_{\min})^{2\alpha}}\right)(2\alpha-1)!\left\|w\right\|^{\alpha-1}_{C^{\alpha}(\Omega)}\sum_{\begin{subarray}{c}b_{1}+\dots+b_{i}=\alpha+i-1\\\
b_{j}\geq 2\end{subarray}}1$
$\displaystyle\leq\left(1+\frac{1}{(w^{\prime}_{\min})^{2\alpha}}\right)(2\alpha-1)!B_{2\alpha-1}(1+\left\|w\right\|^{\alpha-1}_{C^{\alpha}(\Omega)}),$

where $B_{2\alpha-1}$ is the Bell number which gives the count of the possible
partitions of the set $\\{1,\ldots,2\alpha-1\\}$. Since $w^{\prime}_{\min}\leq
1$ the following then achieves the desired result,

$\displaystyle\|\operatorname{inv}(w)\|_{C^{\beta}(\Omega)}$
$\displaystyle=\sum_{\alpha=0}^{\beta}\|\operatorname{inv}(w)^{(\alpha)}\|_{C^{0}(\Omega)}\leq\left(1+\frac{1}{(w^{\prime}_{\min})^{2\beta}}\right)(\beta+1)(2\beta-1)!B_{2\beta-1}\left(1+\left\|w\right\|^{\beta-1}_{C^{\beta}(\Omega)}\right).$

$\square$

In Theorem 3.8, to ensure that $1/\sigma_{j}\in H^{\beta}(\Omega)$, where
$\nu_{j}=\min_{1\leq\ell\leq L}\nu_{\ell}$, it is sufficient that $\sigma_{j}$
is uniformly bounded away from zero, $\sigma_{j}\geq\sigma_{min}>0$. The
following Lemma proves this claim along with two other results that will be
used in section 6. In the context of the mixture kernel
$k_{\operatorname{mix}}^{{\\{\sigma_{\ell},k_{\ell}\\}}_{\ell=1}^{L}}$, where
$\sigma_{\ell}$ represents a mixture coefficient, the assumption that
$\sigma_{\ell}$ is non-negative is natural.

###### Lemma 3.12

Suppose $\Omega\subset\mathbb{R}^{d}$. If $\sigma\in H^{\beta}(\Omega)$ for an
integer $\beta>d/2$ and $F\in C^{\infty}(\mathbb{R})$ such that $F\geq
F_{min}>0$ there exists constants $C_{\ref{lemma:reciprocal has
derivatives}}^{\prime}<\infty$ and $C_{\ref{lemma:reciprocal has
derivatives}}^{\prime\prime}=C_{\ref{lemma:reciprocal has
derivatives}}^{\prime\prime}(F_{\min})$ such that

$\displaystyle\|F(\sigma)\|_{H^{\beta}(\Omega)}$ $\displaystyle\leq
C_{\ref{lemma:reciprocal has
derivatives}}^{\prime}\left(1+\|\sigma\|^{\beta}_{H^{\beta}(\Omega)}\right),$
$\displaystyle\|1/F(\sigma)\|_{H^{\beta}(\Omega)}$ $\displaystyle\leq
C_{\ref{lemma:reciprocal has
derivatives}}^{\prime\prime}\left(1+\|\sigma\|^{\beta}_{H^{\beta}(\Omega)}\right).$

In particular, suppose $\sigma\geq\sigma_{\min}>0$, then there exists a
constant $C_{\ref{lemma:reciprocal has derivatives}}=C_{\ref{lemma:reciprocal
has derivatives}}(\sigma_{\min})<\infty$ such that

$\displaystyle\|1/\sigma\|_{H^{\beta}(\Omega)}\leq C_{\ref{lemma:reciprocal
has derivatives}}\left(1+\|\sigma\|^{\beta}_{H^{\beta}(\Omega)}\right).$

##### Proof:

Since $F\in C^{\infty}(\Omega)$ and that for all $k\geq 1$ we have that
$\|F^{(k)}\|_{C^{0}{([F_{\mathrm{min}},\infty))}}<\infty$ the first result
follows from [35, Lemma 29] by restricting the domain to $\Omega$. To show the
second claim notice that $1/F(\sigma)\in C^{\infty}([F_{\min},\infty])$ as a
function of $F(\sigma)$, and further
$\|(1/F(\sigma))^{(k)}\|_{C^{0}{([F_{\min},\infty))}}<\infty$ for all $k\geq
1$. The second result follows again from [35, Lemma 29] by restricting the
domain to $[F_{\min},\infty)$. The third claim follows directly from the first
two cases. $\square$

The remainder of this section is devoted to the proofs of Theorems 3.7 and
3.8. The following proposition will be used in the proof of Theorem 3.8.

###### Proposition 3.13

[2, Theorem 4.39] The Sobolev space $H^{\beta}(\Omega)$,
$\Omega\subset\mathbb{R}^{d}$, is a Banach algebra whenever $2\beta>d$. In
particular, there exists a constant $C_{\ref{prop:banachalg}}$ such that for
all $g,h\in H^{\beta}(\Omega)$ we have

$\|gh\|_{H^{\beta}(\Omega)}\leq
C_{\ref{prop:banachalg}}\|g\|_{H^{\beta}(\Omega)}\|h\|_{H^{\beta}(\Omega)}.$

##### Proof:

[Of Theorem 3.7] We provide the proof for $\beta=2$, which corresponds to
$k_{s}=k_{\text{Mat}(3/2)}$. The general case can be proved similarly using
Faà di Bruno’s formula [25].

We first show some initial results that will help us later. We use the change
of variables $w(x)=y$, so that
$dx=\frac{dx}{dy}dy=w^{\prime}(\operatorname{inv}(w)(y))dy$. Since
$w(\Omega)=\Omega$ (which follows from $\operatorname{inv}(w)$ being well-
defined on $\Omega$) and $|w^{\prime}|\geq c>0$ a.e., we have

$\displaystyle\|u\circ w\|_{L^{2}(\Omega)}^{2}$
$\displaystyle=\int_{w(\Omega)}(u(y))^{2}\frac{dy}{w^{\prime}(\operatorname{inv}(w)(y))}\leq\frac{1}{c}\|u\|_{L^{2}(\Omega)}^{2},$
$\displaystyle\|(u^{\prime}\circ w)w^{\prime}\|_{L^{2}(\Omega)}^{2}$
$\displaystyle=\int_{w(\Omega)}(u^{\prime}(y)w^{\prime}(\operatorname{inv}(w)(y)))^{2}\frac{dy}{w^{\prime}(\operatorname{inv}(w)(y))}\leq\frac{1}{c}\|w\|^{2}_{W^{1,\infty}(\Omega)}\|u^{\prime}\|_{L^{2}(\Omega)}^{2},$
$\displaystyle\|(u^{\prime\prime}\circ
w)(w^{\prime})^{2}\|_{L^{2}(\Omega)}^{2}$
$\displaystyle\leq\frac{1}{c}\|w\|^{4}_{W^{1,\infty}(\Omega)}\|u^{\prime\prime}\|_{L^{2}(\Omega)}^{2},$
$\displaystyle\|(u^{\prime}\circ w)w^{\prime\prime}\|_{L^{2}(\Omega)}^{2}$
$\displaystyle=\int_{w(\Omega)}((u^{\prime}(y))w^{\prime\prime}(x))^{2}\frac{dy}{w^{\prime}(\operatorname{inv}(w)(y))}\leq\frac{1}{c}{\|w\|^{2}_{W^{2,\infty}(\Omega)}}\|u^{\prime}\|_{L^{2}(\Omega)}^{2}.$

First, we show that
$\mathcal{H}_{k_{\operatorname{warp}}^{w,k_{s}}}(\Omega)\hookrightarrow
H^{2}(\Omega)$, where $\hookrightarrow$ denotes a continuous embedding. By
Proposition 3.3, for arbitrary
$g\in\mathcal{H}_{k_{\operatorname{warp}}^{w,k_{s}}}(\Omega)$ there exists
$u\in\mathcal{H}_{k_{s}}(\Omega)\cong H^{2}{(\Omega)}$ such that $g=u\circ w$
and
$\|g\|_{\mathcal{H}_{k_{\operatorname{warp}}^{w,k_{s}}}(\Omega)}=\|u\|_{\mathcal{H}_{k_{s}}(\Omega)}$.
Using the chain rule, the bounds proven above, and Proposition 3.2, we see
that

$\displaystyle\|g\|^{2}_{H^{2}(\Omega)}$ $\displaystyle=\|u\circ
w\|_{H^{2}(\Omega)}^{2}$
$\displaystyle\leq\frac{1}{c}\max\left\\{1,\|w\|_{W^{2,\infty}(\Omega)}^{2}\right\\}^{2}\|u\|_{H^{2}(\Omega)}^{2}$
$\displaystyle\leq
C_{\text{up}}(k_{s})^{2}\frac{1}{c}\max\left\\{1,\|w\|_{W^{2,\infty}(\Omega)}^{2}\right\\}^{2}\|u\|_{\mathcal{H}_{k_{s}}(\Omega)}^{2}$
$\displaystyle=C_{\text{up}}(k_{s})^{2}\frac{1}{c}\max\left\\{1,\|w\|_{W^{2,\infty}(\Omega)}^{2}\right\\}^{2}\|g\|_{\mathcal{H}_{k_{\operatorname{warp}}^{w,k_{s}}}(\Omega)}^{2}.$

Conversely, to show that that
$H^{2}(\Omega)\hookrightarrow\mathcal{H}_{k_{\operatorname{warp}}^{w,k_{s}}}(\Omega)$
we simply replace $w$ with $\operatorname{inv}(w)$ in the bounds. By
Proposition 3.3 and since $w$ is invertible, we then have

$\displaystyle\|u\|^{2}_{\mathcal{H}_{k_{\operatorname{warp}}^{w,k_{s}}}(\Omega)}$
$\displaystyle=\min\left\\{\|\tilde{u}\|^{2}_{H_{k_{s}}(\Omega)}:u=\tilde{u}\circ
w\right\\}$
$\displaystyle=\|u\circ\operatorname{inv}(w)\|^{2}_{\mathcal{H}_{k_{s}}(\Omega)}$
$\displaystyle\leq(C_{\text{low}}(k_{s}))^{-2}\|u\circ\operatorname{inv}(w)\|_{H^{2}(\Omega)}$
$\displaystyle\leq(C_{\text{low}}(k_{s}))^{-2}\frac{1}{c^{\prime}}\max\left\\{1,\|\operatorname{inv}(w)\|_{W^{2,\infty}(\Omega)}^{2}\right\\}^{2}\|u\|_{H^{2}(\Omega)}^{2}.$

This finishes the proof, with

$\displaystyle C_{\text{low}}({k_{\operatorname{warp}}^{w,k_{s}}})$
$\displaystyle=C_{\text{low}}(k_{s})\sqrt{c^{\prime}}\max\left\\{1,\|\operatorname{inv}(w)\|_{W^{2,\infty}(\Omega)}^{2}\right\\}^{-1},$
$\displaystyle C_{\text{up}}({k_{\operatorname{warp}}^{w,k_{s}}})$
$\displaystyle=C_{\text{up}}(k_{s})\frac{1}{\sqrt{c}}\max\left\\{1,\|w\|_{W^{2,\infty}(\Omega)}^{2}\right\\}.$

$\square$

##### Proof:

[Of Theorem 3.8] We first consider the case $L=1$. To show
$\mathcal{H}_{k_{\operatorname{mix}}^{\sigma_{1},k_{1}}}(\Omega)\hookrightarrow
H^{\beta}(\Omega)$, let
$g\in\mathcal{H}_{k_{\operatorname{mix}}^{\sigma_{1},k_{1}}}(\Omega)$ be
arbitrary and note that by Proposition 3.4 there exists $u\in
H^{\beta}(\Omega)$ such that $g=\sigma_{1}u$ and
$\|g\|_{\mathcal{H}_{k_{\operatorname{mix}}^{\sigma_{1},k_{1}}}(\Omega)}=\|u\|_{\mathcal{H}_{k_{1}}(\Omega)}$.
Propositions 3.13 and 3.2 then give

$\displaystyle\|g\|_{H^{\beta}(\Omega)}$ $\displaystyle\leq
C_{\ref{prop:banachalg}}\|\sigma\|_{H^{\beta}(\Omega)}\|u\|_{H^{\beta}(\Omega)}$
$\displaystyle\leq
C_{\ref{prop:banachalg}}C_{\text{up}}(k_{1})\|\sigma\|_{H^{\beta}(\Omega)}\|u\|_{\mathcal{H}_{k_{1}}(\Omega)}$
$\displaystyle=C_{\ref{prop:banachalg}}C_{\text{up}}(k_{1})\|\sigma\|_{H^{\beta}(\Omega)}\|g\|_{\mathcal{H}_{k_{\operatorname{mix}}^{\sigma_{1},k_{1}}}(\Omega)}.$

To conversely show that
$H^{\beta}(\Omega)\hookrightarrow{\mathcal{H}_{k_{\operatorname{mix}}^{\sigma_{1},k_{1}}}(\Omega)}$,
let $u\in H^{\beta}(\Omega)$ be arbitrary. Propositions 3.4, 3.2, and 3.13
then give

$\displaystyle\|{u}\|_{\mathcal{H}_{k_{\operatorname{mix}}^{\sigma_{1},k_{1}}}(\Omega)}$
$\displaystyle=\|u/\sigma_{1}\|_{\mathcal{H}_{k_{1}}(\Omega)}$
$\displaystyle\leq(C_{\text{low}}(k_{1}))^{-1}\|{u}/\sigma_{1}\|_{H^{\beta}(\Omega)}$
$\displaystyle\leq
C_{\ref{prop:banachalg}}(C_{\text{low}}(k_{1}))^{-1}\|1/\sigma_{1}\|_{H^{\beta}(\Omega)}\|{u}\|_{H^{\beta}(\Omega)}.$

This finishes the proof for $L=1$. For the case $L=2$, without loss of
generality assume that $\nu_{1}<\nu_{2}$ and hence

$H^{\nu_{2}+d/2}(\Omega)\hookrightarrow
H^{\nu_{1}+d/2}(\Omega)=H^{\beta}(\Omega).$

We first show that
$\mathcal{H}_{k^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{2}}_{\operatorname{mix}}}(\Omega)\hookrightarrow
H^{\beta}(\Omega)$. Let
$g\in\mathcal{H}_{k^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{2}}_{\operatorname{mix}}(\Omega)}$
be arbitrary. By Propositions 3.6 and 3.4, we can write
$g=\sigma_{1}g_{1}+\sigma_{2}g_{2}$ for $g_{i}\in\mathcal{H}_{k_{i}}(\Omega)$
such that
$\|g\|^{2}_{\mathcal{H}_{k^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{2}}_{\operatorname{mix}}}(\Omega)}=\|g_{1}\|^{2}_{\mathcal{H}_{k_{1}}(\Omega)}+\|g_{2}\|_{\mathcal{H}_{k_{2}}(\Omega)}^{2}$.
Then the triangle inequality in $H^{\beta}(\Omega)$, together with
Propositions 3.13 and 3.4, gives

$\displaystyle\|g\|_{H^{\beta}(\Omega)}^{2}$
$\displaystyle=\|\sigma_{1}g_{1}+\sigma_{2}g_{2}\|_{H^{\beta}(\Omega)}^{2}$
$\displaystyle\leq
2C_{\ref{prop:banachalg}}^{2}\left(\|\sigma_{1}\|^{2}_{H^{\beta}(\Omega)}\|g_{1}\|^{2}_{H^{\beta}(\Omega)}+\|\sigma_{2}\|^{2}_{H^{\beta}(\Omega)}\|g_{2}\|_{H^{\beta}(\Omega)}^{2}\right)$
$\displaystyle\leq
2C_{\ref{prop:banachalg}}^{2}\max\\{C_{\text{up}}(k_{1})^{2}\|\sigma_{1}\|^{2}_{H^{\beta}(\Omega)},C_{\text{up}}(k_{2})^{2}\|\sigma_{2}\|^{2}_{H^{\beta}(\Omega)}\\}\left(\|g_{1}\|^{2}_{\mathcal{H}_{k_{1}}(\Omega)}+\|g_{2}\|_{\mathcal{H}_{k_{2}}(\Omega)}^{2}\right)$
$\displaystyle=2C_{\ref{prop:banachalg}}^{2}\max\\{C_{\text{up}}(k_{1})^{2}\|\sigma_{1}\|^{2}_{H^{\beta}(\Omega)},C_{\text{up}}(k_{2})^{2}\|\sigma_{2}\|^{2}_{H^{\beta}(\Omega)}\\}\|g\|^{2}_{\mathcal{H}_{k^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{2}}_{\operatorname{mix}}}(\Omega)}.$

To conversely show that
$H^{\beta}(\Omega)\hookrightarrow\mathcal{H}_{k^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{2}}_{\operatorname{mix}}}(\Omega)$,
let $g\in H^{\beta}(\Omega)$ be arbitrary. Then Propositions 3.6, 3.4, 3.2 and
3.13 give

$\displaystyle\|g\|^{2}_{\mathcal{H}_{k^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{2}}_{\operatorname{mix}}}(\Omega)}$
$\displaystyle=\min\\{\|g_{1}\|^{2}_{\mathcal{H}_{k_{1}}(\Omega)}+\|g_{2}\|^{2}_{\mathcal{H}_{k_{2}}(\Omega)}:g=g_{1}+g_{2},g_{i}\in\mathcal{H}_{k_{i}}(\Omega),i=1,2\\}$
$\displaystyle\leq\|g\|_{\mathcal{H}_{k_{1}}(\Omega)}^{2}$
$\displaystyle\leq(C_{\text{low}}(k_{1}))^{-2}\|g/\sigma_{1}\|^{2}_{H^{\beta}(\Omega)}$
$\displaystyle\leq
C_{\ref{prop:banachalg}}^{2}(C_{\text{low}}(k_{1}))^{-2}\|1/\sigma_{1}\|^{2}_{H^{\beta}}\|g\|^{2}_{H^{\beta}(\Omega)}.$

This finishes the proof for $L=2$. For $L\geq 3$, we can use an inductive
argument. Suppose we have already shown that
$\mathcal{H}_{k^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L-1}}_{\operatorname{mix}}}(\Omega)\cong
H^{\beta}(\Omega)$. Then notice we can write

$\displaystyle{k^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}}_{\operatorname{mix}}}(u,u^{\prime})$
$\displaystyle=\tilde{\sigma}(u)\tilde{\sigma}(u^{\prime}){k^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L-1}}_{\operatorname{mix}}}(u,u^{\prime})+\sigma_{L}(u)\sigma_{L}(u^{\prime})k_{L}(u,u^{\prime}),$

with $\tilde{\sigma}\equiv 1$. We then use the argument shown in this proof
for when $L=2$ to show that
$\mathcal{H}_{k^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}}_{\operatorname{mix}}}(\Omega)\cong
H^{\beta}(\Omega)$. In particular, we obtain

$\displaystyle
C_{\text{low}}(k_{\operatorname{mix}}^{{\\{\sigma_{\ell},k_{\ell}\\}}_{\ell=1}^{L}})$
$\displaystyle=C_{\ref{prop:banachalg}}^{-1}C_{\text{low}}(k_{j})\|1/\sigma_{j}\|_{H^{\beta}(\Omega)}^{-1},\qquad\nu_{j}=\min_{1\leq\ell\leq
L}\nu_{\ell},$ $\displaystyle
C_{\text{up}}(k_{\operatorname{mix}}^{{\\{\sigma_{\ell},k_{\ell}\\}}_{\ell=1}^{L}})$
$\displaystyle=\sqrt{2}C_{\ref{prop:banachalg}}\max_{{1\leq\ell\leq
L}}\\{C_{\text{up}}(k_{\ell})\|\sigma_{\ell}\|_{H^{\beta}(\Omega)}\\}.$

$\square$

##  4 Sample path regularity

In this section, we explore the path properties of a Gaussian process,
particularly focusing on the warping and mixture kernels introduced in section
2.2.2. Similar results hold for the convolution kernel, but are omitted for
brevity since they are not used in this work. These will play a crucial role
in the convergence analysis throughout the rest of this paper. We begin by
recalling known results on both the classic regularity of Mátern kernels and
the Sobolev regularity properties of a Gaussian process.

######  Proposition 4.1 ([38, Lemma 11])

Let $k=k_{\text{Mat}(\nu)}$. Then $k\in C^{2\nu}(\Omega\times\Omega)$ if $\nu$
is not an integer and $k\in C^{\alpha}(\Omega\times\Omega)$ for any
$\alpha<2\nu$ if $\nu$ is an integer.

######  Proposition 4.2 ([50, Theorem 2])

Let $k\in C^{p,p}(\Omega\times\Omega$) and $g\sim\mathcal{GP}(0,k)$. Then
$g\in H^{p}(\Omega)$ almost surely.

The remainder of this section establishes results regarding the classic
regularity properties of a Gaussian process. In particular, we present general
outcomes regarding the Hölder regularity properties based solely on the
condition given in Assumption 1 below. These results are detailed in Lemma
4.3, using the Sobolev embedding theorem, and in Lemma 4.8, using mean-square
properties. We provide more specific statements for the non-stationary kernels
discussed in section 2.2.2 in Corollaries 4.4 and 4.9.

While it may appear counterintuitive to employ two distinct methods for
presenting sample properties based on the dimension under consideration, the
Sobolev embedding theorems exhibit dimension dependence, while the mean-square
properties do not. As a result, we choose Sobolev methods when $d=1$ and rely
on mean-square properties when $d\geq 2$.

A Gaussian process $g\sim\mathcal{GP}(0,k)$ is said to be sample continuous if
$g$ is continuous on $\Omega$ almost surely, and sample $q$-partial
differentiable if $D^{q}g$ exists on $\Omega$ almost surely for
$q\in\mathbb{N}_{0}^{d}$.

###### Assumption 1

We have $k\in C^{p,p}(\Omega\times\Omega)$ for some $p>0$.

###### Lemma 4.3

Suppose that Assumption 1 is satisfied and $g\sim\mathcal{GP}(0,k)$. Then
$g\in C^{\alpha-d/2}(\Omega)$ almost surely, where $\alpha=p$ if $p$ is not an
integer and $\alpha<p$ otherwise.

##### Proof:

By Proposition 4.2, we see that $g\in H^{p}(\Omega)$ almost surely. The result
then follows from the Sobolev embedding theorem (see, e.g. [2, Theorem 4.12]).
$\square$

An application of Lemma 4.3 then gives the following result.

###### Corollary 4.4

Suppose that either:

  * •

$k=k_{\operatorname{warp}}^{w,k_{\text{Mat}(\nu)}}$ for $w\in C^{p}(\Omega)$
such that $p\leq\nu$, or

  * •

$k=k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\text{Mat}(\nu_{\ell})}\\}_{\ell=1}^{L}}$
for $\sigma_{\ell}\in C^{p_{\ell}}(\Omega)$ with $p=\min_{0\leq\ell\leq
L}{p_{\ell}}$ such that $p\leq\min_{0\leq\ell\leq L}{\nu_{\ell}}$.

Then for $g\sim\mathcal{GP}(0,k)$ we have that $g\in C^{\alpha-d/2}(\Omega)$
almost surely, where $\alpha=p$ if $p$ is not an integer and $\alpha<p$
otherwise.

##### Proof:

Assumption 1 is satisfied due to Proposition 4.1 and the chain rule or product
rule. $\square$

Sample regularity properties are crucial for our analytical framework, but
establishing the mean-square properties of a Gaussian process is often more
straightforward. Our objective is to infer conclusions about sample regularity
properties through an analysis of mean-square properties. We follow the
methodology outlined in [61, section 3], where Hölder sample path properties
for Gaussian processes with a Mátern kernel are established through the use of
mean-square properties.

We start by introducing the concepts of mean-square continuity and
differentiability. We say that $g\sim\mathcal{GP}(0,k)$ is

  * (i)

_mean-square continuous_ on $\Omega$, if for all $u\in\Omega$

$\mathbb{E}_{g\sim\mathcal{GP}(0,k)}[g(u)]<\infty\quad\text{ and
}\lim_{u^{\prime}\to
u}\mathbb{E}_{{}_{g\sim\mathcal{GP}(0,k)}}[|g(u)-g(u^{\prime})|^{2}]=0.$

  * (ii)

_mean-square $\bar{q}$-partial differentiable_ on $\Omega$ for
$\bar{q}\in\mathbb{N}^{d}_{0}$ with $|\bar{q}|=1$, if there exits a random
field $D^{\bar{q}}_{\operatorname{MS}}g$ with finite second moments such that
with $u_{\hat{q}}$ the unit vector in the direction of $\bar{q}$, we have for
all $u\in\Omega$

$\lim_{h\to
0}\mathbb{E}_{{}_{g\sim\mathcal{GP}(0,k)}}\left[\left|D^{\bar{q}}_{\operatorname{MS}}g-\frac{g(u)-g(u+u_{\bar{q}}h)}{h}\right|^{2}\right]=0.$

For general $q\in\mathbb{N}^{d}_{0}$, we denote $D^{q}_{\operatorname{MS}}g$
as the _mean-square $q$-partial derivative of $g$_, and this is defined
analogously. We refer to [61, section 3.1] for a more comprehensive definition
and discussion.

A modification of $g\sim{\mathcal{GP}}(0,k)$ is a random process $\tilde{g}$
such that for each $u\in\Omega$ we have
$\mathbb{P}_{g\sim\mathcal{GP}(0,k)}(g(u)=\tilde{g}(u))=1$. A Gaussian process
is completely defined by its finite dimensional distributions, thus allowing
us to work with a modification, as the laws of both $g$ and $\tilde{g}$ are
the same.

For a mean-square continuous Gaussian process $g\sim\mathcal{GP}(0,k)$ with
$k\in C^{|q|,|q|}(\Omega\times\Omega)$, we have that (see e.g. [61])

$D_{\operatorname{MS}}^{q}g\sim\mathcal{GP}\left(0,D^{q,q}k\right).$

We now prove the main result on the Hölder regularity of samples from a
Gaussian process through a sequence of lemmas. The final result is given in
Lemma 4.8.

###### Lemma 4.5

Suppose $k$ satisfies Assumption 1. Then $g\sim\mathcal{GP}(0,k)$ is mean-
square $i$-partial differentiable and
$D_{\operatorname{MS}}^{i}\mathcal{GP}(0,k)$ is mean-square continuous for
every $i\in\mathbb{N}_{0}^{d}$ with $0\leq|i|\leq\lceil p-1\rceil$.

##### Proof:

Our proof follows the proof of [61, Proposition 1]. Fix
$i\in\mathbb{N}_{0}^{d}$ with $0\leq|i|\leq\lceil p-1\rceil$. Mean-square
$i$-partial differentiability follows from [53, section 2] since $k\in
C^{p,p}(\Omega\times\Omega)$ by Assumption 1. Mean-square continuity of
$D_{\operatorname{MS}}^{i}\mathcal{GP}(0,k)$ follows since

$\mathbb{E}_{g\sim\mathcal{GP}(0,k)}[|D_{\operatorname{MS}}^{i}g(u)-D_{\operatorname{MS}}^{i}g(u^{\prime})|^{2}]=D^{i,i}(k(u,u)+k(u^{\prime},u^{\prime})-2k(u,u^{\prime}))\to
0\text{ as }u\to u^{\prime},$

since $k\in C^{p,p}(\Omega\times\Omega)$. $\square$

###### Lemma 4.6

Suppose $k$ satisfies Assumption 1. Then for all $i\in\mathbb{N}_{0}^{d}$ with
$0\leq|i|\leq\lceil p-1\rceil$, $D_{\operatorname{MS}}^{i}\mathcal{GP}(0,k)$
has a modification that is sample continuous.

##### Proof:

Since $D^{i,i}k\in C^{p-|i|,p-|i|}(\Omega\times\Omega)$ and since
$p-|i|>\lceil p-1\rceil-|i|>0$ we see that $D^{i,i}k$ is $p-|i|-$Hölder
continuous on $\Omega\times\Omega$ and so for
$u,u^{\prime},\tilde{u}\in\Omega$ we have

$\left|D^{i,i}k(u,u^{\prime})-D^{i,i}k(u,\tilde{u})\right|\leq
C(k)\|u^{\prime}-\tilde{u}\|_{2}^{p-|i|}.$

Hence, for any $u,u^{\prime}\in\Omega$

$D^{i,i}k(u,u)+D^{i,i}k(u^{\prime},u^{\prime})-2D^{i,i}k(u,u^{\prime})\leq
C(k)\|u-u^{\prime}\|_{2}^{p-|i|}.$

Following the proof of [61, Proposition 2, A.3] we see that this inequality,
via Kolmogrov’s Continuity Theorem (see e.g., [61, Theorem 3]), proves the
result. $\square$

###### Proposition 4.7

[61, Corollary 1] Let $k\in C^{\lceil p-1\rceil,\lceil
p-1\rceil}(\Omega\times\Omega)$, for some $p>0$. Suppose that for all
$i\in\mathbb{N}_{0}^{d}$ with $0\leq|i|\leq\lceil p-1\rceil$,
$g\sim\mathcal{GP}(0,k)$ has mean-square derivatives
$D^{i}_{\operatorname{MS}}g$, and $D^{i}_{\operatorname{MS}}g$ is mean-square
continuous and sample path continuous. Then $g$ has continuous sample
derivatives $D^{i}g$ which satisfy $D^{i}g=D_{\operatorname{MS}}^{i}g$ almost
surely, for all $0\leq|i|\leq\lceil p-1\rceil$.

###### Lemma 4.8

Suppose Assumption 1 is satisfied. Then there exists a modification
$\tilde{g}$ of $g\sim\mathcal{GP}(0,k)$ such that $\tilde{g}\in C^{\lceil
p-1\rceil}(\Omega)$ almost surely.

##### Proof:

By Assumption 1 we have $k\in C^{p,p}(\Omega\times\Omega)$ and so also $k\in
C^{\lceil p-1\rceil,\lceil p-1\rceil}(\Omega\times\Omega)$. By Lemma 4.5,
$g\sim\mathcal{GP}(0,k)$ is mean-square $\lceil p-1\rceil$-partial
differentiable and $D^{i}_{\operatorname{MS}}g$ is mean-square continuous for
$0\leq|i|\leq\lceil p-1\rceil$. By Lemma 4.6, we can work with a modification
$\tilde{g}$ of $g$ such that $D^{i}_{MS}\tilde{g}$ is sample continuous for
each $0\leq|i|\leq\lceil p-1\rceil$. The result then follows from Proposition
4.7. $\square$

Notice that in Lemma 4.3, the regularity of sample paths decreases as $d$
increases. In contrast, Lemma 4.8 is dimension-independent, and therefore will
often be the choice we will use in dimensions $d\geq 2$. Thus, in general, we
shall only use Lemma 4.3 when $d=1$. An application of Lemma 4.8 gives the
following result.

###### Corollary 4.9

Suppose that either:

  * •

$k=k_{\operatorname{warp}}^{w,k_{\text{Mat}(\nu)}}$ for $w\in C^{p}(\Omega)$
such that $p\leq\nu$ if $\nu$ is not an integer and $p<\nu$ otherwise, or

  * •

$k=k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\text{Mat}(\nu_{\ell})}\\}_{\ell=1}^{L}}$
for $\sigma_{\ell}\in C^{p_{\ell}}(\Omega)$ with $p=\min_{0\leq\ell\leq
L}{p_{\ell}}$ such that $p\leq\min_{0\leq\ell\leq L}{\nu_{\ell}}$.

Then there exists a modification $\tilde{g}$ of $g\sim\mathcal{GP}(0,k)$ such
that $\tilde{g}\in C^{\lceil p-1\rceil}(\Omega)$ almost surely.

##### Proof:

Assumption 1 is satisfied due to Proposition 4.1 and the chain rule or the
product rule. The result then follows from Lemma 4.8. $\square$

###### Remark 4.10

(Sobolev vs Hölder regularity) For $k\in C^{p,p}(\Omega\times\Omega)$, the
results in this section give sample path regularity $H^{p}(\Omega)$ and
$C^{\tilde{p}}(\Omega)$, for $\tilde{p}<p$. There are results in the
literature (see e.g. [14, Theorem 1(ii)]) which show that the Sobolev and
Hölder sample path regularity of some Gaussian processes is the same, i.e.
$\tilde{p}=p$. However, the proof of these results relies on properties of the
eigendecomposition of the covariance operator of the Gaussian process, and to
the best of our knowledge, there are no such results for general covariance
kernels. Hence, the Hölder regularity proved in Corollaries 4.4 and 4.9 may be
sub-optimal.

##  5 Error analysis of non-stationary GP regression

In this section we discuss the convergence of the GP emulator $f_{N}$, as
defined in (2) with fixed hyper-parameters $\theta$ and (8) with estimated
hyper-parameters $\widehat{\theta}_{N}$, to the true function $f$. In
particular, we show that the mean function $m_{N,\delta^{2}}^{f}$ converges to
the true function $f$ and, in the case of noise-free data, the covariance
function $k_{N,0}$ converges to 0. This type of convergence is required for
example when we want to use GP emulators as surrogate models in inverse
problems [55, 57, 20]. Convergence rates will be given in terms of the fill
distance $h_{U_{N},\Omega}$, which is the maximum distance any point in
$\Omega$ can be from $U_{N}$:

$h_{U_{N},\Omega}:=\sup_{u\in\Omega}\inf_{u_{n}\in U_{N}}\|u-u_{n}\|_{2}.$

This can be translated into convergence rates in terms of the number of
training points $N$ for specific point sets $U_{N}$, see, e.g. the discussion
in [57] and the references therein. The fastest possible rate of convergence
for $h_{U_{N},\Omega}$ is $N^{-1/d}$, and this is obtained for example for
equispaced points on a tensor grid.

We start by analysing the convergence of $m_{N,0}^{f}$ to $f$ (i.e. the noise-
free case) in sections 5.1 and 5.2, under different assumptions on the
function $f$ and kernel $k$ used in the GP prior (1). We analyse the
convergence of $m_{N,\delta^{2}}^{f}$ to $f$ (i.e. the noisy case) in section
5.3. We discuss the effect of estimating the hyper-parameters $\theta$ by
$\widehat{\theta}_{N}$ in an empirical Bayesian approach in section 5.4. In
section 5.5, we present various extensions to the analysis presented in the
previous sections, including convergence of $k_{N,0}$ to 0, faster convergence
of $m_{N,0}^{f}$ and $k_{N,0}$ for separable kernels, and convergence of
$m_{N,0}^{f}$, $m_{N,\delta^{2}}^{f}$ and $k_{N,0}$ in the case of randomly
distributed training points $U_{N}$.

###  5.1 Convergence of $m_{N,0}^{f}$ in the RKHS

In this section we are interested in the general setting where the native
space $\mathcal{H}_{k}$ corresponding to the prior covariance kernel $k$ is
not necessarily known explicitly. A benefit to this is that the results are
applicable to a wide variety of kernels, but this also means that showing that
$f$ satisfies the required assumption that it is an element of
$\mathcal{H}_{k}(\Omega)$ is often difficult to verify. A general convergence
result, based only on the general requirements given in Assumption 2, is
presented in Proposition 5.1 below.

###### Assumption 2

We have $k\in C^{2p}(\Omega\times\Omega)$, for some $p\in\mathbb{N}$, and
$f\in\mathcal{H}_{k}(\Omega)$.

######  Proposition 5.1 (Convergence in mean for $f\in\mathcal{H}_{k}(\Omega)$
[62, Theorem 11.13])

Suppose Assumption 2 holds with $p\in\mathbb{N}$. Then there exist constants
$h_{0}$ and $C$, independent of $f$, $k$ and $N$, such that

$\|D^{\alpha}f-D^{\alpha}m_{N,0}^{f}\|_{C^{0}(\Omega)}\leq
CC_{\ref{thm:Wend11.13}}(k)h^{p-|\alpha|}_{U_{N},\Omega}\|f\|_{\mathcal{H}_{k}(\Omega)},$

for all $\alpha\in\mathbb{N}^{d}_{0}$ with $|\alpha|<p$ and
$h_{U_{N},\Omega}\leq h_{0}$. The constant $C_{\ref{thm:Wend11.13}}$ is given
by

$C_{\ref{thm:Wend11.13}}(k):=\max_{\begin{subarray}{c}(\alpha_{1},\alpha_{2})\in\mathbb{N}^{2d}_{0}\\\
|\alpha_{1}|+|\alpha_{2}|=2p\end{subarray}}\max_{u_{1},u_{2}\in\Omega}\left|D^{\alpha_{1},\alpha_{2}}k(u_{1},u_{2})\right|.$

In Assumption 2, the regularity assumption on $k$ can often be easy to verify,
and this is done for specific examples of the warping, mixture and convolution
kernels from section 2.2.2 below. When the RKHS $\mathcal{H}_{k}(\Omega)$ is
not known explicitly, the regularity assumption on $f$ can be much more
difficult to verify, and we do not address this assumption further in this
work. When the RKHS is known explicitly, for example given by a Sobolev space
as in section 5.2, then faster convergence rates than those in Proposition 5.1
can sometimes be obtained by using the specific structure of the RKHS. Note
that even though the RKHS is not known explicitly, some general properties can
be deduced from Assumption 2. For example, we know that we have
$\mathcal{H}_{k}(\Omega)\hookrightarrow C^{p}(\Omega)$ by [62, Theorem 10.45].

In Lemmas 5.2, 5.3 and 5.4 below, we verify the regularity assumption on $k$
for special cases of the non-stationary kernels presented in section 2.2.2.
The resulting convergence rates are summarised in Corollary 5.5. For the
stationary kernel used in the construction of the non-stationary kernels, we
focus on the Gaussian kernel $k_{\text{Mat}(\infty)}$, for which the RKHS is
defined by rather technical conditions (see e.g. [19, Theorem 1]). Since the
Gaussian kernel is in $C^{\infty}(\Omega\times\Omega)$, the regularity of the
non-stationary kernel depends only on the regularity of the hyper-parameters
$w,\\{\sigma_{\ell}\\}_{\ell=1}^{L}$ and ${\lambda_{a}}$. If we use a
stationary kernel with finite regularity, then the regularity will depend both
on the hyper-parameters $w,\\{\sigma_{\ell}\\}_{\ell=1}^{L}$ and
${\lambda_{a}}$ and the stationary kernel. For instance, for the warping
kernel we obtain regularity $2p=\min\\{2q,2\tilde{q}\\}$ for $w\in
C^{2q}(\Omega)$ and $k_{s}\in C^{2\tilde{q}}(\Omega\times\Omega)$.

######  Lemma 5.2 (Warping kernel)

Let $\Omega\subset\mathbb{R}$. Suppose $(i)$ $k_{s}=k_{\text{Mat}(\infty)}$
and $(ii)$ $w\in C^{2p}(\Omega)$ for some $p\in\mathbb{N}$. Then
$k_{\operatorname{warp}}^{w,k_{s}}\in C^{2p}(\Omega\times\Omega)$, and

$C_{\ref{thm:Wend11.13}}(k_{\operatorname{warp}}^{w,k_{s}})=C_{0}\sigma^{2}\sqrt{(2p)!}B_{2p}\|w\|_{C^{2p}(\Omega)},$

where $C_{0}\leq 1.0866$ is a constant and $B_{2p}$ is the Bell number, which
represents the number of possible partitions of the set $\\{1,\ldots,2p\\}$.

######  Lemma 5.3 (Mixture kernel)

Let $\Omega\subset\mathbb{R}$. Suppose
$\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}$ are such that for all
$\ell=1,\ldots,L$ we have $(i)$ $k_{\ell}=k_{\text{Mat}(\infty)}$, and $(ii)$
$\sigma_{\ell}\in C^{2p}(\Omega)$ for some $p\in\mathbb{N}$. Then
$k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}}\in
C^{2p}(\Omega\times\Omega)$, and with $C_{0}\leq 1.0866$ is a constant,

$C_{\ref{thm:Wend11.13}}(k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}})=C_{0}2p2^{4p}\sqrt{(2p)!}\max_{1\leq\ell\leq
L}\|\sigma_{\ell}\|^{2}_{C^{2p}(\Omega)}.$

######  Lemma 5.4 (Convolution kernel)

Let $\Omega\subset\mathbb{R}$. Suppose $(i)$ $k_{i}=k_{\text{Mat}(\infty)}$,
$(ii)$ ${\lambda_{a}}(u)\geq c>0$ for all $u\in\Omega$, and $(iii)$
${\lambda_{a}}\in C^{2p}(\Omega)$ for some $p\in\mathbb{N}$. Then
$k_{\operatorname{conv}}^{{\lambda_{a}},k_{i}}\in C^{2p}(\Omega\times\Omega)$,
and

$\displaystyle
C_{\ref{thm:Wend11.13}}(k_{\operatorname{conv}}^{{\lambda_{a}},k_{i}})=2^{1/2}\sigma^{2}$
$\displaystyle\left(|5/4-2p|^{2p}\max\\{\|{\lambda_{a}}\|_{C^{0}(\Omega)},1\\}^{1/4-2p}\|{\lambda_{a}}\|_{C^{2p}(\Omega)}\right)^{2}$
$\displaystyle\Bigg{(}\frac{(4p)!}{2^{4p}(2p)!}2C_{0}\sqrt{(2p!)}\|{\lambda_{a}}\|^{2}_{C^{2p}(\Omega)}(2p)!(2c)^{-1-2p})\|{\lambda_{a}}\|_{C^{2p}(\Omega)}\Bigg{)}B_{2p}^{5}$

where $C_{0}\leq 1.0866$ is a constant and $B_{2p}$ is the Bell number, which
represents the number of possible partitions of the set $\\{1,\ldots,2p\\}$.

The proof of Lemmas 5.2, 5.3 and 5.4 can be found in the appendix. For brevity
we only consider the case $\Omega\subset\mathbb{R}$, but the general case can
be proved using similar techniques. An application of Proposition 5.1 then
gives the following result.

###### Corollary 5.5

Let
$k\in\\{k_{\operatorname{warp}}^{w,k_{s}},k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}},k_{\operatorname{conv}}^{{\lambda_{a}},k_{i}}\\}$.
Suppose $(i)$ the assumptions of Lemma 5.2, 5.3 or 5.4 hold with
$p\in\mathbb{N}$, respectively, and $(ii)$ $f\in\mathcal{H}_{k}(\Omega)$. Then
there exists a constant $h_{0}>0$ such that

$\|D^{\mu}f-D^{\mu}m_{N,0}^{f}\|_{C^{0}(\Omega)}\leq
C_{\ref{thm:Wend11.13}}(k)h^{p-|\mu|}_{U_{N},\Omega}\|f\|_{\mathcal{H}_{k}(\Omega)},$

for all $\mu\in\mathbb{N}^{d}_{0}$ with $|\mu|<p$ and $h_{U_{N},\Omega}\leq
h_{0}$.

###  5.2 Convergence of $m_{N,0}^{f}$ in a Sobolev space with noise-free data

The main drawback of the results in section 5.1 is the lack of transparency
concerning the assumption that $f$ needs to be an element of the RKHS of the
non-stationary kernel. In this section, we prove convergence results in the
special case that the RKHS is norm-equivalent to a Sobolev space, as
considered in section 3.2. A general convergence result, based only on the
general requirements given in Assumption 3, is presented in Proposition 5.6
below.

###### Assumption 3

For some $\beta\geq d/2$, we have $\mathcal{H}_{k}(\Omega)\cong
H^{\beta}(\Omega)$ with constants $C_{\text{low}}(k)$ and $C_{\text{up}}(k)$
such that

$C_{\text{low}}(k)\|\cdot\|_{\mathcal{H}_{k}(\Omega)}\leq\|\cdot\|_{H^{\beta}(\Omega)}\leq
C_{\text{up}}(k)\|\cdot\|_{\mathcal{H}_{k}(\Omega)},$

and $f\in H^{\beta}(\Omega)$.

######  Proposition 5.6 (Convergence in mean for $f\in H^{\beta}(\Omega)$ [57,
Theorem 3.5] )

Suppose Assumption 3 holds with $\beta\geq d/2$. Then there exist constants
$h_{0}$ and $C_{\ref{thm:Convergence_in_sob}}$, independent of $f$, $k$ and
$N$, such that

$\|f-m_{N,0}^{f}\|_{H^{\alpha}(\Omega)}\leq
C_{\ref{thm:Convergence_in_sob}}C_{\text{low}}(k)^{-1}C_{\text{up}}(k)\,h_{U_{N},\Omega}^{\beta-\alpha}\,\|f\|_{H^{\beta}(\Omega)},$

for all $\alpha\leq\beta$ and $h_{U_{N},\Omega}\leq h_{0}$.

The specific statement of Proposition 5.6, with the dependence on the
constants $C_{\text{low}}(k)$ and $C_{\text{up}}(k)$ made explicit, is due to
[57]. This will be important in section 5.4, when we estimate hyper-parameters
in $k$. However, the result was originally proved in [33].

An application of Proposition 5.6 then gives the following convergence result.

###### Corollary 5.7

Let
$k\in\\{k_{\operatorname{warp}}^{w,k_{s}},k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}}\\}$.
Suppose $(i)$ the assumptions of Theorem 3.7 or 3.8 hold with $\beta\geq d/2$,
respectively, and $(ii)$ $f\in H^{\beta}(\Omega)$. Then there exist constants
$h_{0}$ and $C_{\ref{thm:Convergence_in_sob}}$, independent of $f$, $k$ and
$N$, such that

$\|f-m_{N,0}^{f}\|_{H^{\alpha}(\Omega)}\leq
C_{\ref{thm:Convergence_in_sob}}C_{\text{low}}(k)^{-1}C_{\text{up}}(k)\,h_{U_{N},\Omega}^{\beta-\alpha}\,\|f\|_{H^{\beta}(\Omega)},$

for all $\alpha\leq\beta$ and $h_{U_{N},\Omega}\leq h_{0}$.

###  5.3 Convergence of $m_{N,\delta^{2}}^{f}$ in a Sobolev space with noisy
or noise-free data

In the case of noisy data, we have the following general convergence result.
The proof is given at the end of this section. Recall that our observations
take the form $y_{N}=\\{u_{n},f(u_{n})+\varepsilon_{n}\\}_{n=1}^{N}$.

###### Theorem 5.8

Suppose Assumption 3 holds with $\beta\geq d/2$. Then there exists a constant
$h_{0}$, independent of $f$, $k$, $\delta^{2}$ and $N$, such that for any
$\varepsilon=\\{\varepsilon_{n}\\}_{n=1}^{N}\in\mathbb{R}^{N}$ we have

$\displaystyle\|f-m_{N,\delta^{2}}^{f}\|_{H^{\alpha}(\Omega)}$
$\displaystyle\leq
h_{U_{N},\Omega}^{\beta-\alpha}\,\left(2C_{\text{low}}(k)^{-1}C_{\text{up}}(k)\|f\|_{H^{\beta}(\Omega)}+\delta^{-1}C_{\text{up}}(k)\|\varepsilon\|_{2}\right)$
$\displaystyle\qquad\qquad+\,h_{U_{N},\Omega}^{d/2-\alpha}\left(2\|\varepsilon\|_{2}+\delta
C_{\text{low}}(k)^{-1}C_{\text{up}}(k))\|f\|_{H^{\beta}(\Omega)}\right),$

for all $\alpha\leq\lfloor\beta\rfloor$ and $h_{U_{N},\Omega}\leq h_{0}$.

If further $f\in H^{\beta}(\overline{\Omega})$,
$0<k_{\mathrm{min}}:=\min_{u\in\overline{\Omega}}k(u,u)$ and $\varepsilon\sim
N(0,\delta^{2}\mathrm{I})$, then there exists a constant
$C_{\ref{thm:conv_sob_noisy}}$, independent of $f$, $k$, $\delta^{2}$ and $N$,
such that

$\displaystyle\mathbb{E}_{\varepsilon}\left[\|f-m_{N,\delta^{2}}^{f}\|_{H^{\alpha}(\Omega)}^{2}\right]^{1/2}$
$\displaystyle\leq
C_{\ref{thm:conv_sob_noisy}}\left(h_{U_{N},\Omega}^{\beta-\alpha}\,\left(C_{\text{low}}(k)^{-1}C_{\text{up}}(k)\|f\|_{H^{\beta}(\Omega)}+C_{\text{up}}(k)N^{1/2}\right)\right.$
$\displaystyle\qquad\left.+\,h_{U_{N},\Omega}^{d/2-\alpha}\left(C_{\text{low}}(k)^{-1+d/2\beta}\|f\|^{1-d/2\beta}_{H^{\beta}(\overline{\Omega})}+C_{\text{up}}(k)^{d/2\beta}\right)N^{d/4\beta}\right),$

for all $\alpha\leq\lfloor\beta\rfloor$, $h_{U_{N},\Omega}\leq h_{0}$ and
$N\geq N^{*}(k_{\mathrm{min}})$.

The first term in the error bounds in Theorem 5.8, involving
$h_{U_{N},\Omega}^{\beta-\alpha}C_{\text{low}}(k)^{-1}C_{\text{up}}(k)\|f\|_{H^{\beta}(\Omega)}$,
converges at the same rate as the error for noise-free data in Proposition
5.6. The other terms quantify the effect of the noise on convergence. Note
that the second bound does not derive from taking the second moment of the
first bound. Applying the equality
$\mathbb{E}_{\varepsilon}[\|\varepsilon\|_{2}^{2}]=\delta^{2}N$ to the first
bound would give a factor of $h_{U_{N},\Omega}^{d/2-\alpha}N^{1/2}$ in the
second bound, which is non-decreasing in $N$ for any $\alpha\geq 0$ since
$h_{U_{N},\Omega}=o(N^{-1/d})$. However, a refined proof technique allows to
conclude on the sharper bound with $N^{d/4\beta}$.

As already discussed in [64], we note that the first error bound in Theorem
5.8 applies to any noise $\varepsilon$, and hence also applies in the
misspecified setting where the assumption $\varepsilon\sim
N(0,\delta^{2}\mathrm{I})$ used to derive $m_{N,\delta^{2}}^{f}$ is not
satisfied. In particular, in the case of noise-free data with $\varepsilon=0$
we have

$\|f-m_{N,\delta^{2}}^{f}\|_{H^{\alpha}(\Omega)}\leq
h_{U_{N},\Omega}^{\beta-\alpha}\,\left(2C_{\text{low}}(k)^{-1}C_{\text{up}}(k)\|f\|_{H^{\beta}(\Omega)}\right)+\,h_{U_{N},\Omega}^{d/2-\alpha}\left(\delta
C_{\text{low}}(k)^{-1}C_{\text{up}}(k))\|f\|_{H^{\beta}(\Omega)}\right).$

Since $\beta\geq d/2$, the second term converges at a generally slower rate
than the first. This can be offset by choosing the noise $\delta$ small, and
in particular, we recover the rate $h_{U_{N},\Omega}^{\beta-\alpha}$ of
$\|f-m_{N,0}^{f}\|_{H^{\alpha}(\Omega)}$ if we choose
$\delta=O(h_{U_{N},\Omega}^{\beta-d/2})$.

An application of Theorem 5.8 gives the following convergence result.

###### Corollary 5.9

Let
$k\in\\{k_{\operatorname{warp}}^{w,k_{s}},k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}}\\}$.
Suppose $(i)$ the assumptions of Theorem 3.7 or 3.8 hold with $\beta\geq d/2$,
respectively, $(ii)$ $0<\sigma_{\mathrm{min}}\leq\sigma_{j}$ for
$\nu_{j}=\min_{1\leq\ell\leq L}\nu_{\ell}$ if
$k=k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}}$, and
$(iii)$ $f\in H^{\beta}(\overline{\Omega})$. Then Assumption 3 holds with
$\beta$, and the conclusions of Theorem 5.8 hold.

##### Proof:

[Of Theorem 5.8] We follow the same proof structure as in [64, Theorem 2 and
Lemma 17], making simplifications/improvements in some steps and keeping track
of how the constants appearing in the error estimates depend on the covariance
kernel $k$. For brevity, we outline the proof and only give detail in the
steps where modifications take place. Note that $(i)$ Assumptions 1-5 of [64]
are satisifed in the setting of this work, $(ii)$ Assumption 6 of [64] is
proved below, and $(iii)$ Assumption 7 of [64] can be relaxed to $f\in
H^{\beta}(\overline{\Omega})$ since we are assuming that the RKHS of $k$ is
$H^{\beta}(\Omega)$.

An application of [3, Theorem 3.2] (i.e. [64, Theorem 12]) gives, for a
constant $C_{1}$ independent of $f$, $k$ and $N$,

$\|f-m_{N,\delta^{2}}^{f}\|_{H^{\alpha}(\Omega)}\leq
C_{1}\left(h_{U_{N},\Omega}^{\beta-\alpha}\,\|f-m_{N,\delta^{2}}^{f}\|_{H^{\beta}(\Omega)}+h_{U_{N},\Omega}^{d/2-\alpha}\,\|f-m_{N,\delta^{2}}^{f}\|_{2,U_{N}}\right),$

for all $\alpha\leq\lfloor\beta\rfloor$ and $h_{U_{N},\Omega}\leq h_{0}$. Here
$\|g\|_{2,U_{N}}=\sqrt{\sum_{n=1}^{N}g(u_{n})^{2}}$ denotes the 2-norm of a
function $g\in H^{\beta}(\Omega)$ evaluated on the discrete set $U_{N}$.
Following the same proof technique as in [64, Theorem 2 and Lemma 17], we
obtain for the first term the bound

$\displaystyle\|f-m_{N,\delta^{2}}^{f}\|_{H^{\beta}(\Omega)}$
$\displaystyle\leq(1+C_{\text{low}}(k)^{-1}C_{\text{up}}(k))\|f\|_{H^{\beta}(\Omega)}+\delta^{-1}C_{\text{up}}(k)\|\varepsilon\|_{2}.$

For the second term, we have by [64, Lemma 17]

$\displaystyle\|f-m_{N,\delta^{2}}^{f}\|_{2,D_{N}}\leq
2\|\varepsilon\|_{2}+\delta
C_{\text{low}}(k)^{-1}C_{\text{up}}(k))\|f\|_{H^{\beta}(\Omega)}.$

This finishes the proof of the first claim.

For the second claim, For the second term, we use [58, Theorem 1] (i.e. [64,
Theorem 11]), which gives, for a constant $C_{2}$ independent of $f$, $k$ and
$N$,

$\mathbb{E}_{\varepsilon}\left[\|f-m_{N,\delta^{2}}^{f}\|_{2,U_{N}}^{2}\right]^{1/2}\leq
C_{2}N^{1/2}(\psi_{f,k}^{-1}(N)).$

Here, $\psi_{f,k}^{-1}$ is defined via

(11) $\displaystyle\psi_{f,k}^{-1}(N)$
$\displaystyle=\sup\\{\epsilon>0:\psi_{f,k}(\epsilon)\geq N\\},$
$\displaystyle\psi_{f,k}(\epsilon)$
$\displaystyle=\frac{\phi_{f,k}(\epsilon)}{\epsilon^{2}},$
$\displaystyle\phi_{f,k}(\epsilon)$
$\displaystyle=\inf_{\begin{subarray}{c}h\in\mathcal{H}_{k}(\overline{\Omega})\\\
\|h-f\|_{L^{\infty}(\overline{\Omega})}<\epsilon\end{subarray}}\frac{1}{2}\|h\|^{2}_{\mathcal{H}_{k}(\overline{\Omega})}-\log\mathbb{P}_{k}\left[\\{g:\|g\|_{L^{\infty}(\overline{\Omega})}\leq\epsilon\\}\right],$

where $\mathbb{P}_{k}$ denotes the measure of the prior $\mathcal{GP}(0,k)$.
We hence proceed to bound the concentration function $\phi_{f,k}(\epsilon)$.
For the first term, we use the simple bound $h=f$ to give

$\inf_{\begin{subarray}{c}h\in\mathcal{H}_{k}(\overline{\Omega})\\\
\|h-f\|_{L^{\infty}(\overline{\Omega})}<\epsilon\end{subarray}}\frac{1}{2}\|h\|^{2}_{\mathcal{H}_{k}(\overline{\Omega})}\leq\frac{1}{2}\|f\|^{2}_{\mathcal{H}_{k}(\overline{\Omega})}\leq\frac{1}{2}C_{\text{low}}(k)^{-2}\|f\|^{2}_{H^{\beta}(\overline{\Omega})},$

where by slight abuse of notation we use $C_{\text{low}}(k)$ to denote the
norm equivalence constant on $\Omega$ as well as $\overline{\Omega}$.

For the second term, we track the dependence on $k$ in the proof of [64, Lemma
24] to obtain the bound

$H(\mathcal{H}^{\prime}_{k}(\overline{\Omega}),\epsilon)\leq
C_{3}C_{\text{up}}(k)^{d/\beta}\epsilon^{-d/\beta},$

where $\mathcal{H}^{\prime}_{k}(\overline{\Omega})$ denotes the unit ball in
$\mathcal{H}_{k}(\overline{\Omega})$ and
$H(\mathcal{H}^{\prime}_{k}(\overline{\Omega}),\epsilon)$ denotes the metric
entropy, i.e. the logarithm of the $\epsilon$-covering number of
$\mathcal{H}^{\prime}_{k}(\overline{\Omega})$. The main observation to make is
that the unit ball in $\mathcal{H}_{k}(\overline{\Omega})$ is contained in the
ball of radius $C_{\text{up}}(k)$ in $H^{\beta}(\overline{\Omega})$. As in
[64, Equation (34)], we then have with
$\phi^{(2)}_{f,k}(\epsilon):=-\log\mathbb{P}_{k}\left[\\{g:\|g\|_{L^{\infty}(\overline{\Omega})}\leq\epsilon\\}\right]$

$\phi^{(2)}_{f,k}(\epsilon)=\log
2+C_{3}C_{\text{up}}(k)^{d/\beta}\epsilon^{-d/\beta}8^{d/2\beta}\phi^{(2)}_{f,k}(\epsilon/2)^{d/2\beta}.$

Now observe that, using the lower bound
$\int_{c}^{\infty}\exp(-t^{2}/2)\mathrm{d}t\geq\frac{c}{c^{2}+1}\exp(-c^{2}/2)$
(see e.g. [18, Equation (2.23)], we have for $u\in\overline{\Omega}$ s.t.
$k(u,u)=\min_{u\in\overline{\Omega}}k(u,u)=k_{\mathrm{min}}$ the bound

$\mathbb{P}_{k}\left[\\{g:\|g\|_{L^{\infty}(\overline{\Omega})}\leq
c\\}\right]\leq\mathbb{P}_{k}\left[\\{g:|g(u)|\leq c\\}\right]\leq
1-\frac{2ck_{\mathrm{min}}}{\sqrt{2\pi}(c^{2}+k_{\mathrm{min}}^{2})}\exp(-\frac{c^{2}}{2k_{\mathrm{min}}^{2}}).$

This means that [64, Assumption 6] is satisfied for any $N\in\mathbb{N}$ with
any $c>0$ and
$\gamma_{N}\equiv\gamma(k,c)=-\log(1-\frac{2ck_{\mathrm{min}}}{\sqrt{2\pi}(c^{2}+k_{\mathrm{min}}^{2})}\exp(-\frac{c^{2}}{2k_{\mathrm{min}}^{2}}))$.

For $\epsilon/2\leq\min\\{c,(\log
2)^{-1}C_{3}8^{d/2\beta}\gamma(k,c)^{d/2\beta})^{\beta/d}\\}=:\epsilon^{*}$
for some $0<c<1$, we then have

$\phi^{(2)}_{f,k}(\epsilon)=C_{3}(C_{\text{up}}(k)^{d/\beta}+1)\epsilon^{-d/\beta}8^{d/2\beta}\phi^{(2)}_{f,k}(\epsilon/2)^{d/2\beta}.$

Taking logarithms and following the iterative procedure as in the proof of
[64, Theorem 2], we arrive at the estimate

$\phi^{(2)}_{f,k}(\epsilon)\leq
C_{4}C_{\text{up}}(k)^{2d/(2\beta-d)}\epsilon^{-2d/(2\beta-d)},$

for all $\epsilon\leq\epsilon^{*}$ and a constant $C_{4}$ independent of $f$,
$k$ and $N$. This gives

$\displaystyle\phi_{f,k}(\epsilon)\epsilon^{-2}$
$\displaystyle\leq\frac{1}{2}C_{\text{low}}(k)^{-2}\|f\|^{2}_{H^{\beta}(\overline{\Omega})}\epsilon^{-2}+C_{4}C_{\text{up}}(k)^{2d/(2\beta-d)}\epsilon^{-2-2d/(2\beta-d)},$
$\displaystyle\leq\left(\frac{1}{2}C_{\text{low}}(k)^{-2}\|f\|^{2}_{H^{\beta}(\overline{\Omega})}+C_{4}C_{\text{up}}(k)^{2d/(2\beta-d)}\right)\epsilon^{-2-2d/(2\beta-d)},$

for all $\epsilon\leq\epsilon^{*}$. For sufficiently large $N\geq
N^{*}(\epsilon^{*})$, we then have

(12)
$\psi_{f,k}^{-1}(N)\leq\left(\frac{1}{2}C_{\text{low}}(k)^{-2}\|f\|^{2}_{H^{\beta}(\overline{\Omega})}+C_{4}C_{\text{up}}(k)^{2d/(2\beta-d)}\right)^{1/2-d/(4\beta)}N^{-1/2+d/4\beta}.$

Combining all steps together and using the equality
$\mathbb{E}_{\varepsilon}[\|\varepsilon\|_{2}^{2}]=\delta^{2}N$, we have for
$\alpha\leq\lfloor\beta\rfloor$, $h_{U_{N},\Omega}^{\beta-\alpha}\leq h_{0}$
and $N\geq N^{*}(k_{\mathrm{min}}$,

$\displaystyle\mathbb{E}_{\varepsilon}\left[\|f-m_{N,\delta^{2}}^{f}\|_{H^{\alpha}(\Omega)}^{2}\right]^{1/2}\leq
C_{1}\left(h_{U_{N},\Omega}^{\beta-\alpha}\,\left((1+C_{\text{low}}(k)^{-1}C_{\text{up}}(k))\|f\|_{H^{\beta}(\Omega)}+C_{\text{up}}(k)N^{1/2}\right)\right.$
$\displaystyle\qquad\qquad\left.+\,h_{U_{N},\Omega}^{d/2-\alpha}\,C_{2}\left(\frac{1}{2}C_{\text{low}}(k)^{-2}\|f\|^{2}_{H^{\beta}(\overline{\Omega})}+C_{4}C_{\text{up}}(k)^{2d/(2\beta-d)}\right)^{-1/2+d/(4\beta)}N^{d/4\beta}\right),$

which completes the proof. $\square$

###  5.4 Hyper-parameters estimated by empirical Bayes’

In this section we study the effect of the estimation of hyper-parameters on
the error bounds presented in sections 5.1, 5.2 and 5.3. Recall from section
2.3 that within the framework of empirical Bayes’, we compute a point estimate
$\widehat{\theta}_{N}$ of the hyper-parameters $\theta$ using the data
$y_{N}$, and then plug this into the prior (1) to be used in the standard GP
regression prior-to-posterior update. For ease of presentation, we will
restrict our attention to the case where the stationary kernel is a Matérn
kernel $k_{\text{Mat}(\nu)}$. For the non-stationary kernels considered in
this work, we are then interested in estimating the following hyper-
parameters:

  * •

For $k_{\operatorname{warp}}^{w,k_{s}}$, we estimate
$\theta=\\{w,\sigma^{2}\\}$ and fix $\\{\lambda,\nu\\}$.

  * •

For $k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}}$, we
estimate
$\theta=\\{L,\\{\sigma_{\ell}\\}_{\ell=1}^{L},\\{\lambda_{\ell}\\}_{\ell=1}^{L}\\}$
and fix $\\{\sigma^{2},\\{\nu_{\ell}\\}_{\ell=1}^{L}\\}.$

  * •

For $k_{\operatorname{conv}}^{{\lambda_{a}},k_{i}}$, we estimate
$\theta=\\{{\lambda_{a}},\sigma^{2}\\}$ and fix $\\{\lambda,\nu\\}$.

For the warping and convolution kernels, the non-stationary functions $w$ and
${\lambda_{a}}$, respectively, are designed to capture length scales. The
parameter $\lambda$ in the Matérn kernel is hence superfluous, and can be kept
fixed (e.g. at 1). Similarly, the marginal variance parameter $\sigma^{2}$ in
the Matérn kernel can be kept fixed (e.g. at 1) for the mixture kernel.

In all cases, it is possible to extend the results to the smoothness parameter
$\nu$ being estimated as well. If the target function $f$ is in the RKHS of
the non-stationary kernel for all possible values of $\widehat{\nu}_{N}$, the
convergence results presented here hold with straightforward modifications. If
$f$ on the other hand has limited regularity, and is not necessarily an
element of the RKHS for all $\widehat{\nu}_{N}$, then the error bounds require
a qualitative change as presented in [57, Theorem 3.5]. Since this would
significantly increase required notation and technicality, and since the focus
of this work is on modelling different length scales, we omitted estimating
$\nu$ in this work for ease of presentation.

####  5.4.1 Convergence of $m_{N,0}^{f}$ in the RKHS with noise-free data

If Assumption 2 holds for $k(\cdot,\cdot;\widehat{\theta}_{N})$, for all
$N\in\mathbb{N}$, then Proposition 5.1 immediately gives the error bound

$\|D^{\mu}f-D^{\mu}m_{N,0}^{f}(\widehat{\theta}_{N})\|_{C^{0}(\Omega)}\leq
C_{\ref{thm:Wend11.13}}(k(\cdot,\cdot;\widehat{\theta}_{N}))\;h^{p-|\mu|}_{U_{N},\Omega}\;\|f\|_{\mathcal{H}_{k(\cdot,\cdot;\widehat{\theta}_{N})}(\Omega)}.$

The convergence rates we obtain hence depend on the behaviour of
$C_{\ref{thm:Wend11.13}}(k(\cdot,\cdot;\widehat{\theta}_{N}))$ and
$\|f\|_{\mathcal{H}_{k(\cdot,\cdot;\widehat{\theta}_{N})}(\Omega)}$. If these
quantities can be bounded uniformly in $N$, we obtain the same convergence
rates as in the case of fixed hyper-parameters, see Corollary 5.11. More
generally, we obtain a convergence rate of $h^{p-|\mu|}_{U_{N},\Omega}\,N^{r}$
if
$C_{\ref{thm:Wend11.13}}(k(\cdot,\cdot;\widehat{\theta}_{N}))\|f\|_{\mathcal{H}_{k(\cdot,\cdot;\widehat{\theta}_{N})}(\Omega)}\leq
CN^{r}$ for some $r\in\mathbb{R}$ for our sequence of estimated hyper-
parameters $\\{\widehat{\theta}_{N}\\}_{N=1}^{\infty}$.

If the RKHS $\mathcal{H}_{k(\cdot,\cdot;\widehat{\theta}_{N})}(\Omega)$ is not
known explicitly, it may be difficult to quantify the behaviour of
$\|f\|_{\mathcal{H}_{k(\cdot,\cdot;\widehat{\theta}_{N})}(\Omega)}$. The
constant $C_{\ref{thm:Wend11.13}}(k(\cdot,\cdot;\widehat{\theta}_{N}))$ is
more tractable, and can be bounded uniformly in $N$ under the assumptions
given in Theorem 5.10 below.

###### Theorem 5.10

Let $p\in\mathbb{N}$ and $\Omega\subset\mathbb{R}$. Suppose either

  * (i)

$k=k_{\operatorname{warp}}^{w,k_{s}}$, $k_{s}=k_{\text{Mat}(\infty)}$ and
$\\{\|\widehat{w}_{N}\|_{C^{2p}(\Omega)},\widehat{\sigma}^{2}_{N}\\}$ lie in a
compact subset of $(0,\infty)^{2}$,

  * (ii)

$k=k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}}$,
$k_{\ell}=k_{\text{Mat}(\infty)}$ for $1\leq\ell\leq\infty$, and
$\\{\widehat{L}_{N},\max_{1\leq\ell\leq\widehat{L}_{N}}\|\widehat{\sigma}_{\ell,N}\|_{C^{2p}(\Omega)}\\}$
lie in a compact subset of $\mathbb{N}\times(0,\infty)$, or

  * (iii)

$k=k_{\operatorname{conv}}^{{\lambda_{a}},k_{i}}$,
$k_{i}=k_{\text{Mat}(\infty)}$ and
$\\{\|\widehat{\lambda_{a,N}}\|_{C^{2p}(\Omega)},\inf_{u\in\Omega}\widehat{\lambda_{a,N}}(u),\widehat{\sigma}^{2}_{N}\\}$
lie in a compact subset of $(0,\infty)^{3}$.

Then there exists a constant $C_{\ref{thm: conditions on hyperparam native
space norm}}$, independent of $N$, such that

$\sup_{\widehat{\theta}_{N}}C_{\ref{thm:Wend11.13}}(k(\cdot,\cdot;\widehat{\theta}_{N}))\leq
C_{\ref{thm: conditions on hyperparam native space norm}}.$

##### Proof:

For $(i)$ the condition on $\widehat{\sigma}^{2}_{N}$ follows directly from
[57, Theorem 3.5], and the condition on $\widehat{w}_{N}$ follows from Lemma
5.2. For $(ii)$, the conditions follow from Lemma 5.3. For $(iii)$, the
conditions on $\widehat{\sigma}^{2}_{N}$ follow directly from [57, Theorem
3.5], and the condition on $\widehat{\lambda}_{a,N}$ and
$\inf_{u\in\Omega}{\lambda_{a,N}}(u)$ follows from Lemma 5.4. $\square$

Note that there are no conditions on
$\\{\widehat{\lambda}_{\ell,N}\\}_{\ell=1}^{L}$ in Theorem 5.10 $(ii)$. This
is because the only terms that depend on
$\\{\widehat{\lambda}_{\ell,N}\\}_{\ell=1}^{L}$ are inside the exponential and
this can always be bounded by $1$, i.e.
$\exp(-\|u-u^{\prime}\|_{2}^{2}/\widehat{\lambda}_{\ell,N})\leq 1$ for all
$u,u^{\prime}\in\Omega$ and $0<\widehat{\lambda}_{\ell,N}<\infty$, see (16) in
the appendix.

###### Corollary 5.11

Suppose the assumptions Theorem 5.10 hold. Then for
$f\in\mathcal{H}_{k(\cdot,\cdot;\widehat{\theta}_{N})}(\Omega)$ there exists a
constant $h_{0}$, independent of $f$ and $N$, such that

$\|D^{\mu}f-D^{\mu}m_{N,0}^{f}(\widehat{\theta}_{N})\|_{C^{0}(\Omega)}\leq
C_{\ref{thm: conditions on hyperparam native space
norm}}\;h^{p-|\mu|}_{U_{N},\Omega}\;\sup_{\widehat{\theta}_{N}}\|f\|_{\mathcal{H}_{k(\cdot,\cdot;\widehat{\theta}_{N})}(\Omega)},$

for all $\mu\in\mathbb{N}_{0}$ with $\mu<p$ and $h_{U_{N},\Omega}\leq h_{0}$.

####  5.4.2 Convergence of $m_{N,\delta^{2}}^{f}$ in a Sobolev space with
noisy or noise-free data

If Assumption 3 holds for $k(\cdot,\cdot;\widehat{\theta}_{N})$, for all
$N\in\mathbb{N}$, Proposition 5.6 and Theorem 5.8 similarly give the error
bounds

$\displaystyle\|f-m_{N,0}^{f}(\widehat{\theta}_{N})\|_{H^{\alpha}(\Omega)}$
$\displaystyle\leq
C_{\ref{thm:Convergence_in_sob}}C_{\text{low}}(k(\cdot,\cdot;\widehat{\theta}_{N}))^{-1}C_{\text{up}}(k(\cdot,\cdot;\widehat{\theta}_{N}))\,h_{U_{N},\Omega}^{\beta-\alpha}\,\|f\|_{H^{\beta}(\Omega)},$
$\displaystyle\|f-m_{N,\delta^{2}}^{f}(\widehat{\theta}_{N})\|_{H^{\alpha}(\Omega)}$
$\displaystyle\leq
h_{U_{N},\Omega}^{\beta-\alpha}\,\left(2C_{\text{low}}(k(\widehat{\theta}_{N}))^{-1}C_{\text{up}}(k(\widehat{\theta}_{N}))\|f\|_{H^{\beta}(\Omega)}+\delta^{-1}C_{\text{up}}(k(\widehat{\theta}_{N}))\|\varepsilon\|_{2}\right)$
$\displaystyle\qquad\qquad+\,h_{U_{N},\Omega}^{d/2-\alpha}\left(2\|\varepsilon\|_{2}+\delta
C_{\text{low}}(k(\widehat{\theta}_{N}))^{-1}C_{\text{up}}(k(\widehat{\theta}_{N}))\|f\|_{H^{\beta}(\Omega)}\right).$

Note that in contrast to section 5.4.1, the norm on $f$ appearing here is
independent of any of the estimated hyper-parameters. The obtained convergence
rates hence depend only on the behaviour of the constants
$C_{\text{low}}(k(\cdot,\cdot;\widehat{\theta}_{N}))^{-1}$ and
$C_{\text{up}}(k(\cdot,\cdot;\widehat{\theta}_{N}))$. In particular, Theorem
5.12 below gives conditions under which moments of these constant can be
bounded uniformly in $N$.

###### Theorem 5.12

Let $\beta\geq d/2$. Suppose either

  * (i)

$k=k_{\operatorname{warp}}^{w,k_{s}}$, $k_{s}=k_{\text{Mat}(\beta-d/2)}$ and
$\\{\|\widehat{w}_{N}\|_{C^{\beta}(\Omega)},\|\operatorname{inv}(\widehat{w}_{N})\|_{C^{\beta}(\Omega)},\min_{u\in\Omega}\widehat{w}_{N}^{\prime}(u),\min_{u\in\Omega}\operatorname{inv}(\widehat{w}_{N})^{\prime}(u),\widehat{\sigma}^{2}_{N}\\}$
lie in a compact subset of $(0,\infty)^{5}$, or

  * (ii)

$k=k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}}$,
$k_{\ell}=k_{\text{Mat}(\nu_{\ell})}$ for $1\leq\ell\leq\infty$, and
$\\{\widehat{L}_{N},\max_{1\leq\ell\leq\widehat{L}_{N}}\\{\|\widehat{\sigma}_{\ell,N}\|_{H^{\beta}(\Omega)},\|1/\widehat{\sigma}_{\ell,N}\|_{H^{\beta}(\Omega)}\\}\\}$
lie in a compact subset of $\mathbb{N}\times(0,\infty)^{2}$.

Then for any $a,b>0$ there exists a constant $C_{\ref{lem:cond on estimates
for Sob bnds}}(a,b)$, independent of $N$, such that

$\sup_{\widehat{\theta}_{N}}C_{\text{low}}(k(\cdot,\cdot;\widehat{\theta}_{N}))^{-a}C_{\text{up}}(k(\cdot,\cdot;\widehat{\theta}_{N}))^{b}\leq
C_{\ref{lem:cond on estimates for Sob bnds}}(a,b).$

##### Proof:

For $(i)$, the condition on $\widehat{\sigma}^{2}_{N}$ follows directly from
[57, Theorem 3.5] and the conditions on $\widehat{w}_{N}$ follow from Theorem
3.7. For $(ii)$, this follows from Theorem 3.8. $\square$

###### Corollary 5.13

Suppose the assumptions of Theorem 5.12 and Corollary 5.7 hold. Then

$\|f-m_{N,0}^{f}(\widehat{\theta}_{N})\|_{H^{\alpha}(\Omega)}\leq
C_{\ref{thm:Convergence_in_sob}}C_{\ref{lem:cond on estimates for Sob
bnds}}(1,1)\,h_{U_{N},\Omega}^{\beta-\alpha}\,\|f\|_{H^{\beta}(\Omega)},$

for all $\alpha\leq\beta$ and $h_{U_{N},\Omega}\leq h_{0}$.

###### Corollary 5.14

Suppose the assumptions of Theorem 5.12 and Corollary 5.9 hold, for
$\sigma_{\mathrm{min}}$ independent of $N$. Then the conclusions of Theorem
5.8 hold.

####  5.4.3 Computing hyper-parameter estimates

The results presented in sections 5.4.1 and 5.4.2 do not make any assumptions
on how the estimated hyper-parameters $\widehat{\theta}_{N}$ are computed. The
methodology to do this efficiently is an important research area of its own,
but is not the focus of this work and so we only briefly discuss this here.

Estimation of the stationary parameters $\\{\lambda,\sigma^{2}\\}$ in the
Matérn kernel, including questions of identifiability, has been studied
extensively (see, e.g. the discussion and references in [57]), and is
frequently done using maximum marginal likelihood estimation or cross
validation. The non-stationary parameters $w$,
$\\{\sigma_{\ell}\\}_{\ell=1}^{L}$ and ${\lambda_{a}}$ are functions, and so
typically a suitable parametrisation is chosen before estimation. In deep
kernel GPs [63], which have garnered considerable interest recently, the
warping function $w$ takes the form of a deep neural network. In [60], the
authors discuss efficient estimation of the mixture coefficients
$\\{\sigma_{\ell}\\}_{\ell=1}^{L}$ and the number of components $L$ using
methods akin to cross validation. A piece-wise linear model is chosen for the
length scale ${\lambda_{a}}$ in [16].

###  5.5 Extensions of error analysis

In this section, we provide several extensions of the error analysis provided
in sections 5.1-5.4. We apply established proof techniques in each context,
and, for brevity, we will omit detailed derivations and refer the reader
instead to the provided references.

#### Convergence of $k_{N,0}$

In the case of noise-free data $y_{N}$, the posterior variance $k_{N,0}$
satisfies the equality in Proposition 5.15 below, allowing us to transfer
convergence results on $m_{N,0}^{f}$ to convergence results on $k_{N,0}$ (see,
e.g. [57]). Given the design points $U_{N}$, we introduce the mapping
$m_{N,0}^{(\cdot)}:\mathcal{H}_{k}(\Omega)\rightarrow\mathcal{H}_{k}(\Omega)$,
which is built on the definition of the posterior mean $m_{N,0}^{f}$ given in
(3). With $g(D_{N}):=[g(u^{1});\dots;g(u^{N})]\in\mathbb{R}^{N}$, we define

(13) $\displaystyle g(u)\mapsto
m_{N,0}^{g}(u)=k(u,D_{N})^{T}K(D_{N},D_{N})^{-1}g(D_{N}).$

######  Proposition 5.15 (e.g. [57, Proposition 3.2])

Suppose $k_{N,0}$ and $m_{N}^{(\cdot)}$ are given by (4) with $\delta^{2}=0$
and (13), respectively. Then for any $u\in\Omega$ we have

$k_{N,0}(u,u)^{\frac{1}{2}}=\sup_{\begin{subarray}{c}{g\in\mathcal{H}_{k}(\Omega)\;\text{s.t.
}}\\\ \|g\|_{\mathcal{H}_{k}(\Omega)}=1\end{subarray}}|g(u)-m^{g}_{N}(u)|.$

An application of Proposition 5.15 then gives us the following convergence
results on $k_{N,0}$. In slight abuse of notation, let use denote by
$k(\cdot;\widehat{\theta}_{N})$ the point-wise posterior variance at $u$, i.e.
$k(u;\widehat{\theta}_{N}):=k(u,u;\widehat{\theta}_{N})$ for $u\in\Omega$.

###### Corollary 5.16

Suppose the assumptions of Proposition 5.1 hold. Then there exists a constant
$h_{0}$, independent of $f$ and $N$, such that

$\|k_{N,0}(\widehat{\theta}_{N})^{1/2}\|_{C^{0}(\Omega)}\leq
C_{\ref{thm:Wend11.13}}(k(\cdot,\cdot;\widehat{\theta}_{N}))h^{p}_{U_{N},\Omega},$

for all $h_{U_{N},\Omega}\leq h_{0}$.

##### Proof:

This follows directly from Propositions 5.1 and 5.15. $\square$

###### Corollary 5.17

Suppose the assumptions of Proposition 5.6 hold. Then there exist constants
$h_{0}$ and $C_{\ref{cor:var_sobolev}}$, independent of $f$, $k$ and $N$, such
that for any $\epsilon>0$

$\|k_{N,0}(\widehat{\theta}_{N})^{1/2}\|_{L^{2}(\Omega)}\leq
C_{\ref{cor:var_sobolev}}C_{\text{low}}(k(\cdot,\cdot;\widehat{\theta}_{N}))^{-1}C_{\text{up}}(k(\cdot,\cdot;\widehat{\theta}_{N}))^{2}\,h_{U_{N},\Omega}^{\beta-d/2-\epsilon},$

for $h_{U_{N},\Omega}\leq h_{0}$.

##### Proof:

This follows from Propositions 5.6 and 5.15, following the same proof
technique as in [57, Theorem 3.8]. $\square$

Under the conditions of Theorems 5.10 and 5.12, respectively, the constants in
the above error bounds can again be bounded uniformly in $N$.

#### Separable covariance kernels

All error bounds derived so far have been expressed in terms of the fill
distance $h_{U_{N},\Omega}$. This can be translated into convergence rates in
$N$ for specific sets of design points $U_{N}$. For example, the fill distance
decays at rate $N^{-1/d}$ for uniform grids. This is the fastest decay rate
that can be obtained by any $U_{N}$, and so there is a strong dimension
dependence in the error bounds presented so far.

For $f\in H^{\beta}(\Omega)$, the rate $N^{-\beta/d}$ given in Proposition 5.6
for $\|f-m_{N,0}^{f}\|_{L^{2}(\Omega)}$ is optimal (see, e.g. the discussion
in [57]) and so cannot be improved. However, by imposing further restrictions
on $f$, we are able to improve the error bounds. This additionally requires a
specific structure on $\Omega$, $U_{N}$ and $k$. In particular, we require
that the kernel $k$ is separable, i.e. that it can be written as the product
of one-dimensional kernels
$k(u,u^{\prime})=\prod_{j=1}^{d}k^{j}(u_{(j)},u^{\prime}_{(j)})$, and that the
domain $\Omega$ is of tensor product structure
$\Omega=\prod_{j=1}^{d}\Omega_{j}$. The training points $U_{N}$ are chosen as
a sparse grid [7], defined as the set of points

$H(q,{{d}}):=\bigcup_{\|i\|_{1}=q}X_{1}^{(i_{1})}\times\cdots
X_{{d}}^{(i_{{d}})},$

where $i\in\mathbb{N}_{0}^{{d}}$, $q\geq{{d}}$, and for $1\leq j\leq{{d}}$ we
choose a sequence $X_{j}^{(p)}:=\\{x_{j,1}^{(p)},\dots,x_{j,m_{p}}^{(p)}\\}$,
$p\in\mathbb{N}$, of nested sets of points in $\Omega_{j}$. We then have the
following result.

###### Corollary 5.18

Suppose $(i)$ $\Omega=\prod_{j=1}^{{d}}\Omega_{j}$, $(ii)$ $U_{N}$ is chosen
as the Smolyak sparse grid $H(q,{{d}})$ with $h_{X_{j}^{(p)},\Omega_{j}}\leq
C_{h}m_{p}^{-r_{h}}$, $(iii)$
$k^{j}\in\\{k_{\operatorname{warp}}^{w,k_{s}},k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}}\\}$,
satisfying the assumptions of Theorem 3.7 or 3.8, respectively, and $(iv)$
$f\in
H^{\\{\beta_{j}\\}}_{\otimes^{{d}}}(\Omega):=\otimes_{j=1}^{{d}}H^{\beta_{j}}(\Omega_{j})$.
Then there exist constants $C_{\ref{cor:conv_sep}}$ and
$C_{\ref{cor:conv_sep}}^{\prime}$, independent of $f$, $k$ and $N$, such that

$\displaystyle\|f-m_{N,0}^{f}(\widehat{\theta}_{N})\|_{H^{\\{\alpha_{j}\\}}_{\otimes^{{d}}}(\Omega)}$
$\displaystyle\leq
C_{\ref{cor:conv_sep}}\left(\prod_{j=1}^{d}C_{\text{low}}^{-1}(k^{j}(\widehat{\theta}_{N}))C_{\text{up}}(k^{j}(\widehat{\theta}_{N}))\right)N^{-\gamma}(\log
N)^{(1+\gamma)(d-1)}\|f\|_{H^{\\{\beta_{j}\\}}_{\otimes^{{d}}}(\Omega)},$
$\displaystyle\|k_{N}(\widehat{\theta}_{N})^{1/2}\|_{L^{2}(\Omega)}$
$\displaystyle\leq
C_{\ref{cor:conv_sep}}^{\prime}\left(\prod_{j=1}^{d}C_{\text{low}}^{-1}(k^{j}(\widehat{\theta}_{N}))C_{\text{up}}(k^{j}(\widehat{\theta}_{N}))^{2}\right)N^{-\gamma+\frac{1}{2}+\epsilon}(\log
N)^{(\frac{3}{2}+\epsilon+\gamma)(d-1)}.$

for any $\epsilon>0$ and $\\{\alpha_{j}\\}\leq\\{\beta_{j}\\}$, where
$\gamma=\min_{1\leq j\leq{{d}}}r_{h}(\beta_{j}-\alpha_{j})$.

##### Proof:

This follows directly from [57, Theorems 3.11 and 3.12], together with
Theorems 3.7 and 3.8. $\square$

Notice that the rate $\gamma$ in Corollary 5.18 is the same rate in $N$ that
we would obtain with Proposition 5.6 in the case $d=1$, provided
$h_{U_{N},\Omega}\leq CN^{-1/d}$ in Proposition 5.6 and
$h_{X_{j}^{(p)},\Omega_{j}}\leq C_{h}m_{p}^{-1}$ in Corollary 5.18 (i.e. the
fill distance decays at the optimal rate in both cases). Up to a log factor,
the convergence rates in Corollary 5.18 are independent of $d$. Under the
conditions of Theorem 5.12 on $k^{j}$, the constants in the above error bounds
can again be bounded uniformly in $N$. Extending Corollary 5.18 to the case of
noisy data is beyond the scope of this work.

#### Random training points $U_{N}$

When the training points $U_{N}$ are chosen randomly, error bounds in terms of
the fill distance for fixed $U_{N}$ are not always meaningful. Corollary 5.19
below gives a result in expectation over $U_{N}$ in the case of noise-free
data.

###### Corollary 5.19

Suppose $U_{N}$ are sampled i.i.d. from a measure $\mu$ with density ${\rho}$
satisfying ${\rho(u)}\geq{\rho_{\mathrm{min}}}>0$ for all
$u\in\overline{\Omega}$, and the assumptions of Corollary 5.13 hold
$\mu$-almost surely. Then there exist constants $C_{\ref{cor:conv_random}}$
and $C_{\ref{cor:conv_random}}^{\prime}$, independent of $f$ and $N$, such
that for all $\alpha\leq\beta$ and $\epsilon>0$ we have

$\displaystyle\mathbb{E}_{\mu}\left[\|f-m_{N,0}^{f}(\widehat{\theta}_{N})\|_{H^{\alpha}(\Omega)}\right]$
$\displaystyle\leq
C_{\ref{cor:conv_random}}\;N^{-\frac{\beta-\alpha}{d}+\epsilon}\,\|f\|_{H^{\beta}(\Omega)},$
$\displaystyle\mathbb{E}_{\mu}\left[\|k_{N,0}^{1/2}(\widehat{\theta}_{N})\|_{L^{2}(\Omega)}\right]$
$\displaystyle\leq
C_{\ref{cor:conv_random}}^{\prime}\;N^{-\frac{\beta}{d}+\frac{1}{2d}+\epsilon}\,\|f\|_{H^{\beta}(\Omega)}.$

##### Proof:

The claim follows from [20, Theorem 3], with $U_{N}=\emptyset$ and
$\pi^{y}(u)\equiv 1$, as well as Proposition 5.15 and the same proof technique
as in [57, Theorem 3.8]. $\square$

The proof of this result is mainly based on the observation that under the
conditions of Theorem 5.12 (holding almost surely in $\mu$) the error bound in
Proposition 5.6 depends on $U_{N}$ only through the fill distance
$h_{U_{N},\Omega}$, together with bounds on
$\mathbb{E}_{\mu}[h_{U_{N},\Omega}^{\beta-\alpha}]$ and
$\mathbb{P}_{\mu}[h_{U_{N},\Omega}>h_{0}]$ proved in [39, 57]. Note that the
rates obtained in Corollary 5.19 are almost optimal, in the sense that the
fill distance $h_{U_{N},\Omega}$ in Corollary 5.13 in the best case decays at
rate $N^{-1/d}$.

For noisy data, we consider the following result for random training points
$U_{N}$. Recall the definition of $\psi_{f,k}^{-1}$ from (11).

###### Corollary 5.20

Suppose $U_{N}$ are sampled i.i.d. from a measure $\mu$ with density ${\rho}$
satisfying ${\rho(u)}\geq{\rho_{\mathrm{min}}}>0$ for all
$u\in\overline{\Omega}$, and the assumptions of Corollary 5.9 hold
$\mu$-almost surely. For $\hat{\beta}=\beta$ for $\beta$ not an integer and
$\hat{\beta}<\beta$ otherwise, define $\alpha(1)=\beta-1/2$ for
$k_{\operatorname{warp}}^{w,k_{s}}$, and define $\alpha(1)=\hat{\beta}-1/2$
and $\alpha(d)=\lceil\hat{\beta}-d/2-1\rceil$, for $d\geq 2$, for
$k=k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}}$. Let
$\psi_{f,k}^{-1}$ be as in (11).

Then for $f\in H^{\beta}(\overline{\Omega})$ there exists a constant
$C_{\ref{cor:conv_random w/ noise}}(k)$, independent of $f$ and $N$, such that
for sufficiently large $N\geq N^{*}(\sigma_{\min})$

$\displaystyle\mathbb{E}_{\varepsilon}\left[\mathbb{E}_{\mu}\left[\|f-m_{N,\delta^{2}}^{f}(\widehat{\theta}_{N})\|^{2}_{L^{2}(\Omega)}\right]\right]^{1/2}$
$\displaystyle\leq\left(\frac{1}{2}C_{\text{low}}(k(\widehat{\theta}_{N}))^{-2}\|f\|^{2}_{H^{\beta}(\overline{\Omega})}+C_{4}C_{\text{up}}(k(\widehat{\theta}_{N}))^{2d/(2\beta-d)}\right)^{1/2-d/(4\beta)}$
$\displaystyle\hskip 128.0374ptC_{\ref{cor:conv_random w/
noise}}(k(\widehat{\theta}_{N}))N^{-r}$

where

$\quad r=\begin{cases}1/2-d/4\beta,\quad&\text{ if }\psi_{f,k}^{-1}(N)\leq
N^{-d/(4\alpha(d)+2d)},\\\ (\alpha(d)+d)/\beta&\text{ if
}\psi_{f,k}^{-1}(N)\geq N^{-d/(4\alpha(d)+2d))}.\end{cases}$

##### Proof:

We break the proof into two cases depending on the kernel $k$.

Case 1, Warping: We have $w\in C^{\beta}(\Omega)$ for $\beta=\nu+1/2$ by
assumption (see Theorem 3.7). Since then $w\in C^{\nu}(\Omega)$, we use
Corollary 4.4 with $p=\nu$ to see that for any
$g\sim\mathcal{GP}(0,k_{\operatorname{warp}}^{w,k_{s}})$ we have $g\in
C^{\alpha(1)}(\Omega)$ almost surely.

Case 2, Mixture: $\sigma_{\ell}\in H^{\beta}(\Omega)$ for $\beta=\nu_{j}+d/2$
with $\nu_{j}=\min_{1\leq\ell\leq L}\nu_{\ell}$ by assumption (see Theorem
3.8). By the Sobolev embedding theorem (see, e.g. [2, Theorem 4.12])
$\sigma_{\ell}\in C^{\hat{\beta}-d/2}(\Omega)$. For $d=1$ (respectively $d\geq
2$) we use Corollary 4.4 (respectively Corollary 4.9) with
$p=\hat{\beta}-d/2\leq{\nu_{j}}$ to see that for any
$g\sim\mathcal{GP}(0,k_{\operatorname{mix}}^{\\{\sigma_{\ell},k_{\ell}\\}_{\ell=1}^{L}})$
we have $g\in C^{\alpha(d)}(\Omega)$ almost surely.

In both cases, using the upper bound for $\psi_{f,k}^{-1}(N)$ in (12) and [58,
Theorem 2], the claim then follows with $C_{\ref{cor:conv_random w/
noise}}(k(\widehat{\theta}_{N}))=\rho_{\min}^{-1}C$, where $C$ is the constant
appearing in [58, Theorem 2]. $\square$ Note that Corollary 5.20 is stated for
a fixed estimate for the hyper-parameters $\widehat{\theta}_{N}$. It is
possible to extend this to a bound uniform in the estimated hyper-parameters,
as in Corollary 5.19, by tracking the dependence on $k$ of the constant $C$
appearing in [58, Theorem 2]. For brevity, we have omitted this analysis.

##  6 Error analysis of deep (and wide) GP regression

In this section, we employ the developed theory to prove convergence rates for
DGP regression. Recall from section 2.3 that in this set-up, we have a deep
Gaussian process $f^{D}$ of depth $D$, constructed using a sequence of
stochastic processes that are conditionally Gaussian [11, 13]:

$\displaystyle f^{0}$ $\displaystyle\sim\mathcal{GP}(0,k_{0}(u,u^{\prime})),$
$\displaystyle f^{n}|f^{n-1}$
$\displaystyle\sim\mathcal{GP}(0,k_{n}(u,u^{\prime};f^{n-1})),\qquad
n=1,\dots,D.$

The initial layer $f^{0}$ is defined as a stationary GP, for example with
Matérn covariance function $k_{0}=k_{\text{Mat}(\nu_{0})}$. In subsequent
layers, the covariance kernel is constructed using a sample from the previous
layer. In this way, we can design a DGP based on non-stationary kernels as in
section 2.2.2, where the non-stationary function(s) $w$,
$\\{\sigma_{\ell}\\}_{\ell=1}^{L}$ or $\lambda_{a}$ rely on a sample from the
previous layer. In particular, we consider the following constructions:

  * •

deep GPs: $k_{0}=k_{\text{Mat}(\nu_{0})}$, and
$k_{n}\in\\{k_{\operatorname{warp}}^{f^{n-1},k_{\text{Mat}(\nu_{n})}},k_{\operatorname{mix}}^{F(f^{n-1}),k_{\text{Mat}(\nu_{n})}}\\}$
with $F$ strictly positive.

All other hyper-parameters appearing in $\\{k_{n}\\}_{n\in\mathbb{N}}$ are
kept fixed. Note that this is a more general construction than the one
considered in [13], since we are allowing the smoothness $\nu_{n}$ to vary
with $n$. In the mixture kernel
$k_{\operatorname{mix}}^{F(f^{n-1}),k_{\text{Mat}(\nu_{n})}}$ we only have one
mixture component in this set-up, and $F(f^{n-1})^{2}$ hence defines a non-
stationary marginal variance for $f^{n}$.

Alternatively, we would like to place prior distributions on the multiple
mixture components $\\{\sigma_{\ell}\\}_{\ell=1}^{L}$. To achieve this, we
construct a DGP width $L$ and depth $D=1$ so that

$\displaystyle f^{0,\ell}$
$\displaystyle\sim\mathcal{GP}(0,k_{0}(u,u^{\prime}))\quad\text{ for
}\ell=1,\ldots,L,\quad i.i.d.,$ $\displaystyle
f^{1,L}|\\{f^{0,\ell}\\}_{\ell=1}^{L}$
$\displaystyle\sim\mathcal{GP}(0,k_{1}(u,u^{\prime};\\{f^{0,\ell}\\}_{\ell=1}^{L})).$

For clarity and to distinguish this construction from the DGP construction
discussed above, we refer to this construction as a wide GP (WGP). We then
consider the following construction:

  * •

wide GPs: $k_{0}=k_{\text{Mat}(\nu_{0})}$, and
$k_{1}=k_{\operatorname{mix}}^{\\{F(f^{0,\ell}),k_{\text{Mat}(\nu_{1})}\\}_{\ell=1}^{L}}$
with $F$ strictly positive.

Again, all other hyper-parameters appearing in $k_{1}$ are kept fixed.

The remainder of this section is devoted to proving convergence rates for DGP
and WGP regression, i.e. regression where the prior is chosen as a DGP or WGP.
Conditioning the DGP or WGP priors on noise-free data becomes very technical
and in general requires working with disintegrations (see e.g. [9] and the
references therein). We hence only consider posterior distributions
conditioned on noisy data
$y_{N}=\\{u_{n},f(u_{n})+\varepsilon_{n}\\}_{n=1}^{N}$, where
$\varepsilon_{n}\sim N(0,\delta^{2})$ i.i.d. with $\delta>0$. We denote the
resulting posteriors on the final layers $f^{D}$ and $f^{1,L}$ by
$\mu^{y_{N},\delta^{2},D}_{D}$ and $\mu^{y_{N},\delta^{2},L}_{1}$,
respectively. By Bayes’ Theorem (see e.g. [54]), the posteriors on the final
layer satisfy

$\frac{\mathrm{d}\mu^{y_{N},\delta^{2},D}_{D}}{\mathrm{d}\mu_{\mathrm{prior}}^{D}}(g)\propto\exp\left(-\frac{1}{2\delta^{2}}\|y_{N}-g(U_{N})\|_{2}^{2}\right),\qquad\frac{\mathrm{d}\mu^{y_{N},\delta^{2},L}_{1}}{\mathrm{d}\mu_{\mathrm{prior}}^{1,L}}(g)\propto\exp\left(-\frac{1}{2\delta^{2}}\|y_{N}-g(U_{N})\|_{2}^{2}\right),$

with $\mu_{\mathrm{prior}}^{D}$ and $\mu_{\mathrm{prior}}^{1,L}$ denoting the
corresponding priors. We denote the resulting posterior distributions on the
other layers by $\mu^{y_{N},\delta^{2},D}_{n}$, for $n=0,\dots,D-1$, and
$\mu^{y_{N},\delta^{2},L}_{0}$, respectively. In the case of noise-free data,
the posteriors derived in this way can be seen as a computable approximation.

In our analysis, we will be interested in the convergence of the means
$\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D}}[f^{D}]$ and
$\mathbb{E}_{\mu^{y_{N},\delta^{2},L}_{1}}[f^{1,L}]$ to the true function $f$.
In contrast to standard GP regression, the posterior distributions can no
longer be found explicitly, and there is no analytical expression for these
means. In practice, one typically uses Markov chain Monte Carlo (MCMC) methods
for sampling from the posterior [13, 29], and the means are then computed as
sample averages.

###  6.1 Convergence of (constrained) DGPs in a Sobolev space

We have the following general convergence result for DGP regression.

###### Assumption 4

Let $D\geq 1$. Suppose

  * (i)

$f$ and $k_{D}(\cdot,\cdot;f^{D-1})$ satisfy Assumption 3 with $\beta>d/2$ for
almost all $f^{D-1}$,

  * (ii)

for some $\tilde{\beta}\in\mathbb{N}_{0}$ and all $q\in\mathbb{N}_{0}^{d}$
with $|q|\leq\tilde{\beta}$, there exists a random variable $Z$ such that
$\sup_{u\in\Omega}|D^{q}(f-m_{N,\delta^{2}}^{f}(u;f^{D-1}))|\leq Z$ almost
surely with $\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D-1}}(Z)<\infty$,

  * (iii)

$\max\left\\{\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D-1}}\left[C_{\text{low}}(k_{D})^{-1}C_{\text{up}}(k_{D})\right],\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D-1}}\left[C_{\text{up}}(k_{D})\right]\right\\}\leq
C_{\ref{assump:dgp_sobolev}}(y_{N})$, with
$C_{\ref{assump:dgp_sobolev}}(y_{N})\leq
C_{\ref{assump:dgp_sobolev}}^{a}<\infty$ if $\varepsilon=0$,
$\mathbb{E}_{\varepsilon}[C_{\ref{assump:dgp_sobolev}}(y_{N})^{2}]^{1/2}\leq
C_{\ref{assump:dgp_sobolev}}^{b}<\infty$ if $\varepsilon\sim
N(0,\delta^{2}\mathrm{I})$, and $C_{\ref{assump:dgp_sobolev}}^{a}$ and
$C_{\ref{assump:dgp_sobolev}}^{b}$ independent of $N$ and $U_{N}$.

###### Theorem 6.1

Suppose Assumption 4 holds. Then, for any $\varepsilon\in\mathbb{R}^{N}$, we
have

$\displaystyle\|f-\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D}}[f^{D}]\|_{H^{\alpha}(\Omega)}$
$\displaystyle\leq
h_{U_{N},\Omega}^{\beta-\alpha}\left(2\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D-1}}\left[C_{\text{low}}(k_{D})^{-1}C_{\text{up}}(k_{D})\right]\|f\|_{H^{\beta}(\Omega)}+\delta^{-1}\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D-1}}\left[C_{\text{up}}(k_{D})\right]\|\varepsilon\|_{2}\right)$
$\displaystyle\qquad+\,h_{U_{N},\Omega}^{d/2-\alpha}\left(2\|\varepsilon\|_{2}+\delta\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D-1}}\left[C_{\text{low}}(k_{D})^{-1}C_{\text{up}}(k_{D})\right]\|f\|_{H^{\beta}(\Omega)}\right),$

for all $\alpha\leq\min\\{\tilde{\beta},\lfloor\beta\rfloor\\}$ and
$h_{U_{N},\Omega}\leq h_{0}$, where $h_{0}$ is as in Theorem 5.8. Hence,

  1. (i)

for $\varepsilon=0$ and $\delta=\delta_{N}\leq
C_{\delta}h_{U_{N},\Omega}^{\beta-d/2}$, we have for all
$\alpha\leq\min\\{\tilde{\beta},\lfloor\beta\rfloor\\}$ and
$h_{U_{N},\Omega}\leq h_{0}$

$\|f-\mathbb{E}_{\mu^{y_{N},\delta_{N}^{2},D}_{D}}[f^{D}]\|_{H^{\alpha}(\Omega)}\leq(2+C_{\delta})C_{\ref{assump:dgp_sobolev}}^{a}h_{U_{N},\Omega}^{\beta-\alpha}\|f\|_{H^{\beta}(\Omega)},$

  2. (ii)

for $\varepsilon\sim N(0,\delta^{2}\mathrm{I})$ and
$\delta=\delta_{N}=C_{\delta}N^{-1/2}$, we have for all
$\alpha\leq\min\\{\tilde{\beta},\lfloor\beta\rfloor\\}$ and
$h_{U_{N},\Omega}\leq h_{0}$

$\displaystyle\mathbb{E}_{\varepsilon}\left[\|f-\mathbb{E}_{\mu^{y_{N},\delta_{N}^{2},D}_{D}}[f^{D}]\|_{H^{\alpha}(\Omega)}\right]\leq\left(2\|f\|_{H^{\beta}(\Omega)}+1\right)C_{\ref{assump:dgp_sobolev}}^{b}\,h_{U_{N},\Omega}^{\beta-\alpha}\,N^{1/2}+\left(2+C_{\ref{assump:dgp_sobolev}}^{b}\|f\|_{H^{\beta}(\Omega)}\right)C_{\delta}\,h_{U_{N},\Omega}^{d/2-\alpha}.$

##### Proof:

We first consider the case $\alpha\in\mathbb{N}_{0}$. Using the definition of
conditional distributions we see that

$\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D}}[f^{D}]=\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D-1}}[\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D|D-1}}[f^{D}]]=\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D-1}}[m_{N,\delta^{2}}^{f}(f^{D-1})],$

where the final equality uses that the conditional distribution
$\mu^{y_{N},\delta^{2},D}_{D|D-1}$ on $f^{D}|f^{D-1},y_{N}$ is Gaussian with
mean $m_{N,\delta^{2}}^{f}({f^{D-1}})$ as defined in (3). The linearity of
expectation then gives

$\displaystyle\|f-\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D}}[f^{D}]\|_{H^{\alpha}(\Omega)}$
$\displaystyle=\|\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D-1}}[f-m_{N,\delta^{2}}^{f}(\cdot;{f^{D-1}})]\|_{H^{\alpha}(\Omega)}.$

To use Theorem 5.8, we now wish to exchange the order of
$\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D-1}}$ and
$\|\cdot\|_{H^{\alpha}(\Omega)}$. Using Assumption 4$(ii)$, we can conclude
that we can do so by the dominated convergence theorem (see, e.g. [2, Theorem
1.50])) . Together with the convexity of $\|\cdot\|_{L^{2}(\Omega)}$ and
Jensen’s inequality, we then have

$\displaystyle\left\|\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D-1}}[f-m_{N,\delta^{2}}^{f}(f^{D-1})]\right\|_{H^{\alpha}(\Omega)}$
$\displaystyle=\sum_{|p|\leq\alpha}\left\|\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D-1}}[D^{p}(f-m_{N,\delta^{2}}^{f}(f^{D-1}))]\right\|_{L^{2}(\Omega)}$
$\displaystyle\leq\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D-1}}\left[\left\|f-m_{N,\delta^{2}}^{f}(f^{D-1})\right\|_{H^{\alpha}(\Omega)}\right].$

Combining the previous bound with Assumption 4$(i)$ and Theorem 5.8, we have

$\displaystyle\|f-\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D}}[f^{D}]\|_{H^{\alpha}(\Omega)}$
$\displaystyle\leq
h_{U_{N},\Omega}^{\beta-\alpha}\left(2\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D-1}}\left[C_{\text{low}}(k_{D})^{-1}C_{\text{up}}(k_{D})\right]\|f\|_{H^{\beta}(\Omega)}+\delta^{-1}\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D-1}}\left[C_{\text{up}}(k_{D})\right]\|\varepsilon\|_{2}\right)$
$\displaystyle\qquad+\,h_{U_{N},\Omega}^{d/2-\alpha}\left(2\|\varepsilon\|_{2}+\delta\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D-1}}\left[C_{\text{low}}(k_{D})^{-1}C_{\text{up}}(k_{D})\right]\|f\|_{H^{\beta}(\Omega)}\right),$

which proves the first claim. The other claims then follow using Assumption
4$(iii)$, Hölder’s inequality with exponents $p=q=1/2$ on the term
$\mathbb{E}_{\varepsilon}\left[\mathbb{E}_{\mu^{y_{N},\delta^{2},D}_{D-1}}\left[C_{\text{up}}(k_{D})\right]\|\varepsilon\|_{2}\right]$,
and the inequality
$\mathbb{E}_{\varepsilon}\left[\|\varepsilon\|_{2}\right]\leq\mathbb{E}_{\varepsilon}\left[\|\varepsilon\|^{2}_{2}\right]^{1/2}=\delta$.

The case of non-integer $\alpha$ follows by interpolation theory in Sobolev
spaces (see e.g. [6]), by considering the linear operator $T$ defined by
$Tf=f-m_{N}^{f}$. The first bound applies to $\lfloor\alpha\rfloor$ and
$\lceil\alpha\rceil$, and hence also to
$\alpha=\lfloor\alpha\rfloor+\theta(\lceil\alpha\rceil-\lfloor\alpha\rfloor)$,
for $0<\theta<1$. $\square$

Similar to Theorem 5.8, the first claim in Theorem 6.1 does not make any
assumptions on the noise $\varepsilon$. The second claim, given in part $(i)$,
shows that with noise-free data and the noise parameter $\delta$ in the
approximate posterior $\mu^{y_{N},\delta^{2},D}_{D}$ chosen appropriately, we
obtain optimal convergence rates $N^{-(\beta+\alpha)/d}$ of the posterior mean
also for DGP regression. The third claim, given in part $(ii)$, shows that
convergence of the posterior mean is guaranteed also in the case of noisy
data, provided that the noise level $\delta$ goes to zero at rate $N^{-1/2}$.
Note however, that the rates are slower than in the case of noise-free data.
For a constant noise level $\delta$, independent of $N$, the error bound in
Theorem 6.1 unfortunately does allow us to conclude on the convergence of the
posterior mean as $N\rightarrow\infty$, cf the discussion after Theorem 5.8.

We now show that Assumption 4 holds for the mixture and warping kernel DGP
constructions, given the conditions specified in Lemmas 6.2 and 6.3. In order
to ensure that Assumption 4 $(ii)$ and $(iii)$ are satisfied, we impose the
condition that sample paths from the $D-1$-th layer are restricted to