and

t1Co-funded by the European Union (ERC, BigBayesUQ, project number:
101041064). Views and opinions expressed are however those of the author(s)
only and do not necessarily reflect those of the European Union or the
European Research Council. Neither the European Union nor the granting
authority can be held responsible for them.

# Uncertainty quantification for sparse spectral variational approximations in
Gaussian process regression

Dennis<EMAIL_ADDRESS>[ Department of Mathematics, Vrije
Universiteit Amsterdam, The Netherlands
Botond<EMAIL_ADDRESS>[ Department of Decision
Sciences and BIDSA, Bocconi University, Italy
Harry van<EMAIL_ADDRESS>[ Department of Mathematics,
Vrije Universiteit Amsterdam, The Netherlands

###### Abstract

We investigate the frequentist properties of the variational sparse Gaussian
Process regression model. In the theoretical analysis we focus on the
variational approach with spectral features as inducing variables. We derive
guarantees and limitations for the frequentist coverage of the resulting
variational credible sets. We also derive sufficient and necessary lower
bounds for the number of inducing variables required to achieve minimax
posterior contraction rates. The implications of these results are
demonstrated for different choices of priors. In a numerical analysis we
consider a wider range of inducing variable methods and observe similar
phenomena beyond the scope of our theoretical findings.

62G20,

62G05, 62G08, 62G15,

Variational Bayes,

Uncertainty quantification,

Nonparametric regression,

Inducing variables method,

Bayesian asymptotics,

###### keywords:

[class=MSC]

###### keywords:

##  1 Introduction

One of the key challenges in Bayesian statistics is to approximate intractable
or computationally infeasible posterior distributions. This problem is
becoming even more pronounced in applications where the amount of available
information is rapidly growing, further increasing the complexity of the
posterior. Variational methods provide a convenient way to overcome such
computational issues in Bayesian statistics. The variational approximation
starts by selecting an appropriate class of distributions, referred to as the
variational class. Then the complex posterior distribution is approximated by
its projection on this variational class with respect to the Kullback-Leibler
divergence. The challenge in choosing the variational class is twofold:
firstly, the variational posterior should reduce computational complexity and
if possible increase the interpretability of the distribution; secondly, for
meaningful inference it is crucial that the variational posterior has good
statistical properties. For an overview of variational Bayes methods we refer
to the review article [5].

Although variational Bayes approximations are routinely used in practice, up
to recently they were considered black box procedures with very limited
theoretical underpinning. In the last few years the asymptotic properties of
the variational posterior have been increasingly investigated. Abstract
results were derived for posterior contraction rates and applied to various
high-dimensional and nonparametric models, considering typically mean-field
variational classes, see for instance [36, 35, 1, 19, 20]. However, almost all
of these results focus on the recovery of the underlying true parameter of
interest and do not address the quality of uncertainty quantification.

In fact, one of the main appeals and strengths of the Bayesian paradigm is
that it provides a probabilistic solution to the statistical problem. The
posterior distribution can be used to quantify the remaining uncertainty about
the parameter of interest. In practice this uncertainty is usually visualised
by plotting credible regions. These are subsets of the parameter space with
prescribed posterior probability (typically $95\%$). In parametric models the
celebrated Bernstein-von Mises theorem [16, 28] provides asymptotic
frequentist coverage guarantees for credible sets under mild assumptions,
meaning that credible sets can be interpreted as frequentist confidence sets.
In high-dimensional and non-parametric frameworks such a strong guarantee does
not automatically hold in general, e.g. [8]. Nevertheless, by now we have a
relatively good understanding of how to tune the prior to achieve asymptotic
confidence guarantees, see for instance [15, 7, 25, 21].

Despite its importance, so far hardly any results are available on the
frequentist reliability of the variational Bayesian uncertainty
quantification, where the real credible sets are replaced by approximate
credible sets derived from the variational posterior. In fact, many of the
available results are rather negative, showing that (mean-field) variational
methods are often over-confident in the sense that they substantially
underestimate the uncertainty of the procedure, see for instance [4, 5] for
some standard examples. There are only a few positive results available. In
[33] a variational version of the Bernstein-von Mises theorem was derived in
parametric models, while in [11] a correction was proposed using linear
response methods to recover the original posterior covariance structure.
However, both of these results consider only parametric models and it is in
general unclear how variational credible sets behave in high- and infinite-
dimensional settings. A related result we are aware of can be found in the
recent paper [27], where confidence statements for a variational posterior are
obtained in a special bandit-like regression setting similar to the results of
[23] for the true posterior. Although interesting, these results and the
underlying techniques do not transfer to the usual nonparametric regression
setting we consider in this paper.

In our analysis we focus on the popular and routinely used Gaussian Process
(GP) regression model. Exact computation of the posterior quickly becomes
infeasible in practice since the computational cost scales cubically in the
sample size. To overcome this problem, a sparse variational approximation
method was proposed in [26]. The variational class is parametrized by so
called inducing variables which are fitted to the posterior. This approach has
become increasingly popular in the machine learning community and has been
applied in various settings including deep Gaussian processes and solving
inverse problems. Recently theoretical guarantees were also derived for it. In
[6] the average Kullback-Leibler distance between the variational and true
posterior was studied, while in [17] minimax posterior contraction rates were
derived. However, non of these results provided guarantees for the frequentist
validity of the resulting uncertainty quantification, which is arguably one of
the main aims in the Bayesian analysis.

We focus in our theoretical analysis on a specific choice of inducing
variables which we call population spectral features, allowing (relatively)
tractable mathematical analysis. In our numerical analysis we also show that
with other choices of inducing variables similar behaviour is obtained. We
observe that in contrast to the simple parametric examples using mean-field
variational approximations, in the nonparametric GP regression framework the
variational posterior provides, from a frequentist perspective, reliable
uncertainty statements for appropriately tuned priors. In fact the good
coverage property does not depend on the number of inducing variables used in
the procedure. Besides the coverage of variational credible sets we also
derive lower bounds for the number of inducing variables one has to use to
achieve minimax contraction rates. This complements the contraction rate
guarantees given in [17], where the sharpness of the lower bound was
conjectured, but not verified. To achieve this we use a different proof
technique than in [17]. We apply a variational version of the standard kernel
ridge regression method [22]. This direct approach provides more control of
posterior properties. Finally, we apply our abstract results to two GP priors
with polynomially and exponentially decaying eigenvalues, respectively.

#### Contributions.

We summarize our contributions below:

  * –

We give an explicit formula for the contraction rate of the variational
posterior in terms of the prior and the true regression function. This gives a
condition on the prior and the minimal number of inducing variables in the
variational approximation needed to obtain minimax contraction rates.

  * –

If the number of inducing variables is too low, the contraction rate is sub-
optimal, regardless of the choice of prior.

  * –

Irrespective of the number of inducing variables, variational credible sets
cover the truths that are at least as smooth as the prior, whereas coverage
may be bad if the prior over-smoothes the true regression function.

#### Outline.

In the next section, we describe the Gaussian Process regression model studied
in this paper. We recall the details of the variational procedure with
inducing variables and derive a connection with kernel ridge regression, used
in the proofs. Lastly we introduce the specific choice of inducing variables
considered in the theoretical analysis. In Section 3, we develop a theory for
contraction rates. Section 4 consists of the theory on uncertainty
quantification. In Section 5 these results are applied to two specific priors,
with polynomially and exponentially decaying eigenvalues, respectively. We
conclude with a numerical analysis, including various inducing variable
methods in Section 6. The proofs are deferred to the Appendix. In Section A we
prove the more abstract results, while the proofs for the examples are given
in Section B.

#### Notation.

For sequences $a_{n},b_{n}$ of non-negative real numbers, we write
$a_{n}\lesssim b_{n}$ if there exists a constant $C>0$ such that $a_{n}\leq
Cb_{n}$ for all $n\in\mathbb{N}$. We write $a_{n}\asymp b_{n}$ if both
$a_{n}\lesssim b_{n}$ and $b_{n}\lesssim a_{n}$ hold. We indicate with an
apostrophe the transpose $A^{\prime}$ of a matrix $A$.

##  2 Variational approximations for Gaussian process regression

Throughout the paper we investigate the nonparametric regression model

$y_{i}=f(x_{i})+\varepsilon_{i},\qquad i=1,\ldots,n,$ (2.1)

with i.i.d. design points $x_{i}$ with respect to some common probability
measure $\mu$ on some $\mathcal{X}\subseteq\mathbb{R}^{d}$, and i.i.d. mean-
zero Gaussian measurement errors
$\varepsilon_{i}\sim\mathcal{N}(0,\sigma^{2})$, for some $\sigma^{2}>0$. We
view the unknown parameter $f$ as an element of the function space
$L^{2}(\mathcal{X},\mu)$. We endow $f$ with a centered Gaussian Process (GP)
prior $\Pi$ determined by its covariance kernel
$k:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$.

The GP prior is conjugate for the regression model (2.1), which means that the
posterior is also a GP. However, the computational and memory costs of
obtaining the posterior are $O(n^{3})$ and $O(n^{2})$, respectively, which is
prohibitive for large data sets. Therefore, in practice various approximation
methods are applied for inference, see [18] for a detailed discussion. One of
the most commonly used approach is the sparse variational approximation using
inducing variables, proposed by [26]. In the next subsection we give a brief
summary of inducing variable variational Bayes methods in general and we will
focus on a specific, analytically convenient version of the method, using
population spectral features.

###  2.1 Variational Bayes with inducing variables

In the variational framework, the posterior distribution is approximated by
projecting it onto an appropriately selected class $\mathcal{Q}$ of
probability measures on $L^{2}(\mathcal{X},\mu)$ with respect to the Kullback-
Leibler divergence. Letting $\Pi(\,\cdot\mid\boldsymbol{x},\boldsymbol{y})$
denote the true posterior, the variational posterior
$\Psi(\,\cdot\mid\boldsymbol{x},\boldsymbol{y})$ is defined as

$\Psi(\,\cdot\mid\boldsymbol{x},\boldsymbol{y})=\mathop{\rm
arg\,min}_{\Psi\in\mathcal{Q}}D_{\mathrm{KL}}(\Psi\|\Pi(\,\cdot\mid\boldsymbol{x},\boldsymbol{y})),$
(2.2)

where $D_{\mathrm{KL}}$ denotes Kullback-Leibler divergence.

In the inducing variables framework, the variational class $\mathcal{Q}$ is
constructed using a collection $\boldsymbol{u}=(u_{1},\ldots,u_{m})$ of
(known, specified) bounded linear functionals evaluated at $f$. The linearity
guarantees that the distribution of $\boldsymbol{u}$ under $\Pi$ is
$m$-dimensional multivariate Gaussian. Furthermore, the prior conditional on
$\boldsymbol{u}$ is another GP law. To obtain a low dimensional optimization
problem and preserve aspects of the prior distribution, [26] proposed to fit a
variational distribution of $\boldsymbol{u}$ to its posterior, while keeping
the conditional prior distribution for $f|\boldsymbol{u}$. More concretely,
the sparse inducing variable variational class consists of the distributions

$\displaystyle\Psi=\int\Pi(\,\cdot\mid\boldsymbol{u})\,d\Psi_{\boldsymbol{u}}(\boldsymbol{u}),$
(2.3)

where $\Psi_{\boldsymbol{u}}$ is any non-degenerate $m$-dimensional Gaussian
distribution, indexing the variational class $\mathcal{Q}$. This way the
posterior information is compressed into the fitted $m$-dimensional
distribution of $\boldsymbol{u}$. Nevertheless, $\Psi$ is still a
nonparametric distribution on $L^{2}(\mathcal{X},\mu)$ which is equivalent to
the prior.

The variational posterior is computed by finding the distribution
$\Psi_{\boldsymbol{u}}$ that minimizes the Kullback-Leibler divergence between
$\Psi$ given in (2.3) and the true posterior. As is shown in [26], a unique
solution exists and can be computed analytically. The corresponding
variational posterior is also GP with respective mean and covariance function

$\displaystyle\hat{f}_{m}^{\textrm{v}}(x)$
$\displaystyle=K_{x\boldsymbol{u}}(K_{\boldsymbol{uu}}+\sigma^{-2}K_{\boldsymbol{uf}}K_{\boldsymbol{fu}})^{-1}K_{\boldsymbol{uf}}\boldsymbol{y},$
(2.4) $\displaystyle\hat{k}_{m}^{\textrm{v}}(x,y)$
$\displaystyle=k(x,y)-K_{x\boldsymbol{u}}K_{\boldsymbol{uu}}^{-1}K_{\boldsymbol{u}y}+K_{x\boldsymbol{u}}(K_{\boldsymbol{uu}}+\sigma^{-2}K_{\boldsymbol{uf}}K_{\boldsymbol{fu}})^{-1}K_{\boldsymbol{u}y}.$
(2.5)

Here $\boldsymbol{y}=(y_{1},\ldots y_{n})$ is the vector of response
variables, $K_{\boldsymbol{uu}}$ is the covariance matrix of $\boldsymbol{u}$
under $\Pi$ (with entries $\operatorname{cov}_{\Pi}(u_{i},u_{j})$), and
similarly $K_{\boldsymbol{fu}}=K_{\boldsymbol{uf}}^{\prime}$ is the $n\times
m$ matrix with entries $\operatorname{cov}_{\Pi}(f(x_{i}),u_{j})=\Pi
f(x_{i})u_{j}$, and
$K_{\boldsymbol{u}x}=K_{x\boldsymbol{u}}^{\prime}=\operatorname{cov}_{\Pi}(\boldsymbol{u},f(x))$.
We note that in the special case $\boldsymbol{u}=\boldsymbol{f}$, the true
posterior is recovered. The above formulas were derived in [26] for the
inducing points method ($u_{j}=f(z_{j})$ for $z_{j}\in\mathcal{X}$,
$j=1,\ldots,m$), but the same computations hold for any choice of the inducing
variables. For completeness we provide the details in Section A.7.

The theoretical properties of this approach have been investigated in [6, 17]
for various choices of inducing variables. The first paper deals with the
accuracy of the variational approximation of the original posterior with
respect to the Kullback-Leibler divergence. In the second paper, upper bounds
were derived for the posterior contraction rate. We extend the latter result
by using a different, kernel ridge regression technique allowing sharper
control of the approximation. This analysis allows us to derive lower bounds
for the contraction rate and to investigate the frequentist coverage
properties of the credible sets resulting from the variational approximation.

###  2.2 Kernel ridge regression

The posterior mean, which is the maximum a posteriori in case of Gaussian
processes, can equivalently be formulated as a kernel ridge regression (KRR)
estimator. Let $\mathbb{H}$ be the reproducing kernel Hilbert space (RKHS)
associated with the GP kernel $k$. Then (see e.g. [14]) the true posterior
mean $\hat{f}_{n}$ solves

$\hat{f}_{n}=\mathop{\rm
arg\,min}_{f\in\mathbb{H}}\sum_{i=1}^{n}(y_{i}-f(x_{i}))^{2}+\sigma^{2}\|f\|_{\mathbb{H}}^{2}.$

It is not difficult to see that the variational posterior mean (2.4) can also
be viewed as a KRR estimator, for an appropriate choice of the RKHS. Since the
inducing variables $u_{j}$ are linear functionals of $f$, the functions

$h_{j}:\mathcal{X}\to\mathbb{R},~{}x\mapsto\operatorname{cov}_{\Pi}(f(x),u_{j})=\Pi
f(x)u_{j}$

are elements of $\mathbb{H}$ (see [30]). Let $\mathbb{H}_{m}$ denote the
linear subspace of $\mathbb{H}$ spanned by the functions $h_{1},\ldots,h_{m}$.
The variational posterior mean (2.4) is an element of $\mathbb{H}_{m}$. The
following lemma states that the inducing variable variational posterior mean
minimizes the same objective function as the mean of the true posterior, but
over the subclass $\mathbb{H}_{m}\subset\mathbb{H}$.

######  Lemma 1.

The variational posterior mean $\hat{f}_{m}^{\textrm{v}}$ given in (2.4)
satisfies

$\hat{f}_{m}^{\textrm{v}}=\mathop{\rm
arg\,min}_{f\in\mathbb{H}_{m}}\sum_{i=1}^{n}(y_{i}-f(x_{i}))^{2}+\sigma^{2}\|f\|_{\mathbb{H}}^{2}.$
(2.6)

For a similar result in context of the more specific Nyström Approximation
method we refer to [34]. Although the above lemma holds for arbitrary choices
of the inducing variables, in the upcoming sections we focus specifically on
the population spectral features approach. We believe that what we infer from
our results holds more generally, as illustrated in our numerical analysis.

###  2.3 Population spectral features

In our theoretical analysis we focus on a choice of inducing variables that
gives the variational posterior the interpretation of a spectral approximation
to the true posterior. We assume that the prior covariance kernel $k$ is
continuous and $k\in L^{2}(\mathcal{X}\times\mathcal{X},\mu\times\mu)$, so
that it has a Mercer decomposition

$k(x,y)=\sum_{j=1}^{\infty}\lambda_{j}\varphi_{j}(x)\varphi_{j}(y),$ (2.7)

where $\lambda_{j}$ is a decreasing, summable sequence of nonnegative numbers
and $(\varphi_{j})_{j\in\mathbb{N}}$ is an orthonormal basis of
$L^{2}(\mathcal{X},\mu)$. Under the current assumptions, the series (2.7)
converges not only in the $L^{2}$ sense but also uniformly on compact subsets
in the support of $\mu\times\mu$ (see [24], Corollary 3.5). We also associate
with $k$ an operator $T:L^{2}(\mathcal{X},\mu)\to L^{2}(\mathcal{X},\mu)$
given by

$Tg(x)=\int_{\mathcal{X}}k(x,y)g(y)\,d\mu(y),$ (2.8)

which is called the covariance operator. The identity (2.7) is equivalent to
the decomposition
$T=\sum_{j=1}^{\infty}\lambda_{j}\langle\cdot,\varphi_{j}\rangle\varphi_{j}$.
We assume that the set of functions $(\varphi_{j})_{j\in\mathbb{N}}$ is
uniformly bounded:

######  Assumption 2.

The functions $\varphi_{j}$ satisfy

$\sup\\{|\varphi_{j}(x)|:j\in\mathbb{N},x\in\mathcal{X}\\}=:C_{\varphi}<\infty.$
(2.9)

In the theoretical part of this paper, we consider the inducing variables

$u_{j}=\langle f,\varphi_{j}\rangle=\int f(x)\varphi_{j}(x)\,d\mu(x),\qquad
j=1,\ldots,m,$

which we refer to as population spectral features. We note that this approach
requires the explicit knowledge of the basis functions $\varphi_{j}$. Let us
write $\varphi_{1:m}(x)=(\varphi_{1}(x),\ldots,\varphi_{m}(x))$, and take
$\Phi$ to be the $n\times m$ matrix whose $i$-th row is
$\varphi_{1:m}(x_{i})^{\prime}$. It follows from Fubini’s theorem that $\Pi
u_{i}u_{j}=\lambda_{j}\delta_{ij}$ and $\Pi
f(x)u_{j}=\lambda_{j}\varphi_{j}(x)$, so in this case

$\displaystyle K_{\boldsymbol{uu}}$
$\displaystyle=\operatorname{diag}(\lambda_{1},\ldots,\lambda_{m})=:\Lambda,$
$\displaystyle K_{x\boldsymbol{u}}$
$\displaystyle=\varphi_{1:m}(x)^{\prime}\Lambda,$ (2.10) $\displaystyle
K_{\boldsymbol{fu}}$ $\displaystyle=\Phi\Lambda.$

In view of (2.4) and (2.5), the variational posterior
$\Psi(\,\cdot\mid\boldsymbol{x},\boldsymbol{y})$ is the law of a Gaussian
process with mean and covariance function

$\displaystyle\hat{f}_{m}^{\textrm{v}}(x)$
$\displaystyle=\varphi_{1:m}(x)^{\prime}(\Lambda^{-1}+\sigma^{-2}\Phi^{\prime}\Phi)^{-1}\Phi^{\prime}\boldsymbol{y}/\sigma^{2},$
(2.11) $\displaystyle\hat{k}_{m}^{\textrm{v}}(x,y)$
$\displaystyle=k(x,y)-\varphi_{1:m}(x)^{\prime}\Big{(}\Lambda-(\Lambda^{-1}+\sigma^{-2}\Phi^{\prime}\Phi)^{-1}\Big{)}\varphi_{1:m}(y),$
$\displaystyle=\varphi_{1:m}(x)^{\prime}(\Lambda^{-1}+\sigma^{-2}\Phi^{\prime}\Phi)^{-1}\varphi_{1:m}(y)+\sum_{j=m+1}^{\infty}\lambda_{j}\varphi_{j}(x)\varphi_{j}(y),$
(2.12)

where the last line follows from (2.7). We note that in this setting
$\mathbb{H}_{m}\subset\mathbb{H}$ is the linear span of the first $m$ basis
functions.

##  3 Contraction rates

First we investigate the asymptotic recovery property of the population
spectral features variational posterior. An upper bound for the contraction
rate was already derived in [17] for general inducing variables methods.
However, the implicit approach based on GP concentration functions do not
directly result in lower bounds for the contraction and therefore does not
imply a lower bound on the number of inducing variables one has to apply to
achieve minimax contraction rates. Therefore, in this article we take a more
direct approach using kernel ridge regression techniques to understand the
limitations of the variational approximation.

In this section we explicitly decompose the contraction rate to bias and
variance terms, which in turn illuminates how optimal contraction imposes a
condition on the minimal dimension $m$ of the approximation. It also shows why
the recovery accuracy does not improve any further by increasing the dimension
$m$ beyond this minimum. Moreover, a converse result to the contraction rate
statement is also given, which says that the variational posterior does not
contract at the optimal rate if $m$ is too low.

We assume that the data are generated according to some true $f_{0}\in
L^{2}(\mathcal{X},\mu)$. Let us denote by $\mathrm{P}_{0}$ the measure on
$\mathcal{X}^{n}\times\mathbb{R}^{n}$ under which
$(\boldsymbol{x},\boldsymbol{y})$ satisfies (2.1) with $f=f_{0}$. First, we
consider the variational posterior mean $\hat{f}_{m}^{\textrm{v}}$ as an
estimator of $f_{0}$. The next lemma decomposes the mean squared error of the
estimator $\hat{f}_{m}^{\textrm{v}}$ into squared bias and variance terms of
the estimator $\hat{f}_{m}^{\textrm{v}}$ under $\mathrm{P}_{0}$. The bias term
consists of two parts. The first part is accounting for the estimation error
in the subspace $\mathbb{H}_{m}$, while the second term equals the squared
norm of the orthogonal projection of $f_{0}$ onto the orthogonal complement
$\mathbb{H}_{m}^{\perp}$ of $\mathbb{H}_{m}$ in $L^{2}(\mathcal{X},\mu)$. The
proof of the lemma is deferred to Section A.6.

######  Lemma 3.

Define

$\nu_{j}=\frac{n\lambda_{j}}{\sigma^{2}+n\lambda_{j}}.$ (3.1)

Let $m=m_{n}$ be such that $m^{2}n^{-1}\log n\to 0$ as $n\to\infty$. Then for
any bounded $f_{0}\in L^{2}(\mathcal{X},\mu)$,

$\mathrm{P}_{0}\|\hat{f}_{m}^{\textrm{v}}-f_{0}\|^{2}\lesssim B_{n}+W_{n},$
(3.2)

where

$B_{n}=\sum_{j=1}^{m}(1-\nu_{j})^{2}\langle
f_{0},\varphi_{j}\rangle^{2}+\sum_{j>m}\langle
f_{0},\varphi_{j}\rangle^{2}\quad\text{and}\quad
W_{n}=\frac{1}{n}\sum_{j=1}^{m}\nu_{j}^{2}.$ (3.3)

######  Remark 4.

The mild technical assumption $m^{2}n^{-1}\log n$ allows for a clean
presentation. It can be weakened, but requires the introduction of a technical
term on the right-hand side of (3.2) as in [12].

The contraction rate of the variational posterior is determined by the terms
$B_{n},W_{n}$ introduced above, and additionally by the term $V_{n}$
introduced below, which characterises the spread of the variational posterior.
The proof of the next theorem is given in Section A.1.

######  Theorem 5.

Let $f_{0}\in L^{2}(\mathcal{X},\mu)$ be a bounded function and
$m=m_{n}\to\infty$ such that $m^{2}n^{-1}\log n\to 0$. Then

$\mathrm{P}_{0}\Psi(\|f-f_{0}\|\geq
M_{n}\epsilon_{n}\mid\boldsymbol{x},\boldsymbol{y})\to 0$ (3.4)

for arbitrary $M_{n}\to\infty$, where

$\epsilon_{n}^{2}=B_{n}+V_{n}+W_{n},$

with $B_{n}$ and $W_{n}$ as in (3.3), and

$V_{n}=\frac{1}{n}\sum_{j=1}^{m}\nu_{j}+\sum_{j>m}\lambda_{j}.$ (3.5)

In the remark below we investigate the bias and variance terms in more detail.
We present an alternative formulation which is more convenient to apply in our
examples and also sheds light on how the rate depends on the dimension $m$ of
the variational approximation.

######  Remark 6.

By considering separately the cases that $n\lambda_{j}\geq 1$ and
$n\lambda_{j}<1$, it follows that

$\nu_{j}\asymp 1\wedge n\lambda_{j},\qquad
1-\nu_{j}=\frac{\sigma^{2}}{n\lambda_{j}}\nu_{j}\asymp
1\wedge(n\lambda_{j})^{-1}.$ (3.6)

Let us introduce

$J_{n}=\max\\{j:n\lambda_{j}\geq 1\\},$ (3.7)

denoting the elbow point of the above quantities. Then the terms in the
contraction rate (3.4) can alternatively be written as

$\displaystyle B_{n}$ $\displaystyle\asymp\sum_{j=1}^{m\wedge
J_{n}}(n\lambda_{j})^{-2}\langle f_{0},\varphi_{j}\rangle^{2}+\sum_{j>m\wedge
J_{n}}\langle f_{0},\varphi_{j}\rangle^{2},$ (3.8) $\displaystyle V_{n}$
$\displaystyle\asymp\frac{m\wedge J_{n}}{n}+\sum_{j>m\wedge
J_{n}}\lambda_{j},\quad W_{n}\asymp\frac{m\wedge
J_{n}}{n}+n\sum_{j=J_{n}}^{m}\lambda_{j}^{2}.$ (3.9)

These identities follow immediately from (3.6). They show that the contraction
rate does not improve any further if $m$ is increased beyond $J_{n}$.

Next we derive lower bounds for the variational posterior contraction rate.
This in turn implies a lower bound on the number of inducing variables needed
to achieve minimax recovery of the truth. Let us assume that $f_{0}$ belongs
to the $\beta$-Sobolev space

$\mathcal{S}^{\beta}:=\\{f\in
L^{2}(\mathcal{X},\mu):\|f\|_{\beta}<\infty\\},\qquad\|f\|_{\beta}^{2}=\sum_{j=1}^{\infty}j^{2\beta/d}\langle
f,\varphi_{j}\rangle^{2},$ (3.10)

for some $\beta>0$. The minimax convergence rate for
$f_{0}\in\mathcal{S}^{\beta}$ is $n^{-\beta/(d+2\beta)}$.

The abstract results of Theorem 5 imply that for appropriately chosen
eigenvalues $\lambda_{j}$ the variational posterior contracts around the truth
with the minimax optimal rate, see Section 5 for two specific examples. In
both of these examples the number of inducing variables has to be at least of
order $n^{d/(d+2\beta)}$, which in fact corresponds to the $L_{2}$-entropy of
the $\beta$-Sobolev ball, also referred to as the “effective dimension” of the
model. The same lower bound on $m$ was observed for various choices of priors
and inducing variable methods in [17]. However, so far there was no
theoretical underpinning available for the sharpness of this threshold. The
following theorem aims to fill this gap by showing that optimal posterior
contraction can not be achieved when $m$ is below this threshold using the
population spectral features variational approach. The proof of the theorem is
given in Section A.2.

######  Theorem 7.

Let $m\asymp n^{r}$, where $0<r<\frac{d}{d+2\beta}$. Then there exists
$f_{0}\in\mathcal{S}^{\beta}$ and $0<p<2\beta/(d+2\beta)$ such that

$\mathrm{P}_{0}\|\hat{f}_{m}^{\textrm{v}}-f_{0}\|^{2}\gtrsim n^{-p},$ (3.11)

irrespective of the choice of prior. Moreover, possibly along a subsequence,
we have

$\mathrm{P}_{0}\Psi(\|f-f_{0}\|^{2}\leq
Mn^{-p}\mid\boldsymbol{x},\boldsymbol{y})\to 0$ (3.12)

for any $M>0$.

##  4 Uncertainty quantification

In this section we present our main results on the frequentist validity of
Bayesian uncertainty quantification resulting in from the variational
approximation. To this end, let us fix $\gamma\in(0,1)$ and consider the ball

$\displaystyle
C_{\gamma}:=\\{f:\|f-\hat{f}_{m}^{\textrm{v}}\|\leq\rho_{n}\\},$ (4.1)

where the radius $\rho_{n}$ is chosen such that $\Psi(f\in
C_{\gamma}\mid\boldsymbol{x},\boldsymbol{y})=1-\gamma$. This set is referred
to as the $(1-\gamma)$-credible ball of the variational posterior. In the next
theorem, we first study the asymptotic size of the radius $\rho_{n}$ under the
frequentist assumption that some $f_{0}\in L^{2}(\mathcal{X},\mu)$ generates
the data. In short, under $\mathrm{P}_{0}$, the asymptotic radius of a
credible set is of the order $V_{n}$, which was defined in (3.5). This
corroborates the remark made earlier that $V_{n}$ characterises the spread of
the variational posterior. The proof is given in Appendix A.3.

######  Theorem 8.

Suppose that $m=m_{n}$ is such that $m^{2}n^{-1}\log n\to 0$ as $n\to\infty$.
Then there exists a positive constant $C$ such that the credible ball
$C_{\gamma}$ defined in (4.1) has radius satisfying

$C^{-1}V_{n}\leq\rho_{n}^{2}\leq CV_{n},$

with $\mathrm{P}_{0}$-probability tending to $1$ for any $f_{0}\in
L^{2}(\mathcal{X},\mu)$.

We now consider the frequentist coverage of the credible set $C_{\gamma}$,
i.e. we are interested in the probability

$\displaystyle\mathrm{P}_{0}(f_{0}\in
C_{\gamma})=\mathrm{P}_{0}(\|\hat{f}_{m}^{\textrm{v}}-f_{0}\|\leq\rho_{n}).$
(4.2)

The next result provides guarantees but also limitations for achieving good
coverage. The proof is deferred to Section A.4.

######  Theorem 9.

For $B_{n},W_{n}$ and $V_{n}$ defined in (3.3) and (3.5), respectively,
consider the ratio

$R_{n}=\frac{B_{n}+W_{n}}{V_{n}}.$ (4.3)

Under the conditions of Theorem 8 and assuming $f_{0}\in
L^{2}(\mathcal{X},\mu)$ is bounded,

  1. 1.

if $R_{n}\to 0$, then
$\mathrm{P}_{0}(\|\hat{f}_{m}^{\textrm{v}}-f_{0}\|\leq\rho_{n})\to 1$;

  2. 2.

if $R_{n}\lesssim 1$, then
$\mathrm{P}_{0}(\|\hat{f}_{m}^{\textrm{v}}-f_{0}\|\leq M_{n}\rho_{n})\to 1$
for any sequence $M_{n}\to\infty$;

  3. 3.

if $R_{n}\to\infty$, then
$\mathrm{P}_{0}(\|\hat{f}_{m}^{\textrm{v}}-f_{0}\|\leq M\rho_{n})\to 0$ for
any $M>0$.

In view of (4.2) the above theorem presents the frequentist coverage
properties of the variational credible ball. It is determined by the relation
of the mean squared error
$\mathrm{P}_{0}\|\hat{f}_{m}^{\textrm{v}}-f_{0}\|^{2}$, studied in Lemma 3,
and the radius $\rho_{n}$ of the credible set, investigated in Theorem 8. The
first two statements are in line with the intuition that good coverage follows
from the credible set’s radius being larger than the loss. In case the radius
and loss are asymptotically comparable, good coverage can be achieved by
(slightly) blowing up the credible set with a growing factor $M_{n}$. The
third statement is the converse, i.e. if the loss exceeds the radius then
coverage will be bad.

We note that in statements 2 and 3, the ratio $R_{n}$ may be replaced by
$B_{n}/V_{n}$, since the variance of the posterior mean $W_{n}$ is always
bounded by the denominator $V_{n}$. The numerator $B_{n}$ represents the
(order of the) squared bias of the variational posterior mean and the
denominator $V_{n}$ corresponds to the variance of the variational posterior.
So Theorem 9 then characterizes coverage by a comparison of bias and variance;
the asymptotic coverage is good if variance dominates bias, and bad if the
bias strictly dominates the variance.

Below we demonstrate in examples that irrespective of the dimension $m$ of the
variational approximation, variational credible sets will cover truths that
are at least as smooth as the prior, and coverage may be bad if the prior
oversmoothes the truth. In this sense the variational posterior has the same
behaviour as the true posterior (see, for example, [15, 12]).

##  5 Examples

We apply our theoretical results from the previous sections to two different,
commonly used choices of the eigenvalues $\lambda_{j}$ of the kernel in (2.7).
First we investigate polynomially decaying eigenvalues. We note that Matérn
covariance kernels (including the Ornstein-Uhlenbeck process) and Riemann-
Liouville processes (including integrated Brownian motions) possess such
covariance structure. As a second example we consider exponentially decaying
eigenvalues. We note that the squared exponential Gaussian processes possess
such exponentially decaying eigenvalues. In our examples we choose the
eigenfunctions $(\varphi_{j})_{j\in\mathbb{N}}$ to meet our (uniform)
boundedness assumption in (2.9). This latter condition is not verified in
general, but only for specific choices of covariance kernels (e.g. the
Ornstein-Uhlenbeck process). Alternatively, it is possible to define a
suitable GP prior using (2.7) by specifying $(\lambda_{j})_{j\in\mathbb{N}}$
and taking any basis $(\varphi_{j})_{j\in\mathbb{N}}$ that satisfies the
boundedness assumption.

###  5.1 Polynomially decaying eigenvalues

In this subsection we consider eigenvalues of the form $\lambda_{j}\asymp
j^{-1-2\alpha/d}$, $j=1,\ldots,$ for some given $\alpha>0$, and eigenfunctions
$(\varphi_{j})_{j\in\mathbb{N}}$ satisfying (2.9). First, by applying Theorems
5 and 7 to the present setting, we derive contraction rates for the
corresponding population spectral features variational Bayes approximations.
The proof of the corollary is deferred to Section B.1.

######  Corollary 10.

Consider the kernel $k$ in (2.7) with $\lambda_{j}\asymp j^{-1-2\alpha/d}$ for
some $\alpha>d/2$ and $(\varphi_{j})_{j\in\mathbb{N}}$ satisfying (2.9). Then,
for $m\asymp n^{r}$ with $1/2>r\geq d/(d+2\alpha)$, the population spectral
features variational posterior contracts around bounded functions
$f_{0}\in\mathcal{S}^{\beta}$ at the rate
$\epsilon_{n}=n^{-(\beta\wedge\alpha)/(d+2\alpha)}$, which is minimax optimal
for $\alpha=\beta$. Furthermore, for $\alpha=\beta$ and $r<d/(d+2\alpha)$
there exists a $f_{0}\in\mathcal{S}^{\beta}$, such that $\epsilon_{n}\gtrsim
n^{-p}$, for some $p<\beta/(d+2\beta)$.

We note that for $1/2>r\geq d/(d+2\alpha)$ the variational approximation
inherits the asymptotic recovery property of the true posterior, i.e. the true
posterior contracts around $f_{0}\in\mathcal{S}^{\beta}$ with rate
$\epsilon_{n}=n^{-(\beta\wedge\alpha)/(d+2\alpha)}$, which is minimax optimal
for $\alpha=\beta$, see e.g. [31, 15]. The characterisation of the contraction
rate in Theorem 5 shows that this rate is achieved uniformly over Sobolev
balls $\\{f\in L^{2}(\mathcal{X},\mu):\|f\|_{\beta}\leq C\\}$. The upper bound
$m\lesssim n^{r}$ with $r<1/2$ is purely technical, allowing simpler
formulation of the general Theorem 5. In [17] the same contraction rate
results were derived using a different proof technique without this upper
bound on $m$. The lower bound $m\gtrsim n^{d/(d+2\alpha)}$, however, is sharp
and of importance. Fewer inducing variables can result in sub-optimal
posterior contraction around the truth, not matching the minimax rate.

Next we focus on the reliability of uncertainty quantification resulting in
from the variational Bayes approximation. We apply our general Theorem 9
providing guarantees but also limitations for the frequentist coverage for the
variational credible sets (4.1). Furthermore, in view of Theorem 8 we derive
an upper bound for the size of the credible ball. The proof of the next
corollary is given in Section B.2.

######  Corollary 11.

Let $k$ be the kernel in (2.7) with $\lambda_{j}\asymp j^{-1-2\alpha/d}$, for
some $\alpha>0$ and eigenfunctions $(\varphi_{j})_{j\in\mathbb{N}}$,
satisfying (2.9). Consider $m\asymp n^{r}$ population spectral inducing
variables for some $r\in(0,1/2)$.

  1. 1.

If $\alpha\leq\beta$, then

$\mathrm{P}_{0}(\|\hat{f}_{m}^{\textrm{v}}-f_{0}\|\leq M_{n}\rho_{n})\to 1$

for any bounded $f_{0}\in\mathcal{S}^{\beta}$ and any sequence
$M_{n}\to\infty$.

  2. 2.

If $\alpha<\beta$ and $r<d/(d+2\alpha)$, then

$\mathrm{P}_{0}(\|\hat{f}_{m}^{\textrm{v}}-f_{0}\|\leq\rho_{n})\to 1$

for any bounded $f_{0}\in\mathcal{S}^{\beta}$.

  3. 3.

If $\alpha>\beta>0$, then there exists $f_{0}\in\mathcal{S}^{\beta}$ such that

$\mathrm{P}_{0}(\|\hat{f}_{m}^{\textrm{v}}-f_{0}\|\leq M\rho_{n})\to 0$

for all $M>0$.

Finally, for $\beta=\alpha$, $m\geq n^{d/(d+2\alpha)}$ and $M>0$ there exists
$C>0$ such that for $\|f_{0}\|_{\beta}\leq M$, we have $P_{0}(\rho_{n}\leq
Cn^{-\beta/(d+2\beta)})\rightarrow 1$.

The first statement of the corollary shows that by considering a prior which
does not oversmooth the true function, the slightly inflated credible sets
provide reliable uncertainty quantification from the frequentist perspective.
In other words, by blowing up the credible balls (4.1) with a multiplicative
factor of $M_{n}\rightarrow\infty$, its frequentist coverage tends to one. We
note that inflating the credible set is in principle equivalent to (slightly)
undersmoothing the prior distribution. In statement 2, we in fact show that
such technical post-processing is not necessary by taking $\alpha<\beta$ and
considering less than $m\leq n^{r}$, with $r<d/(d+2\alpha)$, inducing
variables. We note that in view of Corollary 10, such choices of $\alpha$ and
$m$ result in an inflated contraction rate. In this case the variational
posterior mean is sub-optimally far from the true $f_{0}$ (considering squared
$L^{2}$-loss). At the same time the size of the credible ball is also of
higher magnitude $V_{n}\asymp\sum_{j>m}\lambda_{j}$ which is large enough to
dominate the expected loss of the variational posterior mean. When considering
smoother priors than the true function $f_{0}$, the frequentist coverage of
the variational credible set will tend to zero, independently from the number
of inducing variables applied. Finally, as a follow up of Corollary 10, the
size of the credible set achieves the minimax rate if the regularity of the
prior matches the regularity of $f_{0}$. All these results, in fact, are in
line with the frequentist coverage results derived for the true posterior
distribution [15].

###  5.2 Exponentially decaying eigenvalues

We also study the series prior with (rescaled) exponentially decaying
eigenvalues, motivated by the squared exponential kernel
$k(x,y)=\exp(-b\|x-y\|^{2})$. In this subsection we consider eigenvalues

$\lambda_{n,j}\asymp\exp(-\tau_{n}j^{1/d})\qquad j=1,\ldots,$ (5.1)

in (2.7) for some rescaling sequence $\tau_{n}\to 0$. For fixed $\tau_{n}$,
the decay of the eigenvalues is faster than in the previous example, giving
the associated prior support consisting only of functions that are
substantially smoother than the true function $f_{0}\in\mathcal{S}^{\beta}$.
To compensate for the rapidly diminishing eigenvalues a rescaling factor
$\tau_{n}$ is introduced, which shrinks the trajectories of the prior draws,
making them rougher. Let us consider the rescaling factor $\tau_{n}$ in the
form

$\tau_{n}:=n^{-1/(d+2\alpha)}\log n,$ (5.2)

for some constant $\alpha>0$. Taking $\alpha=\beta$, in fact, results in
minimax contraction rate and reliable uncertainty quantification for the true
posterior distribution, see for instance [29, 2, 13]. The following corollary
states that the contraction rate of the true posterior is inherited by its
variational approximation for sufficiently high $m$. The proof is deferred to
Section B.3.

######  Corollary 12.

Let $k$ be the kernel in (2.7) with eigenvalues satisfying (5.1) with
rescaling sequence (5.2) for some $\alpha>0$ and eigenfunctions $\varphi_{j}$,
$j=1,\ldots$ satisfying (2.9). Then, for $m=m_{n}\geq(\tau_{n}^{-1}\log
n)^{d}=n^{d/(d+2\alpha)}$ such that $m^{2}n^{-1}\log n\to 0$, the variational
posterior contracts around bounded $f_{0}\in\cal H^{\beta}$ at the rate
$\epsilon_{n}=n^{-(\beta\wedge\alpha)/(d+2\alpha)}$, which is minimax optimal
for $\alpha=\beta$. Furthermore, for $\alpha=\beta$ and $m=n^{r}$, with
$r<d/(d+2\alpha)$, there exists an $f_{0}\in\cal H^{\beta}$, for which the
posterior contracts at best at the sub-optimal rate $n^{-p}$, for some
$p<\beta/(d+2\beta)$.

This corollary, on one hand shows that the variational posterior achieves
minimax contraction rates for appropriately chosen rescaling factor and
sufficiently many inducing variables $m\gtrsim n^{d/(d+2\alpha)}$. (The
assumption $m^{2}n^{-1}\log n=o(1)$ is of technical nature, similarly to the
polynomial case.) On the other hand it also implies that we need at least this
many inducing variables, otherwise the variational posterior will provide sub-
optimal recovery for the truth.

Next we discuss the validity of the variational Bayes uncertainty
quantification. The result is slightly different in comparison with the
polynomially decaying eigenvalues. The proof of the next corollary is given in
Section B.4.

######  Corollary 13.

Let $k$ be the kernel in (2.7) with eigenvalues satisfying (5.1) with
rescaling sequence (5.2) for some $\alpha>0$ and eigenfunctions $\varphi_{j}$,
$j=1,\ldots,$ satisfying (2.9). Furthermore, let the approximation have
dimension $m=m_{n}\asymp n^{r}$ for some $r\in(0,1/2)$. Then

  1. 1.

if $r\geq d/(d+2\alpha)$ and $\beta\geq\alpha>0$, then for any bounded
$f_{0}\in\mathcal{S}^{\beta}$ and any $M_{n}\to\infty$,

$\mathrm{P}_{0}(\|\hat{f}_{m}^{\textrm{v}}-f_{0}\|\leq M_{n}\rho_{n})\to 1;$

  2. 2.

if $r<d/(d+2\alpha)$ then for any bounded $f_{0}\in L^{2}(\mathcal{X},\mu)$

$\mathrm{P}_{0}(\|\hat{f}_{m}^{\textrm{v}}-f_{0}\|\leq\rho_{n})\to 1;$

  3. 3.

if $\alpha>\beta>0$ and $m\geq n^{d/(d+2\alpha)}$, then there exists
$f_{0}\in\mathcal{S}^{\beta}$ such that for all $M>0$

$\mathrm{P}_{0}(\|\hat{f}_{m}^{\textrm{v}}-f_{0}\|\leq M\rho_{n})\to 0.$

Note that the only (real) difference between this result and Corollary 11 is
that in the second statement no assumption is made about the (Sobolev)
smoothness of $f_{0}$. If the dimension of the approximation is of order
$m_{n}=o(n^{d/(d+2\alpha)})$, then asymptotically the coverage of the
variational credible set will always be good. The rest of the results are the
same as in Corollary 11, except of the minor difference that in the third
statement the lower bound on $m$ holds exactly and not up to a multiplicative
constant due to the exponentially decaying form of the eigenvalues, hence the
conclusions are the same as well.

##  6 Numerical experiments

In this section we demonstrate on synthetic data sets that the variational
Bayes method proposed by [26] provides reliable uncertainty quantification
(from a frequentist perspective) for appropriately tuned Gaussian process
priors, independently from the number of inducing variables applied. In our
simulation study we go beyond the population spectral features inducing
variable method considered in our theoretical analysis above, and include
other, practically more advantageous and popular approaches as well. We aspire
to extrapolate our findings about the principles that govern the method with
the investigated particular choice of inducing variables to other choices of
inducing variables, for which the theoretical analysis is more complicated.
With this goal in mind, in the current section, we aim not only to illustrate
our findings but also compare various methods, to study empirically if and how
our theoretical results may carry over. Between each of three choices of
inducing variables, we shall point out what are common features and what is
different. Besides the population spectral features, we shall also study
inducing points and what we call ‘sample spectral features’.

In inducing points methods, the choice of inducing variables is
$u_{j}=f(z_{j})$ for a given set of inducing points
$z_{1},\ldots,z_{m}\in\mathcal{X}$. In the literature various approaches were
proposed to choose the inducing points. In our numerical analysis we consider
two specific choices. First we consider the equidistant choice of the inducing
points, i.e. we simply take the $z_{j}$ on an equispaced grid on $[-\pi,\pi]$.
This approach is designed to explore all parts of the signal equally. In case
the design points $(x_{i})_{i=1,\ldots,n}$ are not uniformly distributed the
grid can be modified accordingly to mimic the underlying distribution of the
covariates. We also consider the finite fixed-size determinantal point process
($m$-DPP), see [9].

The other choice of inducing variables is the sample analogue of the
population spectral features $u_{j}=\langle f,\varphi_{j}\rangle$. Instead of
diagonalizing the covariance operator $T$ in (2.8), we decompose the
covariance matrix of the prior at the design points
$K_{\boldsymbol{f\\!f}}=[k(x_{i},x_{j})]_{1\leq i,j\leq n}$. Since
$K_{\boldsymbol{f\\!f}}$ is positive definite, there exists an orthonormal set
of eigenvectors $\boldsymbol{v}_{1},\ldots,\boldsymbol{v}_{n}$ such that
$K_{\boldsymbol{f\\!f}}=\sum_{j=1}^{n}\mu_{j}\boldsymbol{v}_{j}\boldsymbol{v}_{j}^{\prime}$,
where $\mu_{1}\geq\mu_{2}\geq\cdots\geq\mu_{n}\geq 0$ are the eigenvalues of
$K_{\boldsymbol{f\\!f}}$. The sample spectral features are defined
$u_{j}:=\boldsymbol{v}_{j}^{\prime}\boldsymbol{f}=\sum_{i=1}^{n}v_{j,i}f(x_{i})$.
In [17] it was shown that the sample spectral features induce a variational
posterior with similar contraction rate results as the population spectral
features under similar threshold for the number of inducing variables.
Furthermore, the approximation accuracy of the above variational Bayes methods
to the true posterior distribution with respect to the expected Kullback-
Leibler divergence was studied in [6].

We consider the function space $L^{2}([-\pi,\pi],\mu)$, where $\mu$ is taken
as the uniform distribution on $[-\pi,\pi]$. As our underlying true function
(plotted in black in the figures) we take

$f_{0}=\sum_{\ell=1}^{\infty}\varphi_{3\ell}(3\ell)^{-1/2-\beta}/\log(3\ell)\in\mathcal{H}^{\beta}$

with $\beta=0.5$, where $\varphi_{j}$ denote the standard Fourier basis. We
use a synthetic data set with $n=2500$ realisations
$x_{i}\sim_{\textrm{i.i.d.}}\mu$, $i=1,\ldots,n$ and independent data points
$y_{i}\sim_{\text{ind}}\mathcal{N}(f_{0}(x_{i}),\sigma^{2})$ with
$\sigma=0.1$.

We consider both the polynomially and exponentially decaying eigenvalues
$\lambda_{j}$ from Section 5 and take the GP prior with covariance kernel
$k(x,y)=\sum_{j=1}^{\infty}\lambda_{j}\varphi_{j}(x)\varphi_{j}(y).$
Accordingly, in Figure 1 we plot the mean (dark gray) and 95% pointwise
credible sets of the associated true posterior when
$\lambda_{j}=j^{-1-2\beta}$ and $\lambda_{j}=\tau_{n}\exp(-\tau_{n}j/4)$, with
$\tau_{n}=n^{-1/(1+2\beta)}\log n$, respectively.

Figure 1: The true posterior distribution corresponding to the GP prior with
polynomially (left) and exponentially (right) decaying eigenvalues

First we investigate the variational Bayes approximations for the posterior
corresponding to the polynomial GP priors in Figure 2. We plot the population
spectral feature variational approximation in the first line (in blue), the
sample spectral feature approach in the second line (in red), the equidistant
inducing points method in the third line (in green) and the $m$-DPP inducing
points method in the fourth line (in yellow). The black curve stands for the
underlying true function $f_{0}$, the colored curve for the posterior mean,
while the shaded area represents the $95\%$ pointwise credible bands. We
consider two choices for the number of inducing variables. For $m=30$ (left
hand side) we are below, while for $m=60$ (right hand side) we are above the
theoretical threshold $m=n^{1/(1+2\beta)}=50$ obtained in our analysis. One
can observe that for all methods the size of the credible bands are overly
large for $m=30$, while for $m=60$ they closely resemble the true posterior
given on the left hand side of Figure 1 (with perhaps the exception of the
$m$-DPP method). Furthermore, for both choices of $m$ the credible bands
contain the true functional parameter $f_{0}$, illustrating the good
frequentist coverage properties of the variational methods, in line with our
theoretical findings. One can also notice that the approximations from the
population and sample spectral features are fairly similar. For the
equidistant inducing points method the variational posterior means do not seem
to differ significantly of those in the preceding plots, but credible regions
look fairly different. At the same time in the $m$-DPP method the posterior in
certain neighbourhoods can be quite different from the true posterior. In case
of both inducing points methods the intervals are narrow near those points on
the horizontal axis that are close to the inducing points, and widen as the
distance to an inducing point increases. Here, too, the credible regions seem
to be over-conservative, i.e. their width being about equal to that of the
credible set of the true posterior at inducing points and larger elsewhere.

Then in Figure 3 we plot the posterior means and credible sets resulting in
from exponentially decaying eigenvalues $(\lambda_{j})_{j\in\mathbb{N}}$. We
use the same experimental setup as for polynomial eigenvalues, i.e. we take
$m=30$ and $m=60$ inducing variables on the left and right hand side of the
figure, respectively, and consider the above discussed four variational
approximations. The plots are rather similar to what whas shown for the
polynomial eigenvalues in Figure 2, and hence the same conclusions hold for
this prior.

Figure 2: The variational Bayes approximations for GP with polynomially
decaying eigenvalues. for $m=30$ (left) and $m=60$ (right) inducing variables.
The choices of inducing variable methods from top to bottom: population
spectral features (blue), sample spectral features (red), inducing variables
with equidistant design (green) and with $m$-DPP (yellow). In each figure the
true function is plotted in black, the posterior mean is drawn with solid
colored curves, while the shaded area illustrate the $95\%$ point-wise
credible bands.

Figure 3: The variational Bayes approximations for GP with exponentially
decaying eigenvalues. for $m=30$ (left) and $m=60$ (right) inducing variables.
The choices of inducing variable methods from top to bottom: population
spectral features (blue), sample spectral features (red), inducing variables
with equidistant design (green) and with $m$-DPP (yellow). In each figure the
true function is plotted in black, the posterior mean is drawn with solid
colored curves, while the shaded area illustrate the $95\%$ point-wise
credible bands.

##  Appendix A Proof of the general theorems and lemmas

###  A.1 Proof of Theorem 5

In view of Lemma 3 and the inquality

$\mathrm{P}_{0}\Psi(\|f-f_{0}\|^{2}\geq(M_{n}\epsilon_{n})^{2}\mid\boldsymbol{x},\boldsymbol{y})\lesssim\frac{\mathrm{P}_{0}\Psi[\|f-\hat{f}_{m}^{\textrm{v}}\|^{2}|\boldsymbol{x},\boldsymbol{y}]+\mathrm{P}_{0}\|\hat{f}_{m}^{\textrm{v}}-f_{0}\|^{2}}{(M_{n}\epsilon_{n})^{2}}$

(by Markov’s inequality) it suffices to establish
$\mathrm{P}_{0}\Psi[\|f-\hat{f}_{m}^{\textrm{v}}\|^{2}|\boldsymbol{x},\boldsymbol{y}]\lesssim\epsilon_{n}^{2}$.
Using the formula for the posterior covariance (2.12), it follows that

$\displaystyle\mathrm{P}_{0}\Psi[\|f-\hat{f}_{m}^{\textrm{v}}\|^{2}\mid\boldsymbol{x},\boldsymbol{y}]$
$\displaystyle=\mathrm{P}_{0}\int\hat{k}_{m}^{\textrm{v}}(x,x)\,d\mu(x)$
$\displaystyle=\mathrm{P}_{0}\operatorname{tr}((\Lambda^{-1}+\sigma^{-2}\Phi^{\prime}\Phi)^{-1})+\sum_{j=m+1}^{\infty}\lambda_{j}.$
(A.1)

The proof is completed by showing

$\mathrm{P}_{0}\operatorname{tr}((\Lambda^{-1}+\sigma^{-2}\Phi^{\prime}\Phi)^{-1})\lesssim\frac{1}{n}\sum_{j=1}^{m}\nu_{j}.$
(A.2)

To do so, the expectation term is distributed over the event

$\Omega_{m}:=\Big{\\{}\sup_{1\leq j,k\leq
m}\Big{[}\frac{1}{n}\sum_{i=1}^{n}\varphi_{j}(x_{i})\varphi_{k}(x_{i})-\delta_{jk}\Big{]}^{2}>C\frac{\log
n}{n}\Big{\\}}$

and its complement. Note that

$\operatorname{tr}((\Lambda^{-1}+\sigma^{-2}\Phi^{\prime}\Phi)^{-1})\leq\operatorname{tr}\Lambda\leq\sum_{j=1}^{\infty}\lambda_{j}<\infty,$

hence by the union bound and Hoeffding’s inequality,

$\mathrm{P}_{0}1_{\Omega_{m}}\operatorname{tr}((\Lambda^{-1}+\sigma^{-2}\Phi^{\prime}\Phi)^{-1})\lesssim\mathrm{P}_{0}(\Omega_{m})\leq\sum_{j,k=1}^{m}2\exp\Big{(}\frac{-2n^{2}C\log
n}{4C_{\varphi}^{2}n^{2}}\Big{)}=2m^{2}n^{-C/2C_{\varphi}^{2}},$

where $C$ can be chosen arbitrarily large. This term is dominated by the
following upper bound that we establish below:

$\mathrm{P}_{0}1_{\Omega_{m}^{c}}\operatorname{tr}((\sigma^{2}\Lambda^{-1}+\Phi^{\prime}\Phi)^{-1})\lesssim\operatorname{tr}((\sigma^{2}\Lambda^{-1}+nI)^{-1})=\frac{1}{n}\sum_{j=1}^{m}\nu_{j}.$
(A.3)

Using the matrix identity $A^{-1}-B^{-1}=A^{-1}(B-A)B^{-1}$ it follows that

$\displaystyle|\operatorname{tr}((\sigma^{2}\Lambda^{-1}+\Phi^{\prime}\Phi)^{-1})-\operatorname{tr}((\sigma^{2}\Lambda^{-1}+nI)^{-1})|$
$\displaystyle=\operatorname{tr}\Big{(}(\sigma^{2}\Lambda^{-1}+\Phi^{\prime}\Phi)^{-1}(nI-\Phi^{\prime}\Phi)(\sigma^{2}\Lambda^{-1}+nI)^{-1}\Big{)}$
$\displaystyle\leq\|(\sigma^{2}\Lambda^{-1}+\Phi^{\prime}\Phi)^{-1}\|\|nI-\Phi^{\prime}\Phi\|\operatorname{tr}((\sigma^{2}\Lambda^{-1}+nI)^{-1}),$

where $\|A\|$ denotes the operator norm of the matrix $A\in\mathbb{R}^{m\times
m}$ with respect to the Euclidean norm on $\mathbb{R}^{m}$, which for
symmetric $A$ coincides with the largest absolute eigenvalue of $A$. We now
show that on $\Omega_{m}^{c}$, the product of the two norms in the preceding
display is $o(1)$, implying (A.3).

By Lemma 14 below, there exists $j$ between $1$ and $m$ such that the smallest
eigenvalue of $\sigma^{2}\Lambda^{-1}+\Phi^{\prime}\Phi$ (on $\Omega_{m}^{c}$)
is bounded from below by

$\displaystyle\sigma^{2}\lambda_{j}^{-1}+[\Phi^{\prime}\Phi]_{j,j}-\sum_{\begin{subarray}{c}1\leq
k\leq m\\\ k\neq j\end{subarray}}|[\Phi^{\prime}\Phi]_{j,k}|$
$\displaystyle=\sigma^{2}\lambda_{j}^{-1}+\Big{[}\sum_{i=1}^{n}\varphi_{j}(x_{i})^{2}-\sum_{\begin{subarray}{c}1\leq
k\leq m\\\ k\neq
j\end{subarray}}|\sum_{i=1}^{n}\varphi_{k}(x_{i})\varphi_{j}(x_{i})|\Big{]}$
$\displaystyle\geq n-m\sqrt{Cn\log n}\geq cn,$

where the latter inequality being true for some (deterministic) $c>0$ provided
$n$ is large enough, as follows from the assumption $m^{2}n^{-1}\log n\to 0$.
Consequently, the largest eigenvalue of the inverse
$(\sigma^{2}\Lambda^{-1}+\Phi^{\prime}\Phi)^{-1}$ satisfies

$\|(\sigma^{2}\Lambda^{-1}+\Phi^{\prime}\Phi)^{-1}\|\leq\frac{1}{cn}.$

Another application of Lemma 14 similarly shows that

$\|nI-\Phi^{\prime}\Phi\|\leq m\sqrt{Cn\log n}.$

Since $m\sqrt{n^{-1}\log n}\to 0$, the product of the above two norms vanishes
(deterministically) on $\Omega_{m}^{c}$, implying (A.3) and concluding the
proof.

######  Lemma 14 (Gershgorin circle theorem; see [10]).

Let $A\in\mathbb{R}^{n\times n}$ with entries $A_{ij}$. For any eigenvalue
$\lambda$ of $A$, there exists $j\in\\{1,\ldots,n\\}$ such that

$|\lambda-A_{jj}|\leq\sum_{i\neq j}|A_{ij}|.$

###  A.2 Proof of Theorem 7

Take

$f_{0}=\sum_{j=1}^{\infty}j^{-(1+p/r)/2}\varphi_{j},\qquad\text{for some
}\frac{2r\beta}{d}<p<\frac{2\beta}{d+2\beta}$

and note that $f_{0}\in\mathcal{S}^{\beta}$ and

$\sum_{j>m}\langle
f_{0},\varphi_{j}\rangle^{2}\asymp\int_{m}^{\infty}x^{-1-p/r}\,dx\asymp
m^{-p/r}\asymp n^{-p}\gg n^{-2\beta/(d+2\beta)}.$ (A.4)

The left-hand side in this display is a lower bound for
$\|\hat{f}_{m}^{\textrm{v}}-f_{0}\|^{2}$ since
$\langle\hat{f}_{m}^{\textrm{v}},\varphi_{j}\rangle=0$ for $j>m$, so (3.11)
follows.

For the contraction rate statement we fix
$p<p^{\prime}<\frac{2\beta}{d+2\beta}$ and make a case distinction, first
assuming

$\sum_{j>m}\lambda_{j}\lesssim n^{-p^{\prime}}.$ (A.5)

Note that given $M>0$, by (A.4) there exists $C_{0}>0$ such that for $n$ large
enough

$Mn^{-p^{\prime}}\leq(1-C_{0})\sum_{j>m}\langle f_{0},\varphi_{j}\rangle^{2}.$

Combining this with

$\|f-f_{0}\|^{2}=\sum_{j=1}^{\infty}\langle
f-f_{0},\varphi_{j}\rangle^{2}\geq\sum_{j>m}\langle
f_{0},\varphi_{j}\rangle^{2}-2\sum_{j>m}\langle
f_{0},\varphi_{j}\rangle\langle f,\varphi_{j}\rangle$

yields

$\displaystyle\Psi(\|f-f_{0}\|^{2}\leq
Mn^{-p^{\prime}}\mid\boldsymbol{x},\boldsymbol{y})$
$\displaystyle\leq\Psi\Big{(}\sum_{j>m}\langle
f_{0},\varphi_{j}\rangle^{2}-Mn^{-p^{\prime}}\leq 2\sum_{j>m}\langle
f_{0},\varphi_{j}\rangle\langle
f,\varphi_{j}\rangle\,\Big{|}\,\boldsymbol{x},\boldsymbol{y}\Big{)}$
$\displaystyle\leq\Psi\Big{(}\sum_{j>m}\langle
f_{0},\varphi_{j}\rangle^{2}\leq 2C_{0}\sum_{j>m}\langle
f_{0},\varphi_{j}\rangle\langle
f,\varphi_{j}\rangle\,\Big{|}\,\boldsymbol{x},\boldsymbol{y}\Big{)}$
$\displaystyle=\frac{1}{2}\Psi\Big{(}\sum_{j>m}\langle
f_{0},\varphi_{j}\rangle^{2}\leq 2C_{0}\big{|}\sum_{j>m}\langle
f_{0},\varphi_{j}\rangle\langle
f,\varphi_{j}\rangle\big{|}\,\Big{|}\,\boldsymbol{x},\boldsymbol{y}\Big{)},$

where in the last line we used that under the variational posterior the inner
products $\langle f,\varphi_{j}\rangle,j>m$ have a mean-zero Gaussian
distribution. Continuing the preceding display, and using the Markov, triangle
and Cauchy-Schwarz inequalities, it follows that

$\displaystyle\Psi(\|f-f_{0}\|^{2}\leq
Mn^{-p^{\prime}}\mid\boldsymbol{x},\boldsymbol{y})$
$\displaystyle\lesssim\frac{\sum_{j>m}\langle
f_{0},\varphi_{j}\rangle\Psi(|\langle
f,\varphi_{j}\rangle|\mid\boldsymbol{x},\boldsymbol{y})}{\sum_{j>m}\langle
f_{0},\varphi_{j}\rangle^{2}}$ $\displaystyle\lesssim n^{p}\sum_{j>m}\langle
f_{0},\varphi_{j}\rangle\sqrt{\lambda_{j}}$ $\displaystyle\leq
n^{p}\sqrt{\sum_{j>m}\langle
f_{0},\varphi_{j}\rangle^{2}\sum_{j>m}\lambda_{j}}$ $\displaystyle\lesssim
n^{(p-p^{\prime})/2}\to 0,$

where in the last line we used (A.4) and (A.5).

If the assumption (A.5) does not hold, then along a subsequence,

$n^{p^{\prime}}\sum_{j>m}\lambda_{j}\to\infty.$ (A.6)

In this case, we note that for $j>m$ the inner products $\langle
f,\varphi_{j}\rangle\stackrel{{\scriptstyle d}}{{=}}\sqrt{\lambda_{j}}Z_{j}$
with $Z_{j}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,1)$, hence by Markov’s
inequality,

$\displaystyle\Psi(\|f\|^{2}\leq
Mn^{-p^{\prime}}\mid\boldsymbol{x},\boldsymbol{y})$
$\displaystyle\leq\Psi\Big{(}\sum_{j>m}\langle f,\varphi_{j}\rangle^{2}\leq
Mn^{-p^{\prime}}\,\Big{|}\,\boldsymbol{x},\boldsymbol{y}\Big{)}$
$\displaystyle=\Psi\Big{(}\exp(-n^{p^{\prime}}{\textstyle\sum_{j>m}\langle
f,\varphi_{j}\rangle^{2}})\geq\exp(-M)\,\Big{|}\,\boldsymbol{x},\boldsymbol{y}\Big{)}$
$\displaystyle\lesssim\prod_{j>m}\mathbb{E}\exp\Big{(}-n^{p^{\prime}}\lambda_{j}Z_{j}^{2}\Big{)}=\prod_{j>m}(1+2n^{p^{\prime}}\lambda_{j})^{-1/2}$
$\displaystyle\leq\Big{(}1+\sum_{j>m}2n^{p^{\prime}}\lambda_{j}\Big{)}^{-1/2},$

which vanishes along a subsequence by the assumption (A.6). We conclude that
in this case the assertion of the theorem holds for $f_{0}=0$.

###  A.3 Proof of Theorem 8

In view of Markov’s inequality and the definition of $\rho_{n}$

$\rho_{n}^{2}=\rho_{n}^{2}\gamma^{-1}\Psi(\|f-\hat{f}_{m}^{\textrm{v}}\|>\rho_{n}\mid\boldsymbol{x},\boldsymbol{y})\leq\gamma^{-1}\Psi(\|f-\hat{f}_{m}^{\textrm{v}}\|^{2}\mid\boldsymbol{x},\boldsymbol{y}),$
(A.7)

Then the upper bound for $\rho_{n}$ is implied by the preceding display
together with the inequalities (A.1) and (A.2) in the proof of Theorem 5.

Next we establish the lower bound. First note, that in view of Fubini’s
theorem and the expansion (2.12) of the variational posterior covariance,

$\displaystyle\Psi(\langle
f-\hat{f}_{m}^{\textrm{v}},\varphi_{i}\rangle\langle
f-\hat{f}_{m}^{\textrm{v}},\varphi_{j}\rangle\mid\boldsymbol{x},\boldsymbol{y})$
$\displaystyle=\Psi\Big{(}\int(f-\hat{f}_{m}^{\textrm{v}})\varphi_{i}\,d\mu\int(f-\hat{f}_{m}^{\textrm{v}})\varphi_{j}\,d\mu\,\Big{|}\,\boldsymbol{x},\boldsymbol{y}\Big{)}$
$\displaystyle=\int\int\hat{k}_{m}^{\textrm{v}}(u,v)\varphi_{i}(u)\varphi_{j}(v)d\mu(u)d\mu(v)$
$\displaystyle=\begin{cases}[(\Lambda^{-1}+\sigma^{-2}\Phi^{\prime}\Phi)^{-1}]_{i,j}&\text{for
$i\vee j\leq m$,}\\\ \lambda_{j}\delta_{ij}&\text{for $i\vee
j>m$.}\end{cases}$

Therefore, under the variational posterior
$\Psi(\,\cdot\mid\boldsymbol{x},\boldsymbol{y})$ the inner products $\langle
f-\hat{f}_{m}^{\textrm{v}},\varphi_{j}\rangle$ are centered Gaussian random
variables, where the vector of the first $m$ variables has covariance matrix
$(\Lambda^{-1}+\sigma^{-2}\Phi^{\prime}\Phi)^{-1}$ and is independent of the
remaining $j>m$, which are mutually independent with variance $\lambda_{j}$.
Furthermore, note that if $Z$ has an $m$-dimensional normal distribution with
mean zero and covariance matrix $A$ with eigenvalues $a_{1},\ldots,a_{m}$,
then

$\displaystyle\mathbb{E}e^{-v\|Z\|^{2}}$
$\displaystyle=\mathbb{E}e^{-v\sum_{j=1}^{m}Z_{i}^{2}}$
$\displaystyle=\det(I+2vA)^{-1/2}$
$\displaystyle=\prod_{j=1}^{m}(1+2va_{j})^{-1/2}$
$\displaystyle\leq\Big{(}1+2v\sum_{j=1}^{m}a_{j}\Big{)}^{-1/2}$
$\displaystyle=\Big{(}1+2v\operatorname{tr}(A)\Big{)}^{-1/2}.$

Therefore, by using Markov’s inequality,

$\displaystyle 1-\gamma$
$\displaystyle=\Psi(\|f-\hat{f}_{m}^{\textrm{v}}\|\leq\rho_{n}\mid\boldsymbol{x},\boldsymbol{y})$
$\displaystyle=\Psi\big{(}\exp(-\rho_{n}^{-2}\|f-\hat{f}_{m}^{\textrm{v}}\|^{2})\geq\exp(-1)\mid\boldsymbol{x},\boldsymbol{y}\big{)}$
$\displaystyle\leq
e\Psi\Big{(}\exp\Big{(}-\rho_{n}^{-2}\sum_{j=1}^{\infty}\langle
f-\hat{f}_{m}^{\textrm{v}},\varphi_{j}\rangle^{2}\Big{)}\,\Big{|}\,\boldsymbol{x},\boldsymbol{y}\Big{)}$
$\displaystyle=e\Psi\Big{(}\exp\Big{(}-\rho_{n}^{-2}\sum_{j=1}^{m}\langle
f-\hat{f}_{m}^{\textrm{v}},\varphi_{j}\rangle^{2}\Big{)}\,\Big{|}\,\boldsymbol{x},\boldsymbol{y}\Big{)}\cdot\prod_{j>m}\Psi\Big{(}\exp(-\rho_{n}^{-2}\langle
f,\varphi_{j}\rangle^{2})\,\Big{|}\,\boldsymbol{x},\boldsymbol{y}\Big{)}$
$\displaystyle=e\det\Big{(}I+2\rho_{n}^{-2}(\Lambda^{-1}+\sigma^{-2}\Phi^{\prime}\Phi)^{-1}\Big{)}^{-1/2}\cdot\prod_{j>m}(1+2\rho_{n}^{-2}\lambda_{j})^{-1/2}$
$\displaystyle\leq
e\Big{(}1+2\rho_{n}^{-2}\operatorname{tr}\big{(}(\Lambda^{-1}+\sigma^{-2}\Phi^{\prime}\Phi)^{-1}\big{)}+2\rho_{n}^{-2}{\textstyle\sum_{j>m}\lambda_{j}}\Big{)}^{-1/2},$

which in turn implies that

$\rho_{n}^{2}\geq\frac{2}{e^{2}(1-\gamma)^{-2}-1}\Big{(}\operatorname{tr}\big{(}(\Lambda^{-1}+\sigma^{-2}\Phi^{\prime}\Phi)^{-1}\big{)}+{\textstyle\sum_{j>m}\lambda_{j}}\Big{)}.$

The lower bound in the statement of the theorem now follows by the following
inequality (12) from [32]:

$\operatorname{tr}\big{(}(\Lambda^{-1}+\sigma^{-2}\Phi^{\prime}\Phi)^{-1}\big{)}\geq\sum_{j=1}^{m}\frac{\sigma^{2}\lambda_{j}}{\sigma^{2}+\lambda_{j}\sum_{i=1}^{n}\varphi_{j}(x_{i})^{2}}\geq(C_{\varphi}^{-2}\wedge
1)\sigma^{2}\frac{1}{n}\sum_{j=1}^{m}\nu_{j}.$

###  A.4 Proof of Theorem 9

Consider first the statements 1 and 2. Note that, with $K$ as in Theorem 8, by
Markov’s inequality,

$\displaystyle\mathrm{P}_{0}(\|\hat{f}_{m}^{\textrm{v}}-f_{0}\|>M_{n}\rho_{n})$
$\displaystyle\leq\mathrm{P}_{0}(\rho_{n}^{2}<K^{-1}V_{n})+\mathrm{P}_{0}(\|\hat{f}_{m}^{\textrm{v}}-f_{0}\|^{2}\geq
M_{n}^{2}K^{-1}V_{n})$
$\displaystyle\lesssim\mathrm{P}_{0}(\rho_{n}^{2}<K^{-1}V_{n})+M_{n}^{-2}\frac{\mathrm{P}_{0}\|\hat{f}_{m}^{\textrm{v}}-f_{0}\|^{2}}{V_{n}}.$

This can be seen to vanish by applying Theorem 8 to the term on the left in
the upper bound, and by combining Lemma 3 with the assumptions on $R_{n}$ for
the term on the right. Statements 1 and 2 follow, the former by taking
$M_{n}=1$.

Since $\nu_{n}\leq 1$, the assumption $R_{n}\to\infty$ in statement 3 is
equivalent to $B_{n}/V_{n}\rightarrow\infty$. Then letting
$F:L^{2}(\mathcal{X},\mu)\to L^{2}(\mathcal{X},\mu)$ be the operator
$f\mapsto\sum_{j=1}^{m}\nu_{j}\langle f,\varphi_{j}\rangle\varphi_{j}$, and id
the identity operator, the assumption can be further re-written as

$V_{n}/\|(\mathrm{id}-F)f_{0}\|^{2}\to 0.$ (A.8)

Note that by the triangle inequality,

$\displaystyle\mathrm{P}_{0}\big{(}\|\hat{f}_{m}^{\textrm{v}}-f_{0}\|\leq\rho_{n}\big{)}$
$\displaystyle\leq\mathrm{P}_{0}\big{(}\|(\mathrm{id}-F)f_{0}\|\leq\|\hat{f}_{m}^{\textrm{v}}-Ff_{0}\|+M\rho_{n}\big{)}$
$\displaystyle\leq\mathrm{P}_{0}\big{(}\|(\mathrm{id}-F)f_{0}\|\leq
2\|\hat{f}_{m}^{\textrm{v}}-Ff_{0}\|\big{)}+\mathrm{P}_{0}\big{(}\|(\mathrm{id}-F)f_{0}\|\leq
2M\rho_{n}\big{)}.$

The second probability on the right hand side is seen to vanish as
$n\to\infty$ by combining (A.8) with Theorem 8. Regarding the first
probability, (A.10) in the proof of Lemma 3 below, gives

$\mathrm{P}_{0}\|\hat{f}_{m}^{\textrm{v}}-Ff_{0}\|^{2}\lesssim
n^{-1}\sum_{j=1}^{m}\nu_{j}^{2}$

so by Markov’s inequality and (A.8) it follows that

$\displaystyle\mathrm{P}_{0}\big{(}\|(\mathrm{id}-F)f_{0}\|$
$\displaystyle\leq 2\|\hat{f}_{m}^{\textrm{v}}-Ff_{0}\|\big{)}$
$\displaystyle\lesssim\frac{\mathrm{P}_{0}\|\hat{f}_{m}^{\textrm{v}}-Ff_{0}\|^{2}}{\|(\mathrm{id}-F)f_{0}\|^{2}}$
$\displaystyle\lesssim\frac{n^{-1}\sum_{j=1}^{m}\nu_{j}^{2}}{V_{n}}\frac{V_{n}}{\|(\mathrm{id}-F)f_{0}\|^{2}}\to
0.$

###  A.5 Proof of Lemma 1

Any $f\in\mathbb{H}_{m}$ is of the form
$f(x)=\sum_{j=1}^{m}a_{j}h_{j}(x)=K_{x\boldsymbol{u}}\boldsymbol{a}$ for some
$\boldsymbol{a}=(a_{1},\ldots,a_{m})\in\mathbb{R}^{m}$. Moreover, such an $f$
has squared RKHS norm

$\|f\|_{\mathbb{H}}^{2}=\sum_{j}\sum_{k}a_{j}a_{k}\langle
h_{j},h_{k}\rangle_{\mathbb{H}}=\sum_{j}\sum_{k}a_{j}a_{k}\Pi
u_{j}u_{k}=\boldsymbol{a}^{\prime}K_{\boldsymbol{uu}}\boldsymbol{a}.$

It follows that we need to minimise the objective function

$\boldsymbol{a}\mapsto(\boldsymbol{y}-K_{\boldsymbol{fu}}\boldsymbol{a})^{\prime}(\boldsymbol{y}-K_{\boldsymbol{fu}}\boldsymbol{a})+\sigma^{2}\boldsymbol{a}^{\prime}K_{\boldsymbol{uu}}\boldsymbol{a}.$

Setting the derivative equal to zero yields

$-K_{\boldsymbol{uf}}\boldsymbol{y}+(\sigma^{2}K_{\boldsymbol{uu}}+K_{\boldsymbol{uf}}K_{\boldsymbol{fu}})\boldsymbol{a}=0.$
(A.9)

Since the second derivative
$2(\sigma^{2}K_{\boldsymbol{uu}}+K_{\boldsymbol{uf}}K_{\boldsymbol{fu}})$ is
positive definite, the solution
$\boldsymbol{a}^{\star}=(K_{\boldsymbol{uu}}+\sigma^{-2}K_{\boldsymbol{uf}}K_{\boldsymbol{fu}})^{-1}K_{\boldsymbol{uf}}\boldsymbol{y}$
minimises the objective function. It follows that

$\arg\min_{f\in\mathbb{H}_{m}}\sum_{i=1}^{n}(y_{i}-f(x_{i}))^{2}+\sigma^{2}\|f\|_{\mathbb{H}}^{2}=K_{x\boldsymbol{u}}\boldsymbol{a}^{\star}=\hat{f}_{m}^{\textrm{v}}.$

###  A.6 Proof of Lemma 3

The proof follows similar lines of reasoning as [3] and [12] in context of GP
and distributed GP regression, using kernel ridge regression techniques. Here
these standard techniques are adapted to the variational approximation.

Letting $\mathrm{id}$ denote the identity operator and $F$ denote the operator

$F=\sum_{j=1}^{m}\nu_{j}\langle\,\cdot\,,\varphi_{j}\rangle\varphi_{j}$

on $L^{2}(\mathcal{X},\mu)$, we obtain an identity for the sum of the bias
terms

$\sum_{j=1}^{m}(1-\nu_{j})^{2}\langle
f_{0},\varphi_{j}\rangle^{2}+\sum_{j>m}\langle
f_{0},\varphi_{j}\rangle^{2}=\|(\mathrm{id}-F)f_{0}\|^{2}.$

Define $\Delta f_{0}=\hat{f}_{m}^{\textrm{v}}-Ff_{0}$. Since
$\|\hat{f}_{m}^{\textrm{v}}-f_{0}\|^{2}\leq 2\|\Delta
f_{0}\|^{2}+2\|(\mathrm{id}-F)f_{0}\|^{2}$, it suffices to show that

$\mathrm{P}_{0}\|\Delta
f_{0}\|^{2}\lesssim\frac{1}{n}\sum_{j=1}^{m}\nu_{j}^{2}.$ (A.10)

This is done by characterising the variational posterior mean
$\hat{f}_{m}^{\textrm{v}}$ as the root of a “score” function. Let us write
$\hat{f}_{m}^{\textrm{v}}=K_{x\boldsymbol{u}}\boldsymbol{a}^{\star}$ as in the
proof of Lemma 1. Combining (A.9) with the identities (2.10), it follows that
$\hat{f}_{m}^{\textrm{v}}=\varphi_{1:m}(\cdot)^{\prime}\Lambda\boldsymbol{a}^{\star}$
solves the equation

$\displaystyle\boldsymbol{0}$
$\displaystyle=\Lambda\Phi^{\prime}\boldsymbol{y}-(\sigma^{2}+\Lambda\Phi^{\prime}\Phi)\Lambda\boldsymbol{a}^{\star}$
$\displaystyle=\Lambda\Phi^{\prime}\left(\boldsymbol{y}-\begin{bmatrix}\hat{f}_{m}^{\textrm{v}}(x_{1})\\\
\vdots\\\
\hat{f}_{m}^{\textrm{v}}(x_{n})\end{bmatrix}\right)-\sigma^{2}\begin{bmatrix}\langle\hat{f}_{m}^{\textrm{v}},\varphi_{1}\rangle\\\
\vdots\\\ \langle\hat{f}_{m}^{\textrm{v}},\varphi_{m}\rangle\end{bmatrix}$
(A.11)

where we recall that
$\Lambda=\operatorname{diag}(\lambda_{1},\ldots,\lambda_{m})$ and
$\Phi_{i,j}=\varphi_{j}(x_{i})$. Multiplying the $j$-th entry of the vector in
(A.11) with $\varphi_{j}$ gives that $\hat{f}_{m}^{\textrm{v}}$ is the root of
the “score” operator $\hat{S}_{n}$

$\hat{S}_{n}(g)=\sum_{j=1}^{m}\Big{[}\frac{1}{n}\sum_{i=1}^{n}(y_{i}-g(x_{i}))\varphi_{j}(x_{i})\Big{]}\lambda_{j}\varphi_{j}-\frac{\sigma^{2}}{n}g,\qquad
g\in\mathbb{H}_{m},$

that is, $\hat{S}_{n}(\hat{f}_{m}^{\textrm{v}})=0$. Similarly, $Ff_{0}$ is the
root of the “population score”

$S_{n}(g)=\mathrm{P}_{0}\hat{S}_{n}(g)=\sum_{j=1}^{m}\langle
f_{0}-g,\varphi_{j}\rangle\lambda_{j}\varphi_{j}-\frac{\sigma^{2}}{n}g,\qquad
g\in\mathbb{H}_{m}.$

Let $f_{0,m}=\sum_{j=1}^{m}\langle f_{0},\varphi_{j}\rangle\varphi_{j}$ denote
the orthogonal projection of $f_{0}$ onto $\mathbb{H}_{m}$. Then

$S_{n}(\hat{f}_{m}^{\textrm{v}})=Tf_{0,m}-\sum_{j=1}^{m}\frac{n\lambda_{j}+\sigma^{2}}{n}\langle\hat{f}_{m}^{\textrm{v}},\varphi_{j}\rangle=T(f_{0,m}-F^{-1}\hat{f}_{m}^{\textrm{v}}),$

so $-\Delta
f_{0}=Ff_{0}-\hat{f}_{m}^{\textrm{v}}=FT^{-1}S_{n}(\hat{f}_{m}^{\textrm{v}})$
and

$\|\Delta f_{0}\|^{2}=\|FT^{-1}S_{n}(\hat{f}_{m}^{\textrm{v}})\|^{2}\leq
2\|FT^{-1}(\hat{S}_{n}(Ff_{0})+S_{n}(\hat{f}_{m}^{\textrm{v}}))\|^{2}+2\|FT^{-1}\hat{S}_{n}(Ff_{0})\|^{2}.$

We establish at the end of this proof that

$\mathrm{P}_{0}\|FT^{-1}(S_{n}(\hat{f}_{m}^{\textrm{v}})+\hat{S}_{n}(Ff_{0}))\|^{2}\lesssim
m^{4}n^{-D}+o(\mathrm{P}_{0}\|\Delta f_{0}\|^{2}),$ (A.12)

so the preceding two displays together yield

$\mathrm{P}_{0}\|\Delta
f_{0}\|^{2}\lesssim\mathrm{P}_{0}\|FT^{-1}\hat{S}_{n}(Ff_{0})\|^{2}+m^{4}n^{-D}+o(\mathrm{P}_{0}\|\Delta
f_{0}\|^{2}).$ (A.13)

Now we use that $y_{i}=f_{0}(x_{i})+\varepsilon_{i}$ under $\mathrm{P}_{0}$,
so

$\displaystyle\hat{S}_{n}(Ff_{0})$
$\displaystyle=\hat{S}_{n}(Ff_{0})-S_{n}(Ff_{0})$
$\displaystyle=\sum_{j=1}^{m}\Big{[}\frac{1}{n}\sum_{i=1}^{n}(y_{i}-Ff_{0}(x_{i}))\varphi_{j}(x_{i})-\langle(\mathrm{id}-F)f_{0},\varphi_{j}\rangle\Big{]}\lambda_{j}\varphi_{j}$
$\displaystyle=\sum_{j=1}^{m}\Big{[}\frac{1}{n}\sum_{i=1}^{n}(f_{0}-Ff_{0})(x_{i})\varphi_{j}(x_{i})-\langle(\mathrm{id}-F)f_{0},\varphi_{j}\rangle+\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}\varphi_{j}(x_{i})\Big{]}\lambda_{j}\varphi_{j}$

so since the $\varepsilon_{i}$ are independent of the $x_{i}$ and both are
i.i.d.,

$\displaystyle\mathrm{P}_{0}\|FT^{-1}\hat{S}_{n}(Ff_{0})\|^{2}$
$\displaystyle=\sum_{j=1}^{m}\nu_{j}^{2}\mathrm{P}_{0}\Big{[}\frac{1}{n}\sum_{i=1}^{n}[(\mathrm{id}-F)f_{0}](x_{i})\varphi_{j}(x_{i})-\langle(\mathrm{id}-F)f_{0},\varphi_{j}\rangle+\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}\varphi_{j}(x_{i})\Big{]}^{2}$
$\displaystyle=\sum_{j=1}^{m}\nu_{j}^{2}\frac{1}{n}\operatorname{var}_{0}\big{(}[(\mathrm{id}-F)f_{0}](x_{1})\varphi_{j}(x_{1})\big{)}+\frac{\sigma^{2}}{n}\sum_{j=1}^{m}\nu_{j}^{2}\lesssim\frac{1}{n}\sum_{j=1}^{m}\nu_{j}^{2},$

the bottom line following from the boundedness of $f_{0}$. This together with
(A.13) for $D$ large enough yields (A.10), which concludes the proof.

Proof of (A.12). Since
$S_{n}(Ff_{0})=\hat{S}_{n}(\hat{f}_{m}^{\textrm{v}})=0$,

$\displaystyle\hat{S}_{n}(Ff_{0})+S_{n}(\hat{f}_{m}^{\textrm{v}})$
$\displaystyle=\hat{S}_{n}(Ff_{0})-\hat{S}_{n}(\hat{f}_{m}^{\textrm{v}})+S_{n}(\hat{f}_{m}^{\textrm{v}})-S_{n}(Ff_{0})$
$\displaystyle=\sum_{j=1}^{m}\Big{[}\frac{1}{n}\sum_{i=1}^{n}(\Delta
f_{0})(x_{i})\varphi_{j}(x_{i})-\langle\Delta
f_{0},\varphi_{j}\rangle\Big{]}\lambda_{j}\varphi_{j}$

and using $\Delta f_{0}=\sum_{k=1}^{m}\langle\Delta
f_{0},\varphi_{k}\rangle\varphi_{k}$ and the Cauchy-Schwarz inequality, it
follows that

$\displaystyle\mathrm{P}_{0}\|FT^{-1}(\hat{S}_{n}(Ff_{0})+S_{n}(\hat{f}_{m}^{\textrm{v}}))\|^{2}$
$\displaystyle=\mathrm{P}_{0}\Big{\|}\sum_{j=1}^{m}\Big{[}\sum_{k=1}^{m}\langle\Delta
f_{0},\varphi_{k}\rangle\Big{[}\frac{1}{n}\sum_{i=1}^{n}\varphi_{k}(x_{i})\varphi_{j}(x_{i})-\delta_{jk}\Big{]}\Big{]}\nu_{j}\varphi_{j}\Big{\|}$
$\displaystyle=\sum_{j=1}^{m}\nu_{j}^{2}\mathrm{P}_{0}\Big{[}\sum_{k=1}^{m}\langle\Delta
f_{0},\varphi_{k}\rangle\Big{[}\frac{1}{n}\sum_{i=1}^{n}\varphi_{k}(x_{i})\varphi_{j}(x_{i})-\delta_{jk}\Big{]}\Big{]}^{2}$
$\displaystyle\leq\mathrm{P}_{0}\|\Delta
f_{0}\|^{2}\sum_{j=1}^{m}\nu_{j}^{2}\sum_{k=1}^{m}\Big{[}\frac{1}{n}\sum_{i=1}^{n}\varphi_{k}(x_{i})\varphi_{j}(x_{i})-\delta_{jk}\Big{]}^{2}.$
(A.14)

Splitting over the events

$\Omega_{j,m}:=\Big{\\{}\sup_{1\leq k\leq
m}\Big{[}\frac{1}{n}\sum_{i=1}^{n}\varphi_{j}(x_{i})\varphi_{k}(x_{i})-\delta_{jk}\Big{]}^{2}>C\frac{\log
n}{n}\Big{\\}}$ (A.15)

and their complements, and then using Assumption 2, the term in (A.14) is at
most of the order

$\displaystyle\frac{m\log n}{n}\sum_{j=1}^{m}\nu_{j}^{2}\mathrm{P}_{0}\|\Delta
f_{0}\|^{2}+C_{\varphi}^{2}m^{2}\max_{1\leq j\leq
m}\mathrm{P}_{0}1_{\Omega_{j,m}}\|\Delta f_{0}\|^{2}.$

Now (A.12) follows from

$\frac{m\log n}{n}\sum_{j=1}^{m}\nu_{j}^{2}\leq\frac{m^{2}\log n}{n}\to 0$

(by assumption) and, as we prove now

$\max_{1\leq j\leq m}\mathrm{P}_{0}1_{\Omega_{j,m}}\|\Delta
f_{0}\|^{2}\lesssim m^{2}n^{-D}.$ (A.16)

To this end, recall that from (2.11),
$\hat{f}_{m}^{\textrm{v}}(x)=\varphi_{1:m}(x)^{\prime}(\sigma^{2}\Lambda^{-1}+\Phi^{\prime}\Phi)^{-1}\Phi^{\prime}\boldsymbol{y}$,
so

$\displaystyle\|\hat{f}_{m}^{\textrm{v}}\|^{2}=\sum_{k=1}^{m}\langle\hat{f}_{m}^{\textrm{v}},\varphi_{k}\rangle^{2}$
$\displaystyle=\operatorname{tr}\Big{[}(\sigma^{2}\Lambda^{-1}+\Phi^{\prime}\Phi)^{-1}\Phi^{\prime}\boldsymbol{y}\boldsymbol{y}^{\prime}\Phi(\sigma^{2}\Lambda^{-1}+\Phi^{\prime}\Phi)^{-1}\Big{]}$
$\displaystyle\leq\lambda_{1}^{2}\sigma^{-4}\boldsymbol{y}^{\prime}\Phi\Phi^{\prime}\boldsymbol{y}\lesssim\sum_{j=1}^{m}\Big{[}\sum_{i=1}^{n}y_{i}\varphi_{j}(x_{i})\Big{]}^{2}\lesssim
mn\sum_{i=1}^{n}y_{i}^{2},$

and hence, since $y_{i}=f_{0}(x_{i})+\varepsilon_{i}$ for some bounded $f_{0}$
and $\varepsilon_{i}$ independent of $x_{i}$,

$\mathrm{P}_{0}1_{\Omega_{j,m}}\|\hat{f}_{m}^{\textrm{v}}\|^{2}\lesssim
mn\mathrm{P}_{0}1_{\Omega_{j,m}}\Big{[}\sum_{i=1}^{n}f_{0}(x_{i})^{2}+n\sigma^{2}\Big{]}\lesssim
mn^{2}\mathrm{P}_{0}(\Omega_{j,m}).$

This together with $\|Ff_{0}\|\leq\|f_{0}\|<\infty$ yields

$\mathrm{P}_{0}1_{\Omega_{j,m}}\|\Delta f_{0}\|^{2}\leq
2\mathrm{P}_{0}(\Omega_{j,m})\|Ff_{0}\|^{2}+2\mathrm{P}_{0}1_{\Omega_{j,m}}\|\hat{f}_{m}^{\textrm{v}}\|^{2}\lesssim
mn^{2}\mathrm{P}_{0}(\Omega_{j,m}).$

The inequality (A.16) follows from the above by combining a union bound and
Hoeffding’s inequality:

$\mathrm{P}_{0}(\Omega_{j,m})\leq m\cdot 2\exp\Big{(}\frac{-2n^{2}C\log
n}{4C_{\varphi}^{2}n^{2}}\Big{)}=2mn^{-D-2},$

with $D=C/2C_{\varphi}^{2}-2$, where the constant $C$ can be chosen
arbitrarily large.

###  A.7 Variational posterior for general inducing variables

Recall that the posterior is approximated by the variational distribution

$\Psi=\int\Pi(\cdot|\boldsymbol{u})\,d\Psi_{\boldsymbol{u}}(\boldsymbol{u})$
(A.17)

on $L^{2}(\mathcal{X},\mu)$. The above display is equivalent to

$\frac{d\Psi}{d\Pi}(f)=\frac{d\Psi_{\boldsymbol{u}}}{d\Pi_{\boldsymbol{u}}}(\boldsymbol{u}(f)),$
(A.18)

where $\Pi_{\boldsymbol{u}}$ is the distribution of the vector of inducing
variables ${\boldsymbol{u}}$ under the prior $\Pi$ (this follows from
existence of the Radon-Nikodym derivative on the right, which in turn is
guaranteed by the assumption that $\Psi_{\boldsymbol{u}}$ is a non-degenerate
Gaussian). The identity (A.18) shows that the variational posterior $\Psi$ is
parameterised by $f$ through ${\boldsymbol{u}}$ only.

Then, in view of Bayes’ theorem and some additional elementary algebraic
manipulations

$\displaystyle
D_{\mathrm{KL}}(\Psi\|\Pi(\,\cdot\mid\boldsymbol{x},\boldsymbol{y}))$
$\displaystyle=\int\log\frac{d\Psi}{d\Pi(\,\cdot\mid\boldsymbol{x},\boldsymbol{y})}\,d\Psi$
$\displaystyle=\int\log\frac{d\Psi}{d\Pi}\,d\Psi+\int\log\frac{d\Pi}{d\Pi(\,\cdot\mid\boldsymbol{x},\boldsymbol{y})}\,d\Psi$
$\displaystyle=\int\log\frac{d\Psi_{\boldsymbol{u}}}{d\Pi_{\boldsymbol{u}}}\,d\Psi_{\boldsymbol{u}}-\iint\log
p_{f}(X)\,d\Psi(f|{\boldsymbol{u}})\,d\Psi_{\boldsymbol{u}}({\boldsymbol{u}})+\log
p(X).$

In the variational procedure, we aim to minimise this with respect to $\Psi$.
Using that the KL-divergence is non-negative, and
$d\Psi(f|{\boldsymbol{u}})=d\Pi(f|{\boldsymbol{u}})$ by construction, we
obtain

$\log p(X)\geq\iint\log
p_{f}(X)\,d\Pi(f|{\boldsymbol{u}})\,d\Psi_{\boldsymbol{u}}({\boldsymbol{u}})-\int\log\frac{d\Psi_{\boldsymbol{u}}}{d\Pi_{\boldsymbol{u}}}\,d\Psi_{\boldsymbol{u}}.$

The right-hand side is commonly referred to as evidence lower bound (ELBO).
Minimising the KL-divergence is equivalent to maximising the ELBO. Note that
by Jensen’s inequality,

$\displaystyle\mathrm{ELBO}$
$\displaystyle=\int\log\Big{[}\frac{d\Pi_{\boldsymbol{u}}}{d\Psi_{\boldsymbol{u}}}({\boldsymbol{u}})\,{\textstyle\exp\Big{(}\int\log
p_{f}(X)\,d\Pi(f|{\boldsymbol{u}})\Big{)}}\Big{]}\,d\Psi_{\boldsymbol{u}}({\boldsymbol{u}})$
$\displaystyle\leq\log\int{\textstyle\exp\Big{(}\int\log
p_{f}(X)\,d\Pi(f|{\boldsymbol{u}})\Big{)}}\,d\Pi_{\boldsymbol{u}}({\boldsymbol{u}}),$
(A.19)

and the maximum is attained by the distribution $\Psi_{\boldsymbol{u}}^{*}$
defined through

$\frac{d\Psi_{\boldsymbol{u}}^{*}}{d\Pi_{\boldsymbol{u}}}({\boldsymbol{v}})=\frac{\exp\Big{(}\int\log
p_{f}(X)\,d\Pi(f|{\boldsymbol{u}}={\boldsymbol{v}})\Big{)}}{\int\exp\Big{(}\int\log
p_{f}(X)\,d\Pi(f|{\boldsymbol{u}})\Big{)}\,d\Pi_{\boldsymbol{u}}}.$ (A.20)

In general there is no guarantee that the maximizer
$\Psi_{\boldsymbol{u}}^{*}$ of the ELBO is tractable, but in case of the
considered Gaussian Process regression model and Gaussian variational class it
obtains an explicit form. In view of (A.20),

$\displaystyle\frac{d\Psi_{\boldsymbol{u}}^{*}}{d\Pi_{\boldsymbol{u}}}({\boldsymbol{u}})$
$\displaystyle\propto\exp\Big{(}-\frac{1}{2\sigma^{2}}\int\sum_{i=1}^{n}(y_{i}-f(x_{i}))^{2}\,d\Pi(f|{\boldsymbol{u}})\Big{)}$
$\displaystyle=\exp\Big{(}-\frac{1}{2\sigma^{2}}\sum_{i=1}^{n}\Big{[}(y_{i}-\Pi(f(x_{i})|{\boldsymbol{u}}))^{2}+\operatorname{var}_{\Pi(\cdot|{\boldsymbol{u}})}(f(x_{i}))\Big{]}\Big{)}$
$\displaystyle\propto\exp\Big{(}-\frac{1}{2\sigma^{2}}\sum_{i=1}^{n}(y_{i}-\Pi(f(x_{i})|{\boldsymbol{u}}))^{2}\Big{)}$
$\displaystyle=\exp\Big{(}-\frac{1}{2\sigma^{2}}\sum_{i=1}^{n}(y_{i}-K_{x_{i}u}K_{{\boldsymbol{u}}{\boldsymbol{u}}}^{-1}u)^{2}\Big{)},$
(A.21)

where
$K_{x{\boldsymbol{u}}}=\operatorname{cov}_{\Pi}(f(x),{\boldsymbol{u}})=\Pi
f(x){\boldsymbol{u}}^{\text{T}}$ and
$K_{{\boldsymbol{u}}{\boldsymbol{u}}}=\operatorname{cov}_{\Pi}({\boldsymbol{u}},{\boldsymbol{u}})=\Pi{\boldsymbol{u}}{\boldsymbol{u}}^{\text{T}}$,
and we used that

$\operatorname{var}_{\Pi(\cdot|{\boldsymbol{u}})}(f(x))=k(x,x)-K_{x{\boldsymbol{u}}}K_{{\boldsymbol{u}}{\boldsymbol{u}}}^{-1}K_{{\boldsymbol{u}}x},$

which is constant as a function of ${\boldsymbol{u}}$.

Completing the square in (A.21), it follows that the optimal variational
distribution $\Psi_{\boldsymbol{u}}^{*}$ of ${\boldsymbol{u}}$ is Gaussian
with mean
$K_{{\boldsymbol{u}}{\boldsymbol{u}}}(\sigma^{2}K_{{\boldsymbol{u}}{\boldsymbol{u}}}+K_{{\boldsymbol{u}}\boldsymbol{f}}K_{\boldsymbol{f}{\boldsymbol{u}}})^{-1}K_{{\boldsymbol{u}}\boldsymbol{f}}\boldsymbol{y}$
and covariance matrix
$K_{{\boldsymbol{u}}{\boldsymbol{u}}}(K_{{\boldsymbol{u}}{\boldsymbol{u}}}+\sigma^{-2}K_{{\boldsymbol{u}}\boldsymbol{f}}K_{\boldsymbol{f}{\boldsymbol{u}}})^{-1}K_{{\boldsymbol{u}}{\boldsymbol{u}}}$.
By (A.17), the variational distribution $\Psi^{*}$ of $f$ is Gaussian, with
mean function

$\displaystyle\hat{f}_{m}^{\textrm{v}}(x)$
$\displaystyle=\int\Pi(f(x)|{\boldsymbol{u}})\,d\Psi_{\boldsymbol{u}}^{*}({\boldsymbol{u}})$
$\displaystyle=K_{x{\boldsymbol{u}}}K_{{\boldsymbol{u}}{\boldsymbol{u}}}^{-1}\int{\boldsymbol{u}}\,d\Psi_{\boldsymbol{u}}^{*}({\boldsymbol{u}})$
$\displaystyle=K_{x{\boldsymbol{u}}}(\sigma^{2}K_{{\boldsymbol{u}}{\boldsymbol{u}}}+K_{{\boldsymbol{u}}\boldsymbol{f}}K_{\boldsymbol{f}{\boldsymbol{u}}})^{-1}K_{{\boldsymbol{u}}\boldsymbol{f}}\boldsymbol{y}.$

Furthermore, the covariance function is

$\displaystyle(x,y)$
$\displaystyle\mapsto\int\Pi((f-\hat{f}_{m}^{\textrm{v}})(x)(f-\hat{f}_{m}^{\textrm{v}})(y)|{\boldsymbol{u}})\,d\Psi^{*}_{\boldsymbol{u}}({\boldsymbol{u}})$
$\displaystyle=\int\operatorname{cov}_{\Pi(\cdot|{\boldsymbol{u}})}(f(x),f(y))\,d\Psi^{*}_{\boldsymbol{u}}({\boldsymbol{u}})+\operatorname{cov}_{\Psi^{*}_{\boldsymbol{u}}}(\Pi(f(x)|{\boldsymbol{u}}),\Pi(f(y)|{\boldsymbol{u}}))$
$\displaystyle=k(x,y)-K_{x{\boldsymbol{u}}}K_{{\boldsymbol{u}}{\boldsymbol{u}}}^{-1}K_{{\boldsymbol{u}}y}+K_{x{\boldsymbol{u}}}K_{{\boldsymbol{u}}{\boldsymbol{u}}}^{-1}\operatorname{cov}_{\Psi_{\boldsymbol{u}}^{*}}({\boldsymbol{u}},{\boldsymbol{u}})K_{{\boldsymbol{u}}{\boldsymbol{u}}}^{-1}K_{{\boldsymbol{\boldsymbol{u}}}y}$
$\displaystyle=k(x,y)-K_{x{\boldsymbol{u}}}K_{{\boldsymbol{u}}{\boldsymbol{u}}}^{-1}K_{{\boldsymbol{u}}y}+K_{x{\boldsymbol{u}}}(K_{{\boldsymbol{u}}{\boldsymbol{u}}}+\sigma^{-2}K_{{\boldsymbol{u}}\boldsymbol{f}}K_{\boldsymbol{f}{\boldsymbol{u}}})^{-1}K_{{\boldsymbol{u}}y}.$

##  Appendix B Proof of the corollaries

###  B.1 Proof of Corollary 10

Take $d/(d+2\alpha)\leq r<1/2$. Note that $m^{2}n^{-1}\log n=o(1)$, hence the
conditions of Theorem 5 hold. Next note that

$J_{n}=\max\\{j:n\lambda_{j}\geq 1\\}\asymp n^{d/(d+2\alpha)}\lesssim m,$
(B.1)

hence in view of (3.9),

$\displaystyle V_{n}$ $\displaystyle\asymp J_{n}/n+\int_{m\wedge
J_{n}}^{\infty}x^{-1-2\alpha/d}\,dx\asymp J_{n}/n+J_{n}^{-2\alpha/d}\asymp
n^{-2\alpha/(d+2\alpha)}$ $\displaystyle W_{n}$ $\displaystyle\lesssim
J_{n}/n+\int_{J_{n}}^{\infty}x^{-2-4\alpha/d}\,dx\asymp
J_{n}/n+nJ_{n}^{-1-4\alpha/d}\asymp n^{-2\alpha/(d+2\alpha)}.$ (B.2)

To bound the bias term $B_{n}$, we consider (3.8) and deal with the two terms
on the right hand side separately. For the first term note

$\displaystyle\sum_{j=1}^{m\wedge J_{n}}(n\lambda_{j})^{-2}\langle
f_{0},\varphi_{j}\rangle^{2}$
$\displaystyle\lesssim\sum_{j=1}^{J_{n}}(nj^{-1-2\alpha/d})^{-2}\langle
f_{0},\varphi_{j}\rangle^{2}$
$\displaystyle=n^{-2}\sum_{j=1}^{J_{n}}j^{2+4\alpha/d-2\beta/d}j^{2\beta/d}\langle
f_{0},\varphi_{j}\rangle^{2}$ $\displaystyle\lesssim n^{-2}\max_{1\leq
j<J_{n}}j^{2+4\alpha/d-2\beta/d}\|f_{0}\|_{\beta}^{2}$
$\displaystyle\lesssim(n^{-2}\vee
n^{-2\beta/(d+2\alpha)})\|f_{0}\|_{\beta}^{2},$

while for the second term we get

$\sum_{j>m\wedge J_{n}}^{\infty}\langle
f_{0},\varphi_{j}\rangle^{2}<(J_{n}\wedge m)^{-2\beta/d}\sum_{j>m\wedge
J_{n}}^{\infty}j^{2\beta/d}\langle f_{0},\varphi_{j}\rangle^{2}\lesssim
n^{-2\beta/(d+2\alpha)}\|f_{0}\|_{\beta}^{2}.$ (B.3)

hence the contraction rate in (3.4) is
$\epsilon_{n}=n^{-(\beta\wedge\alpha)/(d+2\alpha)}$.

The sub-optimal contraction rate for insufficient amount of inducing variables
is a direct consequence of Theorem 7.

###  B.2 Proof of Corollary 11

Note that $m^{2}n^{-1}\log n=n^{2r-1}\log n\to 0$, hence we can apply Theorem
9 to the present setting.

By (3.9), similarly to (B.2), we get that

$V_{n}\asymp\frac{m\wedge
n^{d/(d+2\alpha)}}{n}+m^{-2\alpha/d}\quad\text{and}\quad
W_{n}\lesssim\frac{m\wedge n^{d/(d+2\alpha)}}{n}.$ (B.4)

Then similarly as in the proof of Corollary 10 we get for
$f_{0}\in\mathcal{S}^{\beta}$, that

$\displaystyle\sum_{j=1}^{m}(1-\nu_{j})^{2}\langle
f_{0},\varphi_{j}\rangle^{2}$ $\displaystyle\lesssim(n^{-2}\vee
n^{-2\beta/(d+2\alpha)})$ $\displaystyle\sum_{j=m+1}^{\infty}\langle
f_{0},\varphi_{j}\rangle^{2}$ $\displaystyle\lesssim
m^{-2\beta/d}\|f_{0}\|_{\beta}^{2}.$

Therefore,

$\displaystyle R_{n}$
$\displaystyle=\frac{B_{n}+W_{n}}{V_{n}}\lesssim\frac{n^{-1}(m\wedge
n^{d/(d+2\alpha)})+n^{-2\beta/(d+2\alpha)}+m^{-2\beta/d}}{n^{-1}(m\wedge
n^{d/(d+2\alpha)})+m^{-2\alpha/d}}.$ (B.5)

In the first case, (B.5) is further bounded from above by a multiple of
$m^{-2\beta/d}/m^{-2\alpha/d}=O(1)$. In the second case, (B.5) is $o(1)$.
Finally, in the third case, let us take $f_{0}\in\mathcal{S}^{\beta}$ given by
$f_{0}=\sum_{j=1}^{\infty}j^{-1/2-q}\varphi_{j}$, for some
$\beta/d<q<\alpha/d$. Then in view of (3.8),

$\sum_{j=1}^{m\wedge J_{n}}(1-\nu_{j})^{2}\langle
f_{0},\varphi_{j}\rangle^{2}+\sum_{j>m\wedge J_{n}}\langle
f_{0},\varphi_{j}\rangle^{2}\gtrsim\sum_{j>m\wedge J_{n}}\langle
f_{0},\varphi_{j}\rangle^{2}\asymp(m\wedge J_{n})^{-2q}$

hence by recalling (B.4),

$R_{n}\gtrsim\frac{m^{-2q}\vee
n^{-2qd/(d+2\alpha)}}{m^{-2\alpha/d}+n^{-2\alpha/(d+2\alpha)}}\to\infty.$

Therefore, in all three cases the statement follows from Theorem 9.

Finally, the upper bound for the radius follows from Theorem 8 and (B.2) with
$\alpha=\beta$.

###  B.3 Proof of Corollary 12

The proof goes similarly to Corollary 10. First assume that $m$ is lower
bounded by

$J_{n}:=(\tau_{n}^{-1}\log n)^{d}=n^{d/(d+2\alpha)}$ (B.6)

and note that $\lambda_{n,J_{n}}\asymp\exp(-\tau_{n}J_{n}^{1/d})=n^{-1}$.
Hence, similarly to (3.8)

$B_{n}\asymp\sum_{j=1}^{m\wedge J_{n}}(n\lambda_{n,j})^{-2}\langle
f_{0},\varphi_{j}\rangle^{2}+\sum_{j>m\wedge J_{n}}\langle
f_{0},\varphi_{j}\rangle^{2}.$

We deal with the two terms on the right hand side of the preceding display
separately.

For any $f_{0}\in\mathcal{S}^{\beta}$, we have

$\displaystyle\sum_{j>m\wedge J_{n}}\langle f_{0},\varphi_{j}\rangle^{2}\leq
J_{n}^{-2\beta/d}\sum_{j>m\wedge J_{n}}j^{2\beta/d}\langle
f_{0},\varphi_{j}\rangle^{2}\lesssim
n^{-2\beta/(d+2\alpha)}\|f_{0}\|_{\beta}^{2}$ (B.7)

and

$\displaystyle\sum_{j=1}^{m\wedge J_{n}}(n\lambda_{n,j})^{-2}\langle
f_{0},\varphi_{j}\rangle^{2}$ $\displaystyle\leq\max_{1\leq i\leq
J_{n}}(n\lambda_{n,i})^{-2}i^{-2\beta/d}\sum_{j=1}^{\infty}j^{2\beta/d}\langle
f_{0},\varphi_{j}\rangle^{2}$ $\displaystyle\lesssim
n^{-2}\Big{(}\exp(2\tau_{n})\vee\frac{\exp(2\tau_{n}J_{n}^{1/d})}{J_{n}^{2\beta/d}}\Big{)}\|f_{0}\|_{\beta}^{2}$
$\displaystyle\lesssim(n^{-2}\vee
n^{-2\beta/(d+2\alpha)})\|f_{0}\|_{\beta}^{2}$ (B.8)

where we used that the function $i\mapsto\exp(2\tau_{n}i)i^{-2\beta/d}$ is
convex in $i$, so the maximum occurs at one of the endpoints.

Next we deal with the variance terms $W_{n}$ and $V_{n}$, defined in (3.3) and
(3.5), respectively. Similarly to (3.9) and (B.2), for $J_{n}$ given in (B.6),

$\displaystyle V_{n}\asymp\frac{m\wedge J_{n}}{n}+\sum_{j>m\wedge
J_{n}}\lambda_{n,j}$ $\displaystyle\asymp\frac{(m\wedge
J_{n})}{n}+\int_{m\wedge J_{n}}e^{-\tau_{n}x^{1/d}}dx$ (B.9)
$\displaystyle\lesssim n^{-\frac{2\alpha}{d+2\alpha}}+\int_{(m\wedge
J_{n})^{1/d}}x^{d-1}e^{-\tau_{n}x}\,dx.$

By partial integration and induction we get that the rightmost integral
satisfies

$\displaystyle\int_{(m\wedge
J_{n})^{1/d}}x^{d-1}e^{-\tau_{n}x}\,dx\asymp(\tau_{n}^{-1}(m\wedge
J_{n})^{1-1/d}\vee\tau_{n}^{-d})\exp(-\tau_{n}(m\wedge J_{n})^{1/d}).$ (B.10)

Since $m\geq J_{n}=(\tau_{n}^{-1}\log n)^{d}$ it is further bounded by

$n^{d/(d+2\alpha)}(\log n)^{-1}\exp(-\log n)=n^{-2\alpha/(d+2\alpha)}/\log n.$

Furthermore, by similar computations

$\displaystyle W_{n}$ $\displaystyle\asymp\frac{m\wedge
J_{n}}{n}+\sum_{j>m\wedge J_{n}}\lambda_{n,j}$ $\displaystyle\lesssim
n^{-\frac{2\alpha}{d+2\alpha}}+n\int_{(m\wedge
J_{n})^{1/d}}x^{d-1}e^{-2\tau_{n}x}\,dx$ $\displaystyle\lesssim
n^{-\frac{2\alpha}{d+2\alpha}}+n(\tau_{n}^{-1}(m\wedge
J_{n})^{1-1/d}\vee\tau_{n}^{-d})\exp(-2\tau_{n}(m\wedge J_{n})^{1/d})$
$\displaystyle\lesssim n^{-\frac{2\alpha}{d+2\alpha}}.$ (B.11)

Our contraction rate result follows from Theorem 5 with
$\epsilon_{n}^{2}=n^{-2(\beta\wedge\alpha)/(d+2\alpha)}$.

The sub-optimal contraction rate for insufficient amount of inducing variables
is a direct consequence of Theorem 7.

###  B.4 Proof of Corollary 13

Note that by assumption we have $m^{2}n^{-1}\log n\to 0$, hence we can apply
Theorem 9. Therefore it is sufficient to investigate the asymptotic behaviour
of the fraction $R_{n}=(B_{n}+W_{n})/V_{n}$, with the terms defined in (3.3)
and (3.5). We also recall the definition of $J_{n}$ given in (B.6) for the
exponentially decaying eigenvalues.

In the second case, following from the bounds (B.7), (B.8) and (B.11), the
numerator of $R_{n}$ is $o(1)$. At the same time, in view of assertions (B.9),
(B.10) and $\tau_{n}m^{1/d}=n^{r/d-1/(d+2\alpha)}\log n\to 0$, the denominator
is bounded from below as

$\displaystyle
V_{n}\gtrsim\sum_{j>m}\lambda_{n,j}\gtrsim\tau_{n}^{-d}\exp(-\tau_{n}m^{1/d})\gtrsim\tau_{n}^{-d}\to\infty,$

and therefore $R_{n}=o(1)$.

Next, in the first case, we have $m\wedge J_{n}\asymp n^{r}\wedge J_{n}\asymp
J_{n}$, hence again in view of (B.7), (B.8), $B_{n}\lesssim n^{-2}\vee
n^{-2\beta/(d+2\alpha)}$. Moreover, following from (B.9) and (B.11), we get
that $V_{n}\gtrsim J_{n}/n=n^{-2\alpha/(d+2\alpha)}$ and $W_{n}\lesssim
n^{-2\alpha/(d+2\alpha)}$, respectively. Combining this upper bounds results
in that $R_{n}=O(1)$.

Finally, in the third case, as in the proof of Corollary 11, let $f_{0}$ be
the function $\sum_{j=1}^{\infty}j^{-1/2-q}\varphi_{j}\in\mathcal{S}^{\beta}$
for some $\beta/d<q<\alpha/d$. In view of (3.8)

$B_{n}\gtrsim\sum_{j>m\wedge J_{n}}\langle
f_{0},\varphi_{j}\rangle^{2}\asymp(m\wedge J_{n})^{-2q}\gtrsim
n^{-2qd/(d+2\alpha)}.$

Furthermore, (B.9) and (B.10) together with $m\geq J_{n}$ imply that,

$V_{n}\lesssim
J_{n}/n+\tau_{n}^{-1}J_{n}^{1-1/d}\exp(-\tau_{n}J_{n}^{1/d})\lesssim
n^{-2\alpha/(d+2\alpha)}.$

Therefore, $R_{n}\gtrsim
n^{-2qd/(d+2\alpha)}/n^{-2\alpha/(d+2\alpha)}\to\infty$. We conclude that, in
all three cases the statements now follow directly from Theorem 9.

Finally, the upper bound for the radius follows from Theorem 8 and the upper
bound for $V_{n}$ in the proof of Corollary B.3 with $\alpha=\beta$.

## References

  * [1] [author] Alquier, PierreP. and Ridgway, JamesJ. (2020). Concentration of tempered posteriors and of their variational approximations. The Annals of Statistics 48 1475–1497. 10.1214/19-AOS1855
  * [2] [author] Bhattacharya, A.A. and Pati, D.D. (2015). Adaptive Bayesian inference in the Gaussian sequence model using exponential-variance priors. Statistics & Probability Letters 103 100–104.
  * [3] [author] Bhattacharya, A.A., Pati, D.D. and Yang, Y.Y. (2017). Frequentist coverage and sup-norm convergence rate in Gaussian process regression. arXiv preprint arXiv:1708.04753.
  * [4] [author] Bishop, Christopher MC. M. and Nasrabadi, Nasser MN. M. (2006). Pattern recognition and machine learning. Springer.
  * [5] [author] Blei, David MD. M., Kucukelbir, AlpA. and McAuliffe, Jon DJ. D. (2017). Variational Inference: A Review for Statisticians. Journal of the American statistical Association 112 859–877.
  * [6] [author] Burt, David R.D. R., Rasmussen, Carl EdwardC. E. and van der Wilk, MarkM. (2020). Convergence of Sparse Variational Inference in Gaussian Processes Regression. Journal of Machine Learning Research 21 1–63.
  * [7] [author] Castillo, IsmaëlI. and Nickl, RichardR. (2014). On the Bernstein–von Mises phenomenon for nonparametric Bayes procedures. The Annals of Statistics 42 1941–1969.
  * [8] [author] Cox, Dennis DD. D. (1993). An analysis of Bayesian inference for nonparametric regression. The Annals of Statistics 21 903–923.
  * [9] [author] Gautier, GuillaumeG., Polito, GuillermoG., Bardenet, RémiR. and Valko, MichalM. (2019). DPPy: DPP Sampling with Python. Journal of Machine Learning Research 20 1–7.
  * [10] [author] Gershgorin, Semyon AranovichS. A. (1931). Uber die Abgrenzung der Eigenwerte einer Matrix. Bulletin de l’Académie des Sciences de l’URSS 6 749–754.
  * [11] [author] Giordano, Ryan JR. J., Broderick, TamaraT. and Jordan, Michael IM. I. (2015). Linear response methods for accurate covariance estimates from mean field variational Bayes. Advances in Neural Information Processing Systems 28\.
  * [12] [author] Hadji, AmineA., Hesselink, TammoT. and Szabo, BotondB. (2022). Optimal recovery and uncertainty quantification for distributed Gaussian process regression. arXiv preprint arXiv:2205.03150.
  * [13] [author] Hadji, AmineA. and Szabó, BotondB. (2021). Can We Trust Bayesian Uncertainty Quantification from Gaussian Process Priors with Squared Exponential Covariance Kernel? SIAM/ASA Journal on Uncertainty Quantification 9 185–230.
  * [14] [author] Kimeldorf, George SG. S. and Wahba, GraceG. (1970). A correspondence between Bayesian estimation on stochastic processes and smoothing by splines. The Annals of Mathematical Statistics 41 495–502.
  * [15] [author] Knapik, BartekB., van der Vaart, Aad WA. W. and van Zanten, JHarryJ. (2011). Bayesian inverse problems with Gaussian priors. The Annals of Statistics 39 2626–2657.
  * [16] [author] Le Cam, LucienL. (2012). Asymptotic methods in statistical decision theory. Springer Science & Business Media.
  * [17] [author] Nieman, DennisD., Szabo, BotondB. and van Zanten, HarryH. (2022). Contraction rates for sparse variational approximations in Gaussian process regression. Journal of Machine Learning Research 23 1–26.
  * [18] [author] Rasmussen, Carl EdwardC. E. and Williams, Christopher KIC. K. (2006). Gaussian processes for machine learning. MIT press.
  * [19] [author] Ray, KolyanK. and Szabó, BotondB. (2022). Variational Bayes for high-dimensional linear regression with sparse priors. Journal of the American Statistical Association 117 1270–1281.
  * [20] [author] Ray, KolyanK., Szabo, BotondB. and Clara, GabrielG. (2020). Spike and slab variational Bayes for high dimensional logistic regression. Advances in Neural Information Processing Systems 33 14423–14434.
  * [21] [author] Rousseau, JudithJ. and Szabo, BotondB. (2020). Asymptotic frequentist coverage properties of Bayesian credible sets for sieve priors. The Annals of Statistics 48 2155–2179.
  * [22] [author] Schölkopf, BernhardB., Smola, Alexander JA. J., Bach, FrancisF. et al. (2002). Learning with kernels: support vector machines, regularization, optimization, and beyond. MIT press.
  * [23] [author] Srinivas, NiranjanN., Krause, AndreasA., Kakade, ShamS. and Seeger, MatthiasM. (2010). Gaussian Process Optimization in the Bandit Setting: No Regret and Experimental Design. International Conference on Machine Learning 1015–1022.
  * [24] [author] Steinwart, IngoI. and Scovel, ClintC. (2012). Mercer’s theorem on general domains: On the interaction between measures, kernels, and RKHSs. Constructive Approximation 35 363–417.
  * [25] [author] Szabó, BotondB., Van Der Vaart, Aad WA. W. and van Zanten, JHJ. (2015). Frequentist coverage of adaptive nonparametric Bayesian credible sets. The Annals of Statistics 43 1391–1428.
  * [26] [author] Titsias, MichalisM. (2009). Variational model selection for sparse Gaussian process regression. Report, University of Manchester, UK.
  * [27] [author] Vakili, SattarS., Scarlett, JonathanJ., Shiu, DashanD. and Bernacchia, AlbertoA. (2022). Improved Convergence Rates for Sparse Approximation Methods in Kernel-Based Learning. arXiv preprint arXiv:2202.04005.
  * [28] [author] Van der Vaart, Aad WA. W. (2000). Asymptotic Statistics. Cambridge University Press.
  * [29] [author] van der Vaart, Aad WA. W. and van Zanten, JHarryJ. (2007). Bayesian inference with rescaled Gaussian process priors. Electronic Journal of Statistics 1 433–448.
  * [30] [author] van der Vaart, Aad WA. W. and van Zanten, JHJ. (2008). Reproducing kernel Hilbert spaces of Gaussian priors. Pushing the Limits of Contemporary Statistics: Contributions in Honor of Jayanta K. Ghosh 200–222.
  * [31] [author] van der Vaart, Aad WA. W. and van Zanten, JHarryJ. (2008). Rates of contraction of posterior distributions based on Gaussian process priors. The Annals of Statistics 36 1435–1463.
  * [32] [author] Vivarelli, FrancescoF. and Opper, ManfredM. (1999). General bounds on Bayes errors for regression with Gaussian processes. Advances in Neural Information Processing Systems 11 302–308.
  * [33] [author] Wang, YixinY. and Blei, David MD. M. (2019). Frequentist consistency of variational Bayes. Journal of the American Statistical Association 114 1147–1161.
  * [34] [author] Wild, VeitV., Kanagawa, MotonobuM. and Sejdinovic, DinoD. (2021). Connections and Equivalences between the Nyström Method and Sparse Variational Gaussian Processes. arXiv preprint arXiv:2106.01121.
  * [35] [author] Yang, YunY., Pati, DebdeepD. and Bhattacharya, AnirbanA. (2020). $\alpha$-variational inference with statistical guarantees. The Annals of Statistics 48 886–905.
  * [36] [author] Zhang, FengshuoF. and Gao, ChaoC. (2020). Convergence rates of variational posterior distributions. The Annals of Statistics 48 2180–2207.