marginparsep has been altered.
topmargin has been altered.
marginparpush has been altered.
The page layout violates the ICML style. Please do not change the page layout,
or include packages like geometry, savetrees, or fullpage, which change it for
you. We’re not able to reliably undo arbitrary changes to the style. Please
remove the offending package(s), or layout-changing commands and try again.

Adaptive Cholesky Gaussian Processes

Simon Bartels 1 2  Kristoffer Stensbo-Smidt 2  Pablo Moreno-Muñoz 2  Wouter
Boomsma 1  Jes Frellsen 2  Søren Hauberg 2

††footnotetext:  1Datalogisk Institut, University of Copenhagen, København,
DENMARK 2DTU Compute, Technical University of Denmark, Kgs. Lyngby, DENMARK.
Correspondence to: Simon Bartels<EMAIL_ADDRESS>

###### Abstract

We present a method to fit exact Gaussian process models to large datasets by
considering only a subset of the data. Our approach is novel in that the size
of the subset is selected on the fly during exact inference with little
computational overhead. From an empirical observation that the log-marginal
likelihood often exhibits a linear trend once a sufficient subset of a dataset
has been observed, we conclude that many large datasets contain redundant
information that only slightly affects the posterior. Based on this, we
provide probabilistic bounds on the full model evidence that can identify such
subsets. Remarkably, these bounds are largely composed of terms that appear in
intermediate steps of the standard Cholesky decomposition, allowing us to
modify the algorithm to adaptively stop the decomposition once enough data
have been observed. Empirically, we show that our method can be directly
plugged into well-known inference schemes to fit exact Gaussian process models
to large datasets.

##  1 Introduction

Figure 1:  The figure shows total log-marginal likelihood as a function of the
size of the training set. The different colors correspond to different
Gaussian process models, using the SE kernel, Equation 14, with length scale
$\ell$ and amplitude $\sigma^{2}$. The inputs were sampled from a zero-mean
Gaussian with variance 100. The targets have mean $-2.5$ and variance
$25$.color=YellowGreen, inline, inline]Kristoffer: I don’t think I understand
how the data look like. It can be seen that the log-marginal likelihood
exhibits a linear trend after sufficiently many observations have been
processed.

It has been observed (Chalupka et al., 2013) that the subset-of-data
approximation can be a hard-to-beat baseline for approximate Gaussian process
inference. However, the question of how to choose the subset is non-trivial to
answer. Here we make an attempt.

The key computational challenge in Gaussian process regression is to evaluate
the log-marginal likelihood of the $N$ observed data points, which is known to
have cubic complexity (Rasmussen & Williams, 2006). In order to arrive at a
computationally less expensive approximation of this log-marginal likelihood,
we first empirically study its behavior as we increase the number of
observations. Figure 1 show this progression for a variety of models. We
elaborate on this figure later, but for now note that after a certain number
of observations, which differ between models, the log-marginal likelihood
starts to progress with a linear trend. This suggest that we may leverage this
near-linearity to estimate the log-marginal likelihood of the full dataset
after having seen only a subset of the data. However, as the point-of-
linearity differs between models, this must be estimated on-the-fly to keep
computations tractable.

In this paper, we approach the problem from a probabilistic numerics
perspective (Hennig et al., 2015). By treating the dataset as a collection of
independent and identically distributed random variables, we provide expected
upper and lower bounds on the log-marginal likelihood, which become tight when
the above-mentioned linear trend arises. We provide a particularly efficient
algorithm for computing the bounds that leverage the intermediate computations
performed by the Cholesky decomposition that is commonly used for evaluating
the log-marginal likelihood. The bounds are therefore practically free to
evaluate. We further show that these bounds allow us to predict _when_ the
linear trend determines the full-data log-marginal likelihood, such that we
can phrase an _optimal stopping problem_ to determine a suitable subset of the
data for a particular model. We call our method _Adaptive Cholesky Gaussian
Process_ (ACGP). Our approach has a complexity of $\mathcal{O}(M^{3})$, where
$M$ is the processed subset-size, inducing an overhead of $\mathcal{O}(M)$ to
the Cholesky. The main difference to previous work is that our algorithm does
_not necessarily_ look at the whole dataset.

##  2 Background

###  2.1 Notation

We use a python-inspired index notation, abbreviating for example
$[y_{1},\ldots,y_{n-1}]^{\top}$ as $\bm{y}_{:n}$—observe that the indexing
starts at 1. With $\operatorname{Diag}$ we define the operator that sets all
off-diagonal entries of a matrix to $0$.

###  2.2 Gaussian Process Regression

We start by briefly reviewing Gaussian process (GP) models and how they are
trained. A GP is a collection of random variables defined in terms of a mean
function, $m(\bm{x})$, and a covariance function or _kernel_ ,
$k(\bm{x},\bm{x}^{\prime})=\operatorname{cov}(f(\bm{x}),f(\bm{x}^{\prime}))$.
We denote such a GP as

$\displaystyle f\sim\mathcal{GP}(m(\bm{x}),k(\bm{x},\bm{x}^{\prime})).$

We consider the training dataset
$\mathcal{D}=\\{\bm{x}_{n},y_{n}\\}^{N}_{n=1}$ with inputs
$\bm{x}_{n}\in\mathbb{R}^{p}$ and outputs $y_{n}\in\mathbb{R}$. The inputs are
collected in the matrix
$\bm{X}=[\bm{x}_{1},\bm{x}_{2},\ldots,\bm{x}_{N}]^{\top}\in\mathbb{R}^{N\times
p}$ and their corresponding function values are denoted as
$\bm{f}=f(\bm{X})\in\mathbb{R}^{N}$. In a similar manner, we denote the test
inputs $\bm{X}_{*}$ and their corresponding function values $\bm{f}_{*}$. The
joint distribution over the training and test outputs, assuming a zero mean
function, $m(\bm{X})=\mathbf{0}$, is then

$\displaystyle p(\bm{f},\bm{f}_{*})=\mathcal{N}\Bigg{(}\begin{matrix}\bm{f}\\\
\bm{f}_{*}\end{matrix}\,\Bigg{|}\,\mathbf{0},\Bigg{[}\begin{matrix}\bm{K}_{\text{ff}}&\bm{K}_{\text{f}*}\\\
\bm{K}_{*\text{f}}&\bm{K}_{**}\end{matrix}\Bigg{]}\Bigg{)},$ (1)

where we have used the shorthand notation
$\bm{K}_{\text{ff}}=k(\bm{X},\bm{X})$ and similarly for the other covariance
matrices. As usual, we will consider the outputs $\bm{y}$ as being noise-
corrupted versions of the function values $\bm{f}$, and we shall parameterize
this corruption through the likelihood function
$p(\bm{y}\operatorname{|}\bm{f})$, which for regression tasks is typically
assumed to be Gaussian,
$p(\bm{y}\operatorname{|}\bm{f})=\mathcal{N}(\bm{f},\sigma^{2}\boldsymbol{I})$.
For such a model, the posterior process can be computed in closed-form:

$\displaystyle p(\bm{f}_{*}\operatorname{|}\bm{y})$
$\displaystyle=\mathcal{N}(\bm{m}_{*},\boldsymbol{S}_{*})$
$\displaystyle\bm{m}_{*}$ $\displaystyle=\bm{K}_{*\text{f}}\bm{K}^{-1}\bm{y}$
$\displaystyle\boldsymbol{S}_{*}$
$\displaystyle=\bm{K}_{**}-\bm{K}_{*\text{f}}\bm{K}^{-1}\bm{K}_{\text{f}*}.$

where $\bm{K}\colonequals\bm{K}_{\text{ff}}+\sigma^{2}\boldsymbol{I}$.

By marginalizing over the function values of the likelihood distribution, we
obtain the marginal likelihood, $p(\bm{y})$, which is typically used for
comparing the performance of models in the Bayesian framework:

$\displaystyle p(\bm{y})=\int
p(\bm{y}\operatorname{|}\bm{f})p(\bm{f})d\bm{f}.$ (2)

While this integral is not tractable in general, it does have a closed-form
solution for Gaussian process regression. Given the GP prior,
$p(\bm{f})=\mathcal{N}(\mathbf{0},\bm{K}_{\text{ff}})$, and the Gaussian
likelihood, the log-marginal likelihood distribution can be found to be

$\displaystyle\log
p(\bm{y})=-\frac{1}{2}\log|\bm{K}|-\frac{1}{2}\bm{y}^{\top}\bm{K}^{-1}\bm{y}-c,$
(3)

where $c=(N/2)\log 2\pi$.

###  2.3 Background on the Cholesky decomposition

Inverting covariance matrices such as $\bm{K}$ is a slow and numerically
unstable procedure. Therefore, in practice, one typically leverages the
Cholesky decomposition of the covariance matrices to compute the inverses. One
contribution of our work is to show that the Cholesky decomposition can be
stopped early to produce a highly efficient numerical scheme for Gaussian
process regression.

The Cholesky decomposition of a symmetric and positive definite matrix
$\boldsymbol{K}$ is the unique, lower111Equivalently, one can define
$\boldsymbol{L}$ to be upper triangular such that
$\boldsymbol{K}=\boldsymbol{L}^{\top}\boldsymbol{L}$. triangular matrix
$\boldsymbol{L}$ such that
$\boldsymbol{K}=\boldsymbol{L}\boldsymbol{L}^{\top}$ (Golub & Van Loan, 2013,
Theorem 4.2.7). The advantage of having such a decomposition is that inversion
with triangular matrices amounts to Gaussian elimination. There are different
options to compute $\boldsymbol{L}$. The Cholesky of a $1\times 1$ matrix is
the square root of the scalar. For larger matrices,

$\displaystyle\operatorname{chol}[\boldsymbol{K}]=\begin{bmatrix}\operatorname{chol}[\boldsymbol{K}_{:s,:s}]&\boldsymbol{0}\\\
\boldsymbol{T}&\operatorname{chol}\left[\boldsymbol{K}_{s:,s:}-\boldsymbol{T}\boldsymbol{T}^{\top}\right]\end{bmatrix},$
(4)

where
$\boldsymbol{T}\colonequals\boldsymbol{K}_{:s,s:}{\operatorname{chol}[\boldsymbol{K}_{:s,:s}]}^{-\top}$
and $s$ is any integer between $1$ and the size of $\boldsymbol{K}$.
color=YellowGreen, inline, ]Kristoffer: I would move the remark about the
python-inspired notation to here and the Diag convention to where it is used
the first time. Hence, extending a given Cholesky to a larger matrix requires
three steps:

  1. 1.

solve the linear equation system $\boldsymbol{T}$,

  2. 2.

apply the down date
$\boldsymbol{K}_{s:,s:}-\boldsymbol{T}\boldsymbol{T}^{\top}$ and

  3. 3.

compute the Cholesky of the down-dated matrix.

An important observation is that
$\boldsymbol{K}_{s:,s:}-\boldsymbol{T}\boldsymbol{T}^{\top}$ is the posterior
covariance matrix $\boldsymbol{S}_{*}+\sigma^{2}\boldsymbol{I}$ when
considering $\boldsymbol{X}_{s:}$ as test points. We will make use of this
observation in Section 3.5. The log-determinant of $\boldsymbol{K}$ can be
obtained from the Cholesky using
$\log|\boldsymbol{K}|=2\sum_{n=1}^{N}\log\boldsymbol{L}_{nn}$. A similar
recursive relationship exists between the quadratic form
$\boldsymbol{y}^{\top}\boldsymbol{K}^{-\\!1}\boldsymbol{y}$ and
$\boldsymbol{L}^{-\\!1}\boldsymbol{y}$ (see, e.g., Rasmussen & Williams
(2006)).

###  2.4 Related work

For a large number of observations $N$, the log-marginal likelihood becomes
computationally demanding to evaluate as the computation of the Cholesky
requires $\mathcal{O}(N^{3})$ operations and $\mathcal{O}(N^{2})$ storage.
Therefore much work has gone into tractable approximations to the log-marginal
likelihood, which is also the focus of our work. Arguably, the most popular
approximation methods for GPs are inducing point methods (Quiñonero-Candela &
Rasmussen, 2005; Snelson & Ghahramani, 2006; Titsias, 2009; Hensman et al.,
2013; 2015; 2017; Shi et al., 2020), where the dataset is approximated through
a set of pseudo-data points (inducing points), summarizing information from
nearby data. Other approaches involve building approximations to
$\boldsymbol{K}$ (Fine & Scheinberg, 2001; Harbrecht et al., 2012; Wilson &
Nickisch, 2015; Rudi et al., 2017; Wang et al., 2019) or aggregating of
distributed local approximations (Gal et al., 2014; Deisenroth & Ng, 2015).
One may also consider separately the approximation of the quadratic form via
linear solvers such as conjugate gradients (Hestenes & Stiefel, 1952; Cutajar
et al., 2016) and the approximation of the log-determinant (Fitzsimons et al.,
2017a; b; Dong et al., 2017). Another line of research is scaling the hardware
(Nguyen et al., 2019).

The work closest in spirit to the present paper is Artemev et al. (2021), who
also propose lower and upper bounds on quadratic form and log-determinant.
There are a number of differences, however. Their bound relies on the method
of conjugate gradients where we work directly with the Cholesky decomposition.
Furthermore, while their bounds are deterministic, ours are probabilistic,
which can make them tighter in certain cases, as they do not need to hold for
all worst-case scenarios.

All above referenced approaches have computational complexity at least
$\mathcal{O}(N)$. However, the size of a dataset is seldom a particularly
chosen value but rather the ad-hoc end of the sampling procedure. The
dependence on the dataset size implies that more available data allows to
spend less computional budget even though more data might not be helpful. This
is the main motivation for our work to derive an algorithm with sub-linear
complexity.

##  3 Methodology

In the following, we will sketch our method. Our main goal is to convey the
idea and intuition. To this end, we use suggestive notation. We refer the
reader to the appendix for a more thorough and formal treatment.

###  3.1 Intuition on the linear extrapolation

The marginal likelihood is typically presented as a joint distribution, but,
using Bayes rule, one can also view it from a cumulative perspective as the
sum of log-conditionals:

$\displaystyle\log p(\bm{y})$ $\displaystyle=\sum_{n=1}^{N}\log
p(y_{n}\operatorname{|}\bm{y}_{:n})\ .$ (5)

With this equation in hand, the phenomena in Figure 1 becomes much clearer.
The figure shows the value of Equation 5 for an increasing number of
observations $N$. When the plot exhibits a linear trend, it is because the
summands $\log p(y_{n}\operatorname{|}\bm{y}_{:n})$ become approximately
constant, implying that the model is not gaining additional knowledge. In
other words, new outputs are conditionally independent given the output
observations seen so far.inline]I’ve removed a footnote by Sören about the
spectra of the matrix. It had raised questions from a reviewer and it is not
entirely correct: the linear trend can also appear when the model is not
generalizing well. We should discuss this because the footnote might be
important for some readers.

The key problem addressed in this paper is how to estimate the full marginal
likelihood, $p(\bm{y})$, from only a subset of $M$ observations. The
cumulative view of the log-marginal likelihood in Equation 5 is our starting
point. color=YellowGreen, inline, ]Kristoffer: I don’t think the following two
sentences are needed anymore. In particular, we will provide bounds, which are
functions of seen observations, on the estimate of the full marginal
likelihood. These bounds will allow us to decide, on the fly, when we have
seen enough observations to accurately estimate the full marginal likelihood.

###  3.2 Stopping strategy

Suppose that we have processed $M$ data points with $(N-M)$ data points still
in the queue. We can then decompose Equation 5 into a sum of terms which have
already been computed and a remaining sum

$\displaystyle\log p(\boldsymbol{y})$
$\displaystyle=\underbrace{\sum_{n=1}^{M}\log
p(y_{n}\mid\boldsymbol{y}_{:n})}_{A:\text{
processed}}+\\!\\!\\!\underbrace{\sum_{n=M+1}^{N}\\!\\!\\!\log
p(y_{n}\mid\boldsymbol{y}_{:n})}_{B:\text{ remaining}}.$

Recall that we consider the $\boldsymbol{x}_{i},y_{i}$ as independent and
identically distributed random variables. Hence, we could estimate $B$ as
$A(N-M)/M$. color=YellowGreen, inline, ]Kristoffer: Would it be clearer if we
replace $A$ and $B$ with, say, $p(\boldsymbol{y}_{A})$ and
$p(\boldsymbol{y}_{B})$? My head keeps thinking of $A$ and $B$ as constants
like $N$ and $M$, which makes this discussion confusing. Yet this is estimator
is biased, since
$(\boldsymbol{x}_{M+1},y_{M+1}),\dots,(\boldsymbol{x}_{N},y_{N})$ interact
non-linearly through the kernel function. Instead, we will derive unbiased
lower and upper bounds, $\mathcal{L}$ and $\mathcal{U}$. To obtain unbiased
estimates, we use the last-$m$ processed points, such that conditioned on the
points up to $s\colonequals M-m$, $\log p(\boldsymbol{y})$ can be sandwiched,

$\displaystyle\mathbb{E}[\mathcal{L}\operatorname{|}\bm{X}_{:s},\bm{y}_{:s}]\leq
A+\mathbb{E}[B\operatorname{|}\bm{X}_{:s},\bm{y}_{:s}]\leq\mathbb{E}[\mathcal{U}\operatorname{|}\bm{X}_{:s},\bm{y}_{:s}],$

and the observations from $s$ to $M$ can be used to estimate $\mathcal{L}$ and
$\mathcal{U}$. We can then detect when the upper and lower bounds are
sufficiently near each other, and stop computations early when the
approximation is sufficiently good. More precisely, given a desired relative
error $r$, we stop when $\mathcal{U}$ and $\mathcal{L}$ have the same sign and
when

$\displaystyle\frac{\mathcal{U}-\mathcal{L}}{2\min(|\mathcal{U}|,|\mathcal{L}|)}<r.$
(6)

If the bounds hold, then the estimator
$\frac{1}{2}\mathcal{L}+\frac{1}{2}\mathcal{U}$ achieves the desired relative
error (Lemma 19 in appendix). This is in contrast to other approximations,
where one specifies a computational budget, rather than a desired accuracy.

The stopping conditions can be checked before or after Step 3 of the Cholesky
decomposition (Section 2.3). Here, we explore the former option since Step 3
is the bottleneck due to being less parallelizable than the other steps. As a
consequence, the term $A$ is available only until index $s$ and our bounds
estimate also the partially processed points. color=YellowGreen, inline,
]Kristoffer: not sure what this last sentence means.

###  3.3 Bounds

From Equation 3, we see that the log-marginal likelihood decomposes into the
log-determinant of the kernel matrix, a quadratic term, and a constant term.
In the following we present upper and lower bounds for both the log-
determinant ($\mathcal{U}_{\text{D}}$ and $\mathcal{L}_{\text{D}}$,
respectively) and the quadratic term ($\mathcal{U}_{\text{Q}}$ and
$\mathcal{L}_{\text{Q}}$).

We will need the posterior equations from Section 2 but for the observations,
that is $p(\boldsymbol{y}_{*}\mid\boldsymbol{y})$, and we will need them as
functions of $\boldsymbol{X}_{*}$. With a slight abuse of notation, denote
$\bm{m}_{*}^{(s)}(\cdot)$ for the posterior mean function conditioned on
$\boldsymbol{y}_{:s}$ and, analogously, define
$\bm{\Sigma}_{*}^{(s)}(\cdot,\cdot\cdot)$ as the posterior covariance function
conditioned on $\boldsymbol{y}_{:s}$. This allows to rewrite Equation 5 as

$\displaystyle\log p(\boldsymbol{y})$
$\displaystyle\propto\sum_{n=1}^{N}\log\bm{\Sigma}_{*}^{(n-1)}(\boldsymbol{x}_{n},\boldsymbol{x}_{n})$
(7)
$\displaystyle+\sum_{n=1}^{N}\frac{(y_{n}-\bm{m}_{*}^{(n-1)}(\boldsymbol{x}_{n}))^{2}}{\bm{\Sigma}_{*}^{(n-1)}(\boldsymbol{x}_{n},\boldsymbol{x}_{n})}$

which reveals that the log-determinant can be written as a sum of posterior
variances and the quadratic form has an expression as normalized square
errors.

Other key ingredients for our bounds are estimates for average posterior
variance and average correlation. Therefore define the short-hands

$\displaystyle\boldsymbol{V}$
$\displaystyle\colonequals\operatorname{Diag}\left[\boldsymbol{\Sigma}_{*}^{(s)}(\boldsymbol{X}_{s:M},\boldsymbol{X}_{s:M})\right]$
and $\displaystyle\boldsymbol{C}$
$\displaystyle\colonequals\sum_{i=1}^{\frac{M}{2}}\boldsymbol{\Sigma}_{*}^{(s)}(\boldsymbol{x}_{s+2i},\boldsymbol{x}_{s+2i-1})\boldsymbol{e}_{2i}\boldsymbol{e}_{2i}^{\top}\
,$

where $\boldsymbol{e}_{j}\in\mathbb{R}^{m}$ is the $j$-th standard basis
vector. The matrix $\boldsymbol{V}$ is simply the diagonal of the posterior
covariance matrix $\boldsymbol{\Sigma}_{*}$. The matrix $\boldsymbol{C}$
consists of every _second_ entry of the first off-diagonal of
$\boldsymbol{\Sigma}_{*}$. These elements are placed on the diagonal with
every second element being $0$. The reason for taking every second element is
of theoretical nature—see Remark 5 in the appendix. In practice we use the
full off-diagonal. color=YellowGreen, inline, ]Kristoffer: If you can use
either the full or a sparse off-diagonal, is it necessary to make these points
here given that we are only conveying the idea and intuition?

####  3.3.1 Bounds on the log-determinant

Both bounds, lower and upper, use that
$\log|\boldsymbol{K}|=\log|\boldsymbol{K}_{:s,:s}|+\log|\boldsymbol{\Sigma}_{*}^{(s)}(\boldsymbol{X}_{s:})|$
which follows from the matrix-determinant lemma. The first term is available
from the already processed datapoints. It is the second addend that needs to
be estimated, which we approach from the perspective of Equation 7. It is
well-established that, for a fixed input, more observations decrease the
posterior variance, and this decrease cannot cross the threshold $\sigma^{2}$
(Rasmussen & Williams, 2006, Question 2.9.4). This remains true when taking
the expectation over the input. Hence, the average of the posterior variances
for inputs $\boldsymbol{X}_{s:M}$ is with high probability an overestimate of
the average posterior variance for inputs with higher index. This motivates
our upper bound on the log-determinant:

$\displaystyle\mathcal{U}_{\text{D}}$
$\displaystyle=\log|\boldsymbol{K}_{:s,:s}|+(N-s)\mu_{D}$ (8)
$\displaystyle\mu_{D}$
$\displaystyle\colonequals\frac{1}{m}\boldsymbol{1}^{\top}\log\left(\boldsymbol{V}\right)\boldsymbol{1}$
$\displaystyle\sslash\text{{\emph{average log posterior variance}}}$

color=YellowGreen, inline, ]Kristoffer: The $\boldsymbol{1}$’s are hideous. I
would prefer to write $\sum\operatorname{diag}(\boldsymbol{V})$ or
$\sum\mathbb{V}[\boldsymbol{f}_{*}]$ or something like that.

To arrive at the lower bound on the log-determinant, we need an argument about
how fast the average posterior variance could decrease which is governed by
the correlation between inputs. The variable $\rho_{D}$ measures the average
correlation, and we show in Theorem 9 in the appendix that this overestimates
the decrease per step with high probability. Since the decrease cannot exceed
$\sigma^{2}$, we introduce $\psi_{D}$ to denote the step which would cross
this threshold.

$\displaystyle\mathcal{L}_{\text{D}}$
$\displaystyle=\log|\boldsymbol{K}_{:s,:s}|+(\psi_{D}-s)\left(\mu_{D}-\frac{\psi_{D}-s-1}{2\sigma^{4}}\rho_{D}\right)$
$\displaystyle\quad+(N-\psi_{D})\log\sigma^{2}.$ (9) $\displaystyle\rho_{D}$
$\displaystyle\colonequals\frac{2}{m}\boldsymbol{1}^{\top}\boldsymbol{C}\boldsymbol{C}\boldsymbol{1}$
$\displaystyle\sslash\text{{\emph{average square cross-correlation}}}$
$\displaystyle\psi_{D}$ $\displaystyle\colonequals\max(N,\lfloor
s-1+\frac{2}{\rho_{D}}\left(\mu_{D}-\log\sigma^{2}\right)\rfloor)$ (10)
$\displaystyle\sslash\text{{\emph{number of steps $\mu_{D}$ can decrease by
$\rho$ before hitting $\log\sigma^{2}$}}}.$

color=YellowGreen, inline, ]Kristoffer: Define $\lfloor\cdot\rfloor$. Both
bounds collapse to the exact solution when $s=N$. The bounds are close when
the average correlation between inputs, $\rho_{D}$, is small. This occurs for
example when the average variance is close to $\sigma^{2}$ since the variance
is an upper bound to the correlation, an example being shown in Figure 2.
color=YellowGreen, inline, ]Kristoffer: I don’t follow how the figure shows
this. Another case where $\rho_{D}$ is small is when points are not correlated
to begin with.

Figure 2: Prediction performance of ACGP vs. exact GP regression. The $5000$
inputs are distributed standard normal multiplied by a factor 5.
color=YellowGreen, inline, ]Kristoffer: Are you sure they are not just evenly
spaced? The targets stem from a periodic function with linear trend corrupted
by independent zero-mean Gaussian noise with variance $2.25$. The kernel is a
Matérn$\frac{5}{2}$ with parameters determined by marginal likelihood
maximization. ACGP stopped early using only a subset of $1536$ data points,
saving more than $50\%$ of the computational cost but getting similar results
as in the exact case.

####  3.3.2 Bounds on the quadratic term

Denote with
$\boldsymbol{e}_{*}\colonequals\boldsymbol{y}_{s:}-\boldsymbol{m}_{*}^{(s)}(\boldsymbol{X}_{s:})$
the prediction error, when considering the first $s$ points as training set
and the remaining inputs as test set. color=YellowGreen, inline, ]Kristoffer:
$\boldsymbol{e}$ has already been used! Could we use something else here or
redefine $\boldsymbol{C}$? Maybe use $\boldsymbol{r}$. Analogous to the bounds
on the log-determinant, one can show with the matrix inversion lemma that
$\boldsymbol{y}^{\top}\boldsymbol{K}^{-\\!1}\boldsymbol{y}=\boldsymbol{y}_{:s}^{\top}\boldsymbol{K}_{:s,:s}^{-\\!1}\boldsymbol{y}_{:s}+\boldsymbol{e}_{*}^{\top}\left(\boldsymbol{\Sigma}_{*}^{(s)}(\boldsymbol{X}_{s:})\right)^{-\\!1}\boldsymbol{e}_{*}$.
Again, the first term will turn out to be already computed. With a slight
abuse of notation let
$\boldsymbol{e}_{*}\colonequals\boldsymbol{y}_{s:M}-\boldsymbol{m}_{*}^{(s)}(\boldsymbol{X}_{s:M})$,
that is, we consider only the first $m$ entries. Our lower bound arises from
another well-known lower bound:
$\boldsymbol{a}^{\top}\boldsymbol{A}^{-\\!1}\boldsymbol{a}\geq
2\boldsymbol{a}^{\top}\boldsymbol{b}-\boldsymbol{b}^{\top}\boldsymbol{A}\boldsymbol{b}$
for all $\boldsymbol{b}$ (see for example Kim & Teh (2018); Artemev et al.
(2021)). We choose $\boldsymbol{b}\colonequals\alpha\boldsymbol{1}$ where
$\alpha$ is chosen to maximize the bound. The result, after some
cancellations, is the following lower bound on the quadratic term:

$\displaystyle\mathcal{L}_{\text{Q}}$
$\displaystyle=\boldsymbol{y}_{:s}^{\top}\boldsymbol{K}_{:s,:s}^{-\\!1}\boldsymbol{y}_{:s}+(N-s)\alpha\left(2\mu_{Q}-\alpha\rho_{Q}\right)$
(11) $\displaystyle\mu_{Q}$
$\displaystyle\colonequals\frac{1}{m}\boldsymbol{e}_{*}^{\top}\boldsymbol{e}_{*}$
$\displaystyle\sslash\text{{\emph{average square error}}}$
$\displaystyle\rho_{Q}$
$\displaystyle\colonequals\frac{1}{m}\boldsymbol{e}_{*}^{\top}\boldsymbol{e}_{*}+\frac{N-s-1}{m}\sum_{j=\frac{s+2}{2}}^{\frac{M}{2}}\boldsymbol{e}_{*,2j}\boldsymbol{e}_{*,2j-1}\boldsymbol{C}_{2j,2j}$

The $\alpha$ maximizing above bound is $\frac{\mu_{Q}}{\rho_{Q}^{2}}$, which
is the value we chose in our implementation. Note however, that
$\mathcal{L}_{Q}$ is an expected lower bound only if $\alpha$ depends on
variables with index smaller than $s$.

Our upper bound arises from the element-wise perspective of Equation 7. We
assume that the expected mean square error
$(y_{n}-\bm{m}_{*}^{(n-1)}(\boldsymbol{x}_{n}))^{2}$ decreases with more
observations. However, though mean square error and variance decrease, their
expected ratio may increase or decrease depending on the choice of kernel,
dataset and number of processed points. Using the average error calibration
with a correction for the decreasing variance, we arrive at our upper bound on
the quadratic term:

$\displaystyle\mathcal{U}_{\text{Q}}$
$\displaystyle=\boldsymbol{y}_{:s}^{\top}\boldsymbol{K}_{:s,:s}^{-\\!1}\boldsymbol{y}_{:s}+(N-s)\left(\mu_{Q}^{\prime}+\rho_{Q}^{\prime}\right)$
(12) $\displaystyle\mu_{Q}^{\prime}$
$\displaystyle\colonequals\frac{1}{m}\boldsymbol{e}_{*}^{\top}\boldsymbol{V}^{-\\!1}\boldsymbol{e}_{*}$
$\displaystyle\sslash\text{{\emph{average error calibration}}}$
$\displaystyle\rho_{Q}^{\prime}$
$\displaystyle\colonequals\frac{N-s-1}{m}\frac{1}{\sigma^{4}}\boldsymbol{e}_{*}^{\top}\boldsymbol{C}\boldsymbol{V}^{-\\!1}\boldsymbol{C}\boldsymbol{e}_{*}$
$\displaystyle\sslash\text{{\emph{average increase in error calibration}}}$

In our implementation we use a slightly different upper bound. The estimate of
the possible decrease of the variance uses the same technique as the lower
bound for the log-determinant. Therefore we can define an analogue to Equation
10 determining the step when the variance estimate falls below $\sigma^{2}$.
In our implementation, addends of the quadratic after this step are estimated
by the more conservative $\sigma^{-2}\mu_{Q}$.

Again, the bounds collapse to the true quantity when $s=N$. The bounds will
give good estimates when the average correlation between inputs, represented
by the matrix $\boldsymbol{C}$, is low or when the model can predict new data
well, that is, when $\boldsymbol{e}_{*}$ is close to $0$.

###  3.4 Result

color=YellowGreen, inline, ]Kristoffer: Result on what? I like “Main
theoretical results” or perhaps something like “Theoretical guarantees of the
bounds”?

######  Assumption 1.

Assume that

$\displaystyle\mathbb{E}\left[f(\boldsymbol{x},\boldsymbol{x}^{\prime})(y_{j}-\boldsymbol{m}_{*}^{(j-1)}(\boldsymbol{x}))^{2}\mid\boldsymbol{X}_{:s},\boldsymbol{y}_{:s}\right]$
$\displaystyle\leq\mathbb{E}\left[f(\boldsymbol{x},\boldsymbol{x}^{\prime})(y_{j}-\boldsymbol{m}_{*}^{(s)}(\boldsymbol{x}))^{2}\mid\boldsymbol{X}_{:s},\boldsymbol{y}_{:s}\right]$

for all $s\in\\{1,\dots,N\\}$ and for all $s<j\leq N$, where
$f(\boldsymbol{x},\boldsymbol{x}^{\prime})$ is either
$\frac{1}{\boldsymbol{\Sigma}_{*}^{(s)}(\boldsymbol{x},\boldsymbol{x})}$ or
$\frac{\boldsymbol{\Sigma}_{*}^{(s)}(\boldsymbol{x},\boldsymbol{x}^{\prime})^{2}}{\sigma^{4}\boldsymbol{\Sigma}_{*}^{(s)}(\boldsymbol{x},\boldsymbol{x})}$.

This assumption is required for the upper bound on the quadratic form, and it
expresses the intuition that the mean square error should decrease with more
data, which, empirically, appears to be true.

######  Theorem 2.

Assume that $(\boldsymbol{x}_{1},y_{1}),\dots,(\boldsymbol{x}_{N},y_{N})$ are
independent and identically distributed, assume that 1 holds, and assume that
$\alpha$ in the definition of $\mathcal{L}_{Q}$ depends only on
$\boldsymbol{x}_{1},y_{1},\dots\boldsymbol{x}_{s},y_{s}$.

For any $s\in\\{1,\dots,N\\}$, the bounds defined in Equations 9, 8, 12 and 11
hold in expectation:

$\displaystyle\mathbb{E}[\mathcal{L}_{D}\mid\mathcal{F}_{s}]\leq\mathbb{E}[\log(\det[\boldsymbol{K}])\mid\mathcal{F}_{s}]\leq\mathbb{E}[\mathcal{U}_{D}\mid\mathcal{F}_{s}]\text{
and}$
$\displaystyle\mathbb{E}[\mathcal{L}_{Q}\mid\mathcal{F}_{s}]\leq\mathbb{E}[\boldsymbol{y}^{\top}\boldsymbol{K}^{-\\!1}\boldsymbol{y}\mid\mathcal{F}_{s}]\leq\mathbb{E}[\mathcal{U}_{Q}\mid\mathcal{F}_{s}]\
,$

where $\mathcal{F}_{s}$ is the $\sigma$-algebra generated by
$\boldsymbol{x}_{1},y_{1},\dots,\boldsymbol{x}_{s},y_{s}$.

Proof and a proof sketch can be found in the appendix.

######  Theorem 3.

Let $r>0$ be a desired relative error and set
$\mathcal{U}\colonequals\mathcal{U}_{D}+\mathcal{U}_{Q}$ and
$\mathcal{L}\colonequals\mathcal{L}_{D}+\mathcal{L}_{Q}$. If the stopping
conditions hold, that is,
$\operatorname{sign}(\mathcal{U})=\operatorname{sign}(\mathcal{L})$ and
Equation 6 is true, then $\log p(\boldsymbol{y})$ can be estimated from
$\frac{1}{2}(\mathcal{U}+\mathcal{L})$ such that, under the condition
$\mathcal{L}_{D}\leq\log(\det[\boldsymbol{K}])\leq\mathcal{U}_{D}\text{ and
}\mathcal{L}_{Q}\leq\boldsymbol{y}^{\top}\boldsymbol{K}^{-\\!1}\boldsymbol{y}\leq\mathcal{U}_{Q}$,
the relative error is smaller than $r$, formally:

$\displaystyle\left|\frac{\log
p(\boldsymbol{y})-\frac{1}{2}(\mathcal{U}+\mathcal{L})}{\log
p(\boldsymbol{y})}\right|\leq r\ .$ (13)

The proof follows from Lemma 19 in the appendix.

Theorem 2 is a first step to obtain a probabilistic statement for Equation 13,
that is, a statement of the form $\mathbb{P}\left(\left|\frac{\log
p(\boldsymbol{y})-\frac{1}{2}(\mathcal{U}+\epsilon_{\mathcal{U},\delta}+\mathcal{L}-\epsilon_{\mathcal{L},\delta})}{\log
p(\boldsymbol{y})}\right|>r\right)\leq\delta$. Theoretically, we can obtain
such a statement using standard concentration inequalities and a union bound
over $s$. In practice, the error guarding constants $\epsilon$ would render
the result trivial. A union bound can be avoided using Hoeffding’s inequality
for martingales (Fan et al., 2012). However, this requires to replace
$s\colonequals M-m$ by a stopping time independent of $M$, which we regard as
future work.

###  3.5 Practical implementation

The proposed bounds turn out to be surprisingly cheap to compute. If we set
the block-size of the Cholesky decomposition to be $m$, the matrix
$\boldsymbol{\Sigma}_{*}^{(s)}$ is exactly the down-dated matrix in step 2 of
the Cholesky decomposition algorithm outlined in Section 2.3. Similarly, the
expressions for the bounds on the quadratic form appear while solving the
linear equation system $\boldsymbol{L}^{-\\!1}\boldsymbol{y}$. A slight
modification to the Cholesky algorithm is enough to compute these bounds on
the fly during the decomposition with little overhead.

Note that the definition of the bounds does not involve variables
$\boldsymbol{x},y$ which have not been processed. This allows an on-the-fly
construction of the kernel matrix, avoiding potentially expensive kernel
function evaluations. Furthermore, it is _not_ necessary to allocate
$\mathcal{O}(N^{2})$ memory in advance—a user can specify a maximal amount of
processed datapoints, hoping that stopping occurs before hitting that limit.
We provide the pseudo-code for this modified algorithm, our key algorithmic
contribution, in supplementary. Additionally, we provide a Python
implementation of our modified Cholesky decomposition and scripts to replicate
the experiments of this paper.222The code is available at the following
repository: anonymized inline]we must not forget to put an actual URL here

##  4 Experiments

In this section we examine bounds and stopping strategy for ACGP. When running
experiments without GPU support, all linear algebra operations are substituted
for direct calls to the openblas library (Wang et al., 2013), for efficient
realization of in-place operations. To still benefit from automatic
differentiation, we used pytorch (Paszke et al., 2019) with a custom backward
function for $\log p(\boldsymbol{y})$. We found that evaluating gradients with
pytorch costs roughly a factor six of the standard evaluation of the log-
marginal. We were able to reduce this computation time to a factor of two.
color=YellowGreen, inline, ]Kristoffer: Maybe add “but it will still affect
the timings in the following sections compared to the OpenBLAS implementation”
or something like that. Otherwise, it is unclear why it is mentioned.

###  4.1 Performance on synthetic data

Figure 2 shows an evaluation of the performance of ACGP in an illustrative
manner with synthetic data. The dataset consists of $5000$ inputs distributed
standard normal multiplied by a factor 5. The targets stem from a periodic
function with linear trend corrupted by independent zero-mean Gaussian noise
with variance $2.25$. The kernel is a Matérn$\frac{5}{2}$ with parameters
determined by marginal likelihood maximization. With $m\colonequals 256$ and
$r\colonequals 0.1$, ACGP uses only $1536$ data points, saving more than
$50\%$ of the computational cost but getting similar results as in the exact
case. color=YellowGreen, inline, ]Kristoffer: This section has been copied to
the caption of Figure 2. It seems weird to have it here again.

###  4.2 Bound quality

In this section we examine the bounds presented in Section 3 and compare them
to those proposed by Artemev et al. (2021, Lemma 2 and Lemma 3). Specifically,
for the determinant we compare to their $\mathcal{O}(N)$ upper bound (Artemev
et al., 2021, Eq. 11) and $\log(\det[\boldsymbol{Q}])$ as lower bound.

We set the number of inducing inputs $M$ for CGLB to 512, 1024 and 2048. For
ACGP, we define $m{}\colonequals 40\times 256=10240$ which is the number of
cores times the default openblas block size for our machines. We compare both
methods using squared exponential kernel (SE) and the Ornstein-Uhlenbeck
kernel (OU)

$\displaystyle k_{\text{SE}}(\boldsymbol{x},\boldsymbol{z})$
$\displaystyle\colonequals\theta\exp\left(-\frac{\|\boldsymbol{x}-\boldsymbol{z}\|^{2}}{2\ell^{2}}\right)$
(14) $\displaystyle k_{\text{OU}}(\boldsymbol{x},\boldsymbol{z})$
$\displaystyle\colonequals\theta\exp\left(-\frac{\|\boldsymbol{x}-\boldsymbol{z}\|}{\ell}\right).$
(15)

where we fix $\sigma^{2}\colonequals 10^{-3}$ and $\theta\colonequals 1$, and
we vary $\ell$ as $\log\ell\in\\{-1,0,1,2,3\\}$. As benchmarking datasets we
use the two datasets consisting of more than 20 000 instances used by Artemev
et al. (2021): kin40k and protein. We further consider two additional datasets
from the UCI repository (Dua & Graff, 2019):
metro333http://archive.ics.uci.edu/ml/datasets/Metro+Interstate+Traffic+Volume
no citation request and
pm25444http://archive.ics.uci.edu/ml/datasets/Beijing+PM2.5+Data (Liang et
al., 2015). We chose these datasets in addition as they are of similar size,
they are marked as regression tasks and no data points are missing.

Figure 3: Upper and lower bounds on the quadratic term for ACGP and CGLB on
the protein dataset using the OU kernel, Equation 15, with a length scale of
$\log\ell=-1$. The black line indicates the result obtained using exact GP
regression with points above it marking the upper bounds and points below
marking the lower bounds. The experiment was repeated a number of times with
different seeds to illustrate the variability in the computation time, shown
here as multiple points. For ACGP the number near the points shows $M$, the
size of the used subset; for CGLB it is the number of inducing inputs. The
color of the points reflects these numbers to help discern the size of the
subset or number of inducing inputs from the repeated experiments.  Figure 4:
Upper and lower bounds on the log-determinant for ACGP and CGLB on the pm25
dataset using the SE kernel, Equation 14, with a length scale of $\log\ell=0$.

The typical picture we observe is that CGLB provides better estimates for the
quadratic term, whereas ACGP is faster to identify the log-determinant.
Figures 3 and 4 show typical examples. Note that Figure 3 is on a log-scale:
for the quadratic form, the upper bounds tend to be less tight than the lower
lower bounds. Generally, there is no clear winner—sometimes ACGP estimates
both quantities faster and sometimes CGLB. For other results, see the
appendix.

The reason why CGLB has more difficulties to approximate the log-determinant
is that the bound involves
$\operatorname{trace}[\boldsymbol{K}-\boldsymbol{Q}]$ where $\boldsymbol{Q}$
is a low rank approximation to $\boldsymbol{K}$. If $\bm{K}_{\text{ff}}$ is of
high rank, the gap in trace can be large. For CGLB the time to compute the
bounds is dominated by the pivoted Cholesky decomposition to select the
inducing inputs. This overhead becomes irrelevant for the following hyper-
parameter tuning experiments, since the selection is computed only once in the
beginning. One conclusion from these experiments is to keep in mind that when
high precision is required, simply computing the exact solution can be a hard
to beat baseline.

###  4.3 Application in hyper-parameter tuning

We repeat the hyper-parameter tuning experiments performed by Artemev et al.
(2021) using the same set-up.We use the same kernel function, a
Matérn$\frac{3}{2}$, and the same optimizer: L-BFGS-B (Liu & Nocedal, 1989)
with SciPy (Virtanen et al., 2020) default parameters if not specified
otherwise. Each configuration has been measured five times with a different
shuffle of the dataset. All algorithms are stopped the latest after 2000
optimization steps, after 12 hours of compute time or when optimization has
failed three times. For CGLB and the exact Cholesky, the L-BFGS-B convergence
criterion “relative change in function value” (ftol) is set to 0. For ACGP, we
successively decrease this tolerance as $(2/3)^{\text{restart}+1}$ and we set
the same value for $r$. Furthermore, we use a block size of $m\colonequals
40\cdot 256=10240$ where $256$ is the standard openblas blocksize on our
machines and $40$ is the number of cores. Artemev et al. (2021) report their
best results using $M=2048$ inducing inputs. We refer to _op. cit._ for a
comparison to more baselines.

We explore two different computing environments. For datasets smaller than 20
000 data points, we ran our experiments on a single GPU. This is the same
setup as in Artemev et al. (2021) with the difference that we use a Titan RTX
whereas they have used a Tesla V100. The results can be summarized in one
paragraph: all methods converge the latest after two minutes. The time
difference between methods is less than twenty seconds. Exact Gaussian process
regression is fastest, more often than not. Figure 5 shows exemplary the
result for the elevators dataset. All other results can be found in the
appendix. We conclude that in an environment with significantly more
processing resources than memory, approximation may just cause overhead.

Figure 5:  The evolution of the root mean square error while tuning parameters
on elevators. When training on a GPU and the kernel matrix fits fully into
memory, exact inference can be a hard-to-beat baseline.

For datasets larger than 20 000 datapoints, our setup differs from Artemev et
al. (2021). We use only CPUs on machines where the kernel matrix still fits
fully into memory. Specifically, we used machines running Ubuntu 18.04 with 50
Gigabytes of RAM and two Intel Xeon E5-2670 v2 CPUs. The datasets are the same
as in Section 4.2. On all datasets, ACGP is essentially exhibiting the same
optimization behavior as the exact Gaussian process regressor, just stretched
out. ACGP can provide results faster than exact optimization but may be slower
in convergence as Figure 6 shows for the protein data-set. This observation is
as expected. However, approximation can also hinder fast convergence as Figure
7 reveals on for the metro dataset. CGLB benefits from caching the chosen
inducing inputs and reusing the solution from the last solved linear equation
system. The algorithm is faster, though it often plateaus at higher objective
function values. The results for kin40k are similar to protein and the results
for pm25 are similar to metro. These results and the evolution of the root
mean square error over time can be found in the appendix. Again, when the
available memory permits, the exact computation is a hard-to-beat baseline.
However, the Cholesky as a standard numerical routine has been engineered over
decades, whereas for the implementations of CGLB and ACGP there is opportunity
for improvement.

Figure 6: The evolution of the exact log marginal likelihood
$p(\boldsymbol{y})$ while tuning parameters for the protein dataset. The
iteratively increasing precision may allow ACGP to provide better solutions
faster than exact optimization at the price of later convergence.  Figure 7:
The evolution of the exact log marginal likelihood $p(\boldsymbol{y})$ while
tuning parameters for the metro dataset. Function evaluations with CGLB are
generally the fastest at the cost of plateauing at higher objective function
values.

##  5 Conclusions

In this paper we have revisited the use of Cholesky decompositions in Gaussian
process regression. We have shown that the Cholesky decomposition almost
computes expected lower and upper bounds on the marginal log-likelihood
associated with GP regression. With fairly limited modifications to this
classic matrix decomposition, we can use these bounds to stop the
decomposition before all observations have been processed. This has the
practical benefit that the kernel matrix $\boldsymbol{K}$ does not have to
computed prior to performing the decomposition, but can rather be computed on-
the-fly.

Empirical results indicate that the approach carries significant promise, but
no clear ‘winner’ can be crowned from our experiments. In general, we find
that exact GP inference leads to better behaved optimization than
approximations such as CGLB and inducing point methods, and that a well-
optimized Cholesky implementation is surprisingly competitive in terms of
performance. An advantage of our approach is that it is essentially parameter-
free. The user has to specify a requested numerical accuracy and the
computational demands will be scaled accordingly. Finally, we note that ACGP
is complementary to much existing work, and should be seen as an addition to
the GP toolbox, rather than a substitute for all existing tools.

#### Acknowledgements

We are grateful for the valuable feedback of all anonymous reviewers who saw
the different iterations of this article.

This work was funded in part by the Novo Nordisk Foundation through the Center
for Basic Machine Learning Research in Life Science (NNF20OC0062606). It also
received funding from the European Research Council (ERC) under the European
Union’s Horizon 2020 research, innovation programme (757360), from a research
grant (15334) from VILLUM FONDEN, from the Danish Ministry of Education and
Science, and from Digital Pilot Hub and Skylab Digital.

## References

  * Artemev et al. (2021) Artemev, A., Burt, D. R., and van der Wilk, M.  Tighter bounds on the log marginal likelihood of gaussian process regression using conjugate gradients.  In Meila, M. and Zhang, T. (eds.), _Proceedings of the 38th International Conference on Machine Learning_ , volume 139 of _Proceedings of Machine Learning Research_ , pp. 362–372, 2021.
  * Bartels (2020) Bartels, S.  _Probabilistic Linear Algebra_.  PhD thesis, University of Tübingen, 2020.
  * Bartels et al. (2021) Bartels, S., Boomsma, W., Frellsen, J., and Garreau, D.  Kernel-Matrix Determinant Estimates from stopped Cholesky Decomposition.  _ArXiv e-prints_ , 2107.10587, 2021.  2107.10587.
  * Chalupka et al. (2013) Chalupka, K., Williams, C. K. I., and Murray, I.  A framework for evaluating approximation methods for Gaussian process regression.  _Journal of Machine Learning Research_ , 14(1):333–350, 2013.
  * Cutajar et al. (2016) Cutajar, K., Osborne, M., Cunningham, J., and Filippone, M.  Preconditioning kernel matrices.  In Balcan, M. F. and Weinberger, K. Q. (eds.), _Proceedings of The 33rd International Conference on Machine Learning_ , volume 48 of _Proceedings of Machine Learning Research_ , pp. 2529–2538, 2016.
  * Deisenroth & Ng (2015) Deisenroth, M. and Ng, J. W.  Distributed Gaussian processes.  In _International Conference on Machine Learning (ICML)_ , pp. 1481–1490, 2015.
  * Dong et al. (2017) Dong, K., Eriksson, D., Nickisch, H., Bindel, D., and Wilson, A. G.  Scalable log determinants for gaussian process kernel learning.  In _Advances in Neural Information Processing Systems_ , pp. 6330–6340, 2017.
  * Dua & Graff (2019) Dua, D. and Graff, C.  UCI machine learning repository, 2019.  URL http://archive.ics.uci.edu/ml.
  * Fan et al. (2012) Fan, X., Grama, I., and Liu, Q.  Hoeffding’s inequality for supermartingales.  _Stochastic Processes and their Applications_ , 122(10):3545–3559, 2012.
  * Fine & Scheinberg (2001) Fine, S. and Scheinberg, K.  Efficient svm training using low-rank kernel representations.  _Journal of Machine Learning Research_ , 2:243–264, 2001\.
  * Fitzsimons et al. (2017a) Fitzsimons, J., Cutajar, K., Osborne, M., Roberts, S., and Filippone, M.  Bayesian inference of log determinants.  In Elidan, G., Kersting, K., and Ihler, A. T. (eds.), _Thirty-Third Conference on Uncertainty in Artificial Intelligence, UAI 2017, August 11-15, 2017, Sydney, Australia_ , 2017a.
  * Fitzsimons et al. (2017b) Fitzsimons, J., Granziol, D., Cutajar, K., Osborne, M., Filippone, M., and Roberts, S.  Entropic trace estimates for log determinants.  In Ceci, M., Hollmén, J., Todorovski, L., Vens, C., and Džeroski, S. (eds.), _Machine Learning and Knowledge Discovery in Databases_ , pp. 323–338, 2017b.
  * Gal et al. (2014) Gal, Y., Van Der Wilk, M., and Rasmussen, C. E.  Distributed variational inference in sparse gaussian process regression and latent variable models.  _arXiv preprint arXiv:1402.1389_ , 2014.
  * George et al. (1986) George, A., Heath, M. T., and Liu, J.  Parallel cholesky factorization on a shared-memory multiprocessor.  _Linear Algebra and its Applications_ , 77:165–187, 1986\.
  * Golub & Van Loan (2013) Golub, G. and Van Loan, C.  _Matrix computations_.  Johns Hopkins Univ Pr, 4 edition, 2013.
  * Harbrecht et al. (2012) Harbrecht, H., Peters, M., and Schneider, R.  On the low-rank approximation by the pivoted Cholesky decomposition.  _Applied Numerical Mathematics_ , 62(4):428–440, 2012.
  * Hennig et al. (2015) Hennig, P., Osborne, M., and Girolami, M.  Probabilistic numerics and uncertainty in computations.  _Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences_ , 471(2179), 2015.
  * Hensman et al. (2013) Hensman, J., Fusi, N., and Lawrence, N. D.  Gaussian processes for big data.  In _Uncertainty in Artificial Intelligence (UAI)_ , pp. 282–290, 2013.
  * Hensman et al. (2015) Hensman, J., Matthews, A., and Ghahramani, Z.  Scalable variational Gaussian process classification.  In _Artificial Intelligence and Statistics (AISTATS)_ , pp. 351–360, 2015.
  * Hensman et al. (2017) Hensman, J., Durrande, N., Solin, A., et al.  Variational fourier features for gaussian processes.  _J. Mach. Learn. Res._ , 18(1):5537–5588, 2017\.
  * Hestenes & Stiefel (1952) Hestenes, M. and Stiefel, E.  Methods of conjugate gradients for solving linear systems.  _Journal of Research of the National Bureau of Standards_ , 49(6):409–436, 1952.
  * Kim & Teh (2018) Kim, H. and Teh, Y. W.  Scaling up the automatic statistician: Scalable structure discovery using gaussian processes.  In Storkey, A. and Perez-Cruz, F. (eds.), _Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics_ , volume 84 of _Proceedings of Machine Learning Research_ , pp. 575–584, 2018.
  * Liang et al. (2015) Liang, X., Zou, T., Guo, B., Li, S., Zhang, H., Zhang, S., Huang, H., and Chen, S. X.  Assessing beijing’s $pm_{2.5}$ pollution: severity, weather impact, apec and winter heating.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 471(2182):20150257, 2015.
  * Liu & Nocedal (1989) Liu, D. C. and Nocedal, J.  On the limited memory bfgs method for large scale optimization.  _Mathematical Programming_ , 45(1):503–528, 1989\.
  * Mnih et al. (2008) Mnih, V., Szepesvári, C., and Audibert, J.-Y.  Empirical Bernstein stopping.  pp. 672–679, 2008.
  * Nguyen et al. (2019) Nguyen, D.-T., Filippone, M., and Michiardi, P.  Exact gaussian process regression with distributed computations.  In _Proceedings of the 34th ACM/SIGAPP Symposium on Applied Computing_ , pp. 1286–1295, 2019.
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S.  Pytorch: An imperative style, high-performance deep learning library.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems 32_ , pp. 8024–8035. 2019.
  * Quiñonero-Candela & Rasmussen (2005) Quiñonero-Candela, J. and Rasmussen, C.  A unifying view of sparse approximate Gaussian process regression.  _J of Machine Learning Research_ , 6:1939–1959, 2005.
  * Rasmussen & Williams (2006) Rasmussen, C. and Williams, C.  _Gaussian Processes for Machine Learning_.  MIT, 2006.
  * Rudi et al. (2017) Rudi, A., Carratino, L., and Rosasco, L.  Falkon: An optimal large scale kernel method.  In Guyon, I., Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 30, 2017.
  * Shi et al. (2020) Shi, J., Titsias, M., and Mnih, A.  Sparse orthogonal variational inference for gaussian processes.  In _International Conference on Artificial Intelligence and Statistics_ , pp. 1932–1942. PMLR, 2020.
  * Snelson & Ghahramani (2006) Snelson, E. and Ghahramani, Z.  Sparse gaussian processes using pseudo-inputs.  _Advances in neural information processing systems_ , 18:1257, 2006.
  * Titsias (2009) Titsias, M.  Variational learning of inducing variables in sparse gaussian processes.  In _Artificial intelligence and statistics_ , pp. 567–574. PMLR, 2009.
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., Carey, C. J., Polat, İ., Feng, Y., Moore, E. W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E. A., Harris, C. R., Archibald, A. M., Ribeiro, A. H., Pedregosa, F., van Mulbregt, P., and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  _Nature Methods_ , 17:261–272, 2020.
  * Wang et al. (2019) Wang, K., Pleiss, G., Gardner, J., Tyree, S., Weinberger, K. Q., and Wilson, A. G.  Exact gaussian processes on a million data points.  _Advances in Neural Information Processing Systems_ , 32:14648–14659, 2019.
  * Wang et al. (2013) Wang, Q., Zhang, X., Zhang, Y., and Yi, Q.  AUGEM: Automatically generate high performance Dense Linear Algebra kernels on x86 CPUs.  In _SC ’13: Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis_ , pp. 1–12, 2013.
  * Wilson & Nickisch (2015) Wilson, A. and Nickisch, H.  Kernel interpolation for scalable structured gaussian processes (kiss-gp).  In _International Conference on Machine Learning_ , pp. 1775–1784. PMLR, 2015.

##  Appendix A ADDITIONAL RESULTS

In this section, we report additional results for both the hyper-parameter
tuning experiments (section A.1) as well as plots to show the quality of the
bounds on both the log-determinant term and the quadratic term (sections
A.2–A.4).

##### Hyper-parameter tuning.

We repeat the experiments performed by Artemev et al. (2021). They randomly
split each dataset into a training set consisting of 2/3 of examples, and a
test set consisting of the remaining third. Inputs and targets are normalized
to have mean 0 and variance 1 within the training set. We use the same kernel
function, a Matérn$\frac{3}{2}$, and the same optimizer. For ACGP, we
successively decrease the optimizer’s tolerance as $(2/3)^{\text{restart}+1}$
and we set the same value for $r$. The block size is set to
40$\cdot$256=10192. Artemev et al. (2021) report their best results using
$M=2048$ inducing inputs.

We explore two different computing environments. For datasets smaller than 20
000 data points, we ran our experiments on a single GPU. This is the same
setup as in Artemev et al. (2021) with the difference that we use a Titan RTX
whereas they have used a Tesla V100. For datasets larger than 20 000
datapoints, our setup differs from Artemev et al. (2021). We use only CPUs on
machines where the kernel matrix still fits fully into memory. Specifically,
we used machines running Ubuntu 18.04 with 50 Gigabytes of RAM and two Intel
Xeon E5-2670 v2 CPUs.

The plots for the hyper-parameter optimization are shown in figures 8–19. Each
point in the plots corresponds to one accepted optimization step for the given
methods. Each point thus corresponds to a particular set of hyper-parameters
during the optimization. In figures 12–19, we show the root-mean-square error,
RMSE, that each methods obtains on the test set at each optimisation step. In
figures 8–11, we show the log-marginal likelihood, $\log p(\bm{y})$, that an
exact GP would have achieved with the specific set of hyper-parameters at each
optimization step for each method.

##### Bound quality plots.

For the bound quality plots, we compare ACGP to CGLB proposed by Artemev et
al. (2021). For CGLB, we compute the bounds with varying number of inducing
inputs $M\colonequals\\{512,1024,2048\\}$ and measure the time it takes to
compute the bounds. The three values for $M$ are connected with a line in the
plots. For ACGP, we define $m{}\colonequals 10192$. We measure the elapsed
time every time a block of data points is added to the processed dataset and
the bounds are recomputed.

We compare both methods using squared exponential kernel (SE) and the
Ornstein-Uhlenbeck kernel (OU).

$\displaystyle k_{SE}(\boldsymbol{x},\boldsymbol{z})$
$\displaystyle\colonequals\theta\exp\left(-\frac{\|\boldsymbol{x}-\boldsymbol{z}\|^{2}}{2\ell^{2}}\right)$
(16) $\displaystyle k_{OU}(\boldsymbol{x},\boldsymbol{z})$
$\displaystyle\colonequals\theta\exp\left(-\frac{\|\boldsymbol{x}-\boldsymbol{z}\|}{\ell}\right).$
(17)

where we fix $\sigma^{2}\colonequals 10^{-3}$ and $\theta\colonequals 1$, and
we vary $\ell$ as $\log\ell\in\\{-1,0,1,2,3\\}$.

###  A.1 Additional results for hyper-parameter tuning

Table 1 summarizes the results reporting for each data-set the averages over
the last points.

DATASET | ALGORITHM | RMSE | EXACT NEGATIVE MLL
---|---|---|---
protein | CGLB | $0.57\pm 0.00$ | $27620.75\pm 35.57$
StoppedCholesky | $0.60\pm 0.03$ | $28399.84\pm 3299.38$
ExactGPR | $0.56\pm 0.00$ | $23735.92\pm 152.59$
wilson-kin40k | CGLB | $0.09\pm 0.00$ | $-16235.72\pm 66.26$
StoppedCholesky | $0.07\pm 0.00$ | $-20843.01\pm 59.27$
ExactGPR | $0.07\pm 0.00$ | $-20837.51\pm 41.83$
metro | CGLB | $0.37\pm 0.06$ | $23974.65\pm 1467.86$
StoppedCholesky | $0.46\pm 0.31$ | $10402.82\pm 23328.20$
ExactGPR | $0.30\pm 0.13$ | $1747.71\pm 17614.45$
pm25 | CGLB | $0.44\pm 0.04$ | $22605.36\pm 1878.74$
StoppedCholesky | $0.55\pm 0.09$ | $24752.24\pm 4790.72$
ExactGPR | $0.44\pm 0.01$ | $19349.50\pm 31.93$
Table 1: Result overview table. Figure 8: Log-marginal likelihood over time
while optimizing hyper-parameters for the metro dataset. Figure 9: Log-
marginal likelihood over time while optimizing hyper-parameters for the pm25
dataset. Figure 10: Log-marginal likelihood over time while optimizing hyper-
parameters for the protein dataset. Figure 11: Log-marginal likelihood over
time while optimizing hyper-parameters for the kin40k dataset. Figure 12: RMSE
over time while optimizing hyper-parameters for the metro dataset. Figure 13:
RMSE over time while optimizing hyper-parameters for the pm25 dataset. Figure
14: RMSE over time while optimizing hyper-parameters for the protein dataset.
Figure 15: RMSE over time while optimizing hyper-parameters for the kin40k
dataset. Figure 16: RMSE over time while optimizing hyper-parameters for the
bike dataset. Figure 17: RMSE over time while optimizing hyper-parameters for
the elevators dataset. Figure 18: RMSE over time while optimizing hyper-
parameters for the pole dataset. Figure 19: RMSE over time while optimizing
hyper-parameters for the pumadyn32nm dataset.

###  A.2 Bounds for experiments on metro

(a) SE kernel, $\log\ell=-1$

(b) SE kernel, $\log\ell=-1$

(c) SE kernel, $\log\ell=0$

(d) SE kernel, $\log\ell=0$

(e) SE kernel, $\log\ell=1$

(f) SE kernel, $\log\ell=1$

(g) SE kernel, $\log\ell=2$

(h) SE kernel, $\log\ell=2$

Figure 20: Upper and lower bounds on the log-determinant term (left column)
and the quadratic term (right column) for the metro dataset when using a
squared exponential (SE) kernel..

(a) OU kernel, $\log\ell=-1$

(b) OU kernel, $\log\ell=-1$

(c) OU kernel, $\log\ell=0$

(d) OU kernel, $\log\ell=0$

(e) OU kernel, $\log\ell=1$

(f) OU kernel, $\log\ell=1$

(g) OU kernel, $\log\ell=2$

(h) OU kernel, $\log\ell=2$

Figure 21: Upper and lower bounds on the log-determinant term (left column)
and the quadratic term (right column) for the metro dataset using an Ornstein-
Uhlenbeck (OU) kernel.

###  A.3 Bounds for experiments on pm25

(a) SE kernel, $\log\ell=-1$

(b) SE kernel, $\log\ell=-1$

(c) SE kernel, $\log\ell=0$

(d) SE kernel, $\log\ell=0$

(e) SE kernel, $\log\ell=1$

(f) SE kernel, $\log\ell=1$

(g) SE kernel, $\log\ell=2$

(h) SE kernel, $\log\ell=2$

Figure 22: Upper and lower bounds on the log-determinant term (left column)
and the quadratic term (right column) for the pm25 dataset.

(a) OU kernel, $\log\ell=-1$

(b) OU kernel, $\log\ell=-1$

(c) OU kernel, $\log\ell=0$

(d) OU kernel, $\log\ell=0$

(e) OU kernel, $\log\ell=1$

(f) OU kernel, $\log\ell=1$

(g) OU kernel, $\log\ell=2$

(h) OU kernel, $\log\ell=2$

Figure 23: Upper and lower bounds on the log-determinant term (left column)
and the quadratic term (right column) for the pm25 dataset using an Ornstein-
Uhlenbeck (OU) kernel.

###  A.4 Bounds for experiments on protein

(a) SE kernel, $\log\ell=-1$

(b) SE kernel, $\log\ell=-1$

(c) SE kernel, $\log\ell=0$

(d) SE kernel, $\log\ell=0$

(e) SE kernel, $\log\ell=1$

(f) SE kernel, $\log\ell=1$

(g) SE kernel, $\log\ell=2$

(h) SE kernel, $\log\ell=2$

Figure 24: Upper and lower bounds on the log-determinant term (left column)
and the quadratic term (right column) for the protein dataset.

(a) OU kernel, $\log\ell=-1$

(b) OU kernel, $\log\ell=-1$

(c) OU kernel, $\log\ell=0$

(d) OU kernel, $\log\ell=0$

(e) OU kernel, $\log\ell=1$

(f) OU kernel, $\log\ell=1$

(g) OU kernel, $\log\ell=2$

(h) OU kernel, $\log\ell=2$

Figure 25: Upper and lower bounds on the log-determinant term (left column)
and the quadratic term (right column) for the protein dataset using an
Ornstein-Uhlenbeck (OU) kernel.

##  Appendix B Notation

We use a python-inspired index notation, abbreviating for example
$[y_{1},\ldots,y_{n}]^{\top}$ as $\bm{y}_{:n}$—observe that the indexing
starts at 1. inline]check that that’s what we need Indexing binds before any
other operation such that $\boldsymbol{K}_{:s,:s}^{-\\!1}$ is the inverse of
$\boldsymbol{K}_{:s,:s}$ and _not_ all elements up to $s$ of
$\boldsymbol{K}^{-\\!1}$. Differing from the main paper, we assume a
heteroskedastic noise model such that we exchange
$\sigma^{2}\in\mathbb{R}^{+}$ for a function of the inputs,
$\sigma^{2}:\mathbb{X}\rightarrow\mathbb{R}^{+}$. With $\sigma^{2}$ we will
refer to $\inf_{\boldsymbol{x}\in\mathbb{X}}\sigma^{2}(\boldsymbol{x})$, which
we assume to be strictly larger than $0$. For $s\in\\{1,\dots,N\\}$ define
$\mathcal{F}_{s}\colonequals\sigma(\boldsymbol{x}_{1},y_{1},\dots,\boldsymbol{x}_{s},y_{s})$
to be the $\sigma$-algebra generated by
$\boldsymbol{x}_{1},y_{1},\dots,\boldsymbol{x}_{s},y_{s}$. With respect to the
main article, we change the letter $M$ to $t$. The motivation for the former
notation is to highlight the role of the variable as a subset size, whereas in
this part, the focus is on $M$ as a stopping time.

##  Appendix C Proof Sketch

In this section of the appendix, we provide additional details, proofs and
theorems on the proposed formulation. The principal equations included in Sec.
3 of the main manuscript are also included here for a better comprehension.

###  C.1 The cumulative perspective

The key issue this paper is concerned with is how to estimate the full
marginal likelihood, $p(\bm{y})$, given only a subset of $n$ observations and
their combined marginal likelihood, $p(\bm{y}_{:n})$. In particular, we will
derive bounds, which are functions of seen observations, on this estimate.
These bounds will allow us to decide, on the fly, when we have seen enough
observations to accurately estimate the full marginal likelihood.

We can write $\log p(\boldsymbol{y})$ equivalently as

$\displaystyle\log p(\bm{y})=\sum_{n=1}^{N}\log
p(y_{n}\mid\boldsymbol{y}_{:n-1}).$ (18)

With this equation in hand, the phenomena shown in Fig. 1 of the main
manuscript becomes much clearer: The figure shows the value of Equation 18 for
an increasing number of observations $N$. When the plot exhibits a linear
trend it is because the summands $\log p(y_{n}\mid\boldsymbol{y}_{:n-1})$
become approximately constant, implying that the model is not gaining
additional knowledge after the $n$th observation.555An alternative way of
understanding the linear trend is that the spectrum of the covariance matrix
$\boldsymbol{K}$ typically drop to $\sigma^{2}$ at some point; since the log-
determinant is the sum of the log-eigenvalues then the linear trend comes from
additional $2\log\sigma$ terms in the sum. From this perspective, we can craft
an approximation by an _optimal stopping problem_ : after processing
observation $n$, we may decide whether to continue processing the sum or
whether to stop and to estimate the remaining $N\\!-\\!n$ terms.

###  C.2 Extrapolation

For each potential stopping point $t$ we can decompose Equation 18 into a sum
of terms which have already been computed and a remaining sum

$\displaystyle\log p(\boldsymbol{y})$
$\displaystyle=\underbrace{\sum_{n=1}^{t}\log
p(y_{n}\mid\boldsymbol{y}_{:n-1})}_{A:\text{
processed}}+\\!\\!\\!\underbrace{\sum_{n=t+1}^{N}\\!\\!\\!\log
p(y_{n}\mid\boldsymbol{y}_{:n-1})}_{B:\text{ remaining}}.$

It is tempting to estimate $B$ as $\frac{N-t}{t}A$, yet this is estimator is
biased. In the following, we will derive lower and upper bounds,
$\mathcal{L}_{t}$ and $\mathcal{U}_{t}$, such that conditioned on the points
already processed, $B$ can be sandwiched,

$\displaystyle\mathbb{E}[\mathcal{L}_{t}\mid\boldsymbol{x}_{1},y_{1},\dots\boldsymbol{x}_{t},y_{t}]\leq\mathbb{E}[B\mid\boldsymbol{x}_{1},y_{1},\dots\boldsymbol{x}_{t},y_{t}]\leq\mathbb{E}[\mathcal{U}_{t}\negthickspace\mid\boldsymbol{x}_{1},y_{1},\dots\boldsymbol{x}_{t},y_{t}].$
(19)

These bounds tighten as we increase the number of observations, which allow us
to monitor convergence of the approximation. We can then detect when the upper
and lower bounds are sufficiently near each other, and stop computations early
when the approximation is sufficiently good. This is in contrast to other
approximations, where one specifies a computational budget, rather than a
desired accuracy.

###  C.3 General bounds

A more practical-minded reader may safely skip this section and continue in
Section C.6 where we show how to use the bounds to obtain a stopped Cholesky
decomposition. Recall that we want to detect the case when the log-marginal
likelihood starts to behave linearly with more processed datapoints as shown
in Fig. 1 in the main manuscript. That is to say, the bounds presented in the
following are valid in general but useful only in the linear setting.

The posterior of the $n$th observation conditioned on the previous is Gaussian
with

$\displaystyle p(y_{n}\mid\boldsymbol{y}_{:n-1})=$
$\displaystyle\mathcal{N}(m_{n-1}(\boldsymbol{x}_{n}),k_{n-1}(\boldsymbol{x}_{n},\boldsymbol{x}_{n})+\sigma^{2}(\boldsymbol{x}_{n}))$
$\displaystyle m_{n-1}(\boldsymbol{x}_{n})\colonequals$ $\displaystyle
k(\boldsymbol{x}_{n},\boldsymbol{X}_{:n-1})\boldsymbol{K}_{n-1}^{-\\!1}\boldsymbol{y}_{:n-1}$
$\displaystyle k_{n-1}(\boldsymbol{x}_{n},\boldsymbol{x}_{n})\colonequals$
$\displaystyle
k(\boldsymbol{x}_{n},\boldsymbol{x}_{n})-k(\boldsymbol{x}_{n},\boldsymbol{X}_{:n-1})\boldsymbol{K}_{n-1}^{-\\!1}k(\boldsymbol{X}_{:n-1},\boldsymbol{x}_{n}),$

where we assumed (w.l.o.g) that $\mu_{0}(\boldsymbol{x})\colonequals 0$.
Inspecting these expressions one finds that

$\displaystyle\log\det{\boldsymbol{K}_{N}}$
$\displaystyle=\sum_{n=1}^{N}\log\left(k_{n-1}(\boldsymbol{x}_{n},\boldsymbol{x}_{n})+\sigma^{2}(\boldsymbol{x}_{n})\right),$
(20)
$\displaystyle\boldsymbol{y}^{\top}\boldsymbol{K}_{N}^{-\\!1}\boldsymbol{y}$
$\displaystyle=\sum_{n=1}^{N}\frac{\left(y_{n}-m_{n-1}(\boldsymbol{x}_{n})\right)^{2}}{k_{n-1}(\boldsymbol{x}_{n},\boldsymbol{x}_{n})+\sigma^{2}(\boldsymbol{x}_{n})}.$
(21)

These expressions are permutation invariant. This allows us to prove that
these terms cannot be too far from their expected values using a _Hoeffding’s
inequality for super-martingales_ by Fan et al. (2012). This observation also
holds for the conditional case, that is, after having observed
$(\boldsymbol{x}_{1},y_{1}),\dots,(\boldsymbol{x}_{t},y_{t})$. Taking for
granted that we can bound $\mathbb{P}\left(|\log
p(\boldsymbol{y})-\mathbb{E}[\log
p(\boldsymbol{y})\mid\boldsymbol{x}_{1},y_{1},\dots,\boldsymbol{x}_{\tau},y_{\tau}]|>\epsilon\right)$
for stopping times $\tau$, we proceed with the estimation of the expectation.

Recall that $\frac{N-t}{t}\log p(\boldsymbol{y}_{:t})$ is _not_ an unbiased
estimator for $\mathbb{E}[\log
p(\boldsymbol{y})\mid\boldsymbol{x}_{1},y_{1},\dots,\boldsymbol{x}_{t},y_{t}]$,
due to the interaction of
$\boldsymbol{x}_{t+1},y_{t+1},\dots,\boldsymbol{x}_{N},y_{N}$. Our strategy is
to find function families $u$ (and $l$) which upper (and lower) bound the
expectation

$\displaystyle l_{n,t}^{d}\leq_{E}\log
k_{n-1}(\boldsymbol{x}_{n},\boldsymbol{x}_{n})\leq_{E}u_{n,t}^{d}$
$\displaystyle
l_{n,t}^{q}\leq_{E}\frac{\left(y_{n}-m_{n-1}(\boldsymbol{x}_{n})\right)^{2}}{k_{n-1}(\boldsymbol{x}_{n},\boldsymbol{x}_{n})+\sigma^{2}(\boldsymbol{x}_{n})}\leq_{E}u_{n,t}^{q},$

where $\leq_{E}$ denotes that the inequality holds in expectation. We will
choose the function families such that the unseen variables interact only in a
_controlled_ manner. More specifically,

$\displaystyle
f^{x}_{n,t}(\boldsymbol{x}_{n},y_{n},\dots,\boldsymbol{x}_{1},y_{1})=\sum_{j=s+1}^{n}g_{t}^{f,x}(\boldsymbol{z}_{n},\boldsymbol{z}_{j};\boldsymbol{z}_{1},\dots\boldsymbol{z}_{s}),$

with $f\in\\{u,l\\}$ and $x\in\\{d,q\\}$. The effect of this restriction
becomes apparent when taking the expectation. The sum over the bounds becomes
the sum of only two terms: variance and covariance, formally:

$\displaystyle\mathbb{E}\left[\sum_{n=s+1}^{N}f^{x}_{n,t}(\boldsymbol{z}_{n},\dots,\boldsymbol{z}_{1})\mid\sigma(\boldsymbol{z}_{1},\dots,\boldsymbol{z}_{s})\right]$
(22)
$\displaystyle=\left(N-s\right)\mathbb{E}\left[g(\boldsymbol{z}_{s+1},\boldsymbol{z}_{s+1},\boldsymbol{z}_{1}\dots,\boldsymbol{z}_{n})\mid\sigma(\boldsymbol{z}_{1},\dots,\boldsymbol{z}_{s})\right]$
$\displaystyle+\left(N-s\right)\frac{N-s+1}{2}\mathbb{E}\left[g(\boldsymbol{z}_{s+1},\boldsymbol{z}_{s+2},\boldsymbol{z}_{1}\dots,\boldsymbol{z}_{n})\mid\sigma(\boldsymbol{z}_{1},\dots,\boldsymbol{z}_{s})\right].$
(23)

We can estimate this expectation from the observations we obtained between $s$
and $t$.

$\displaystyle\approx\frac{N-t}{t-s}\sum_{n=s+1}^{t}g(\boldsymbol{z}_{n},\boldsymbol{z}_{n},\boldsymbol{z}_{1}\dots,\boldsymbol{z}_{s})$
(24)
$\displaystyle+\frac{2(N-t)}{t-s}\frac{N-s+1}{2}\sum_{i=1}^{\frac{t-s}{2}}g(\boldsymbol{z}_{s+2i},\boldsymbol{z}_{s+2i-1},\boldsymbol{z}_{1}\dots,\boldsymbol{z}_{s}).$

As mentioned before, we will present a way of choosing $t$ in Section C.6.

###  C.4 Bounds on the log-determinant

Since the posterior variance of a Gaussian process can never increase with
more data, the average of the (log) posterior variances is an estimator for an
upper bound on the log-determinant (Bartels et al., 2021). Hence in this case,
we simply ignore the interaction between the remaining variables. We set
$g(\boldsymbol{x}_{n},\boldsymbol{x}_{i})\colonequals\delta_{ni}\log\left(k_{s}(\boldsymbol{x}_{n},\boldsymbol{x}_{n})+\sigma^{2}(\boldsymbol{x}_{n})\right)$
where $\delta_{ni}$ denotes Kronecker’s $\delta$.

To obtain a lower bound we use that for $c>0$ and $a\geq b\geq 0$, one can
show that $\log\left(c+a-b\right)\geq\log\left(c+a\right)-\frac{b}{c}$ where
the smaller $b$ the better the bound. In our case
$c=\sigma^{2}(\boldsymbol{x}_{n})$,
$a=k_{s}(\boldsymbol{x}_{n},\boldsymbol{x}_{n})$ and
$b=k_{s}(\boldsymbol{x}_{n},\boldsymbol{X}_{s+1:n-1})\left(k_{s}(\boldsymbol{X}_{s+1:n-1},\boldsymbol{X}_{s+1:n-1})+\sigma^{2}(\boldsymbol{X}_{s+1:n-1})\right)^{-\\!1}k_{s}(\boldsymbol{X}_{s+1:n-1},\boldsymbol{x}_{n})$.
Underestimating the eigenvalues of
$k_{s}(\boldsymbol{X}_{s+1:n-1},\boldsymbol{X}_{s+1:n-1})$ by 0 we obtain a
lower bound, where each quantity can be estimated. Formally, for any $s\leq
t$,

$\displaystyle\log\left(k_{n-1}(\boldsymbol{x}_{n},\boldsymbol{x}_{n})+\sigma^{2}(\boldsymbol{x}_{n})\right)\geq\left(\log\left(k_{s}(\boldsymbol{x}_{n},\boldsymbol{x}_{n})+\sigma^{2}(\boldsymbol{x}_{n})\right)-\sum_{i=s+1}^{n-1}\frac{k_{s}(\boldsymbol{x}_{n},\boldsymbol{x}_{i})^{2}}{\sigma^{2}(\boldsymbol{x}_{n})\sigma^{2}(\boldsymbol{x}_{i})}\right).$
(25)

This bound can be worse than the deterministic lower bound
$\min_{n^{\prime}}\log\sigma^{2}(\boldsymbol{x}_{n^{\prime}})$. It depends on
how large $n$ is, how large the average correlation is and how small
$\sigma^{2}(\cdot)$ is. We can determine the number of steps $n-s$ that this
bound is better by solving a quadratic equation. Denote with $\mu$ the
estimator for the left addend and with $\rho$ the estimator for the second
addend. The tipping point $\psi$ is the solution of
$(\psi-s)\left(\mu-\frac{\psi-s+1}{2}\rho\right)\leq(\psi-s)\min\log\sigma^{2}(\cdot)$.
One solution is $\psi=s$, the other is

$\displaystyle\psi\colonequals\lfloor
s-1+\frac{2}{\rho}\left(\mu-\min\log\sigma^{2}(\cdot)\right)\rfloor.$ (26)

Hence, for $n>\psi$ we set $u^{d}_{n}\colonequals\min\log\sigma^{2}(\cdot)$.

Observe that, the smaller $k_{s}(j,j+1)^{2}$ the closer the bounds. This term
represents the correlation of datapoints conditioned on the $s$ datapoints
observed before. Thus, our bounds come together, when incoming observations
become independent conditioned on what was already observed. Essentially, that
$k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j+1})^{2}=0$ is the basic assumption
of inducing input approximations (Quiñonero-Candela & Rasmussen, 2005).

###  C.5 Bounds on the quadratic form

For an upper bound on the quadratic form we apply a similar trick:

$\displaystyle\frac{x}{c+a-b}\leq\frac{x(c+b)}{c(c+a)},$ (27)

where $x\geq 0$. Further we assume that in expectation the mean square error
improves with more data. inline]the bound below holds only in expectation!
Formally,

$\displaystyle\frac{\left(y_{j}-m_{j-1}(\boldsymbol{x}_{j})\right)^{2}}{k_{j-1}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j})}\leq_{E}\frac{\left(y_{j}-m_{s}(\boldsymbol{x}_{j})\right)^{2}}{\sigma^{2}(\boldsymbol{x}_{j})\left(k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j})\right)}\left(\sigma^{2}(\boldsymbol{x}_{j})+\sum_{i=s+1}^{j-1}\frac{\left(k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{i})\right)^{2}}{\sigma^{2}(\boldsymbol{x}_{i})}\right)$
(28)

For a lower bound observe that

$\displaystyle\boldsymbol{y}^{\top}\boldsymbol{K}^{-\\!1}\boldsymbol{y}=\boldsymbol{y}_{:s}^{\top}\boldsymbol{K}_{s}^{-\\!1}\boldsymbol{y}_{:s}+\left(\boldsymbol{y}_{s+1:N}-m_{s}(\boldsymbol{X}_{t+1:N})\right)^{\top}\boldsymbol{Q}_{s+1:N}^{-\\!1}\left(\boldsymbol{y}_{s+1:N}-m_{s}(\boldsymbol{X}_{t+1:N})\right)$
(29)

where $\boldsymbol{Q}_{s+1:j}\colonequals
k_{s}(\boldsymbol{X}_{s+1:j},\boldsymbol{X}_{s+1:j})+\sigma^{2}(\boldsymbol{X}_{s+1:j})$
with $j\geq s+1$ for the posterior covariance matrix of
$\boldsymbol{X}_{s+1:j}$ conditioned on $\boldsymbol{X}_{1:s}$. We use a trick
we first encountered in Kim & Teh (2018):
$\boldsymbol{y}^{\top}\boldsymbol{A}^{-\\!1}\boldsymbol{y}\geq
2\boldsymbol{y}^{\top}\boldsymbol{b}-\boldsymbol{b}^{\top}\boldsymbol{A}\boldsymbol{b}$,
for any $\boldsymbol{b}$. Applying this inequality directly would result in a
poor lower bound. For brevity introduce
$\boldsymbol{e}\colonequals\boldsymbol{y}_{s+1:N}-m_{s}(\boldsymbol{X}_{t+1:N})$.
We rewrite the second term as

$\displaystyle\boldsymbol{e}^{\top}\operatorname{diag}{\left(\boldsymbol{e}\right)}\left(\operatorname{diag}{\left(\boldsymbol{e}\right)}\boldsymbol{Q}_{s+1:N}\operatorname{diag}{\left(\boldsymbol{e}\right)}\right)^{-\\!1}\operatorname{diag}{\left(\boldsymbol{e}\right)}\boldsymbol{e}$
(30)

Now applying the inequality with
$\boldsymbol{b}\colonequals\alpha\boldsymbol{1}$, we obtain

$\displaystyle
2\alpha\sum_{n=s+1}^{N}\left(y_{n}-m_{s}(\boldsymbol{x}_{n})\right)^{2}-\alpha^{2}\sum_{n,n^{\prime}=s+1}^{N}\left(y_{n}-m_{s}(\boldsymbol{x}_{n})\right)[\boldsymbol{Q}_{s+1:N}]_{nn^{\prime}}\left(y_{n^{\prime}}-m_{s}(\boldsymbol{x}_{n^{\prime}})\right)$
(31)

which is now in the form of Section C.3. After taking the expectation, the
optimal value of $\alpha$ is

$\displaystyle\frac{\mathbb{E}[\left(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1})\right)^{2}\mid\mathcal{F}_{s}]}{\alpha_{\text{den}}}$
(32)

where

$\alpha_{\text{den}}=\mathbb{E}\left[\left(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1})\right)^{2}k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+1})+\sigma^{2}(\boldsymbol{x}_{s+1})\right]\\\
+\mathbb{E}\left[(N-s-1)\left(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1})\right)\left(y_{s+2}-m_{s}(\boldsymbol{x}_{s+2})\right)\left(k_{s}(s+1,s+2)\right)\mid\mathcal{F}_{s}\right],$
(33)

which is $\mathcal{F}_{s}$-measureable and we can estimate it. inline]Any
volunteers to shorten above equation? :D

Observe that, the smaller the square error
$(y_{j}-m_{s}(\boldsymbol{x}_{j}))^{2}$, the closer the bounds. That is, if
the model fit is good, the quadratic form can be easily identified.

###  C.6 Using the Bounds for Stopping the Cholesky

We will use the same stopping strategy as Mnih et al. (2008); Bartels (2020):
when the difference between bounds becomes sufficiently small and their
absolute value is far away from zero. More precisely, when having
deterministic bounds $\mathcal{L}\leq x\leq\mathcal{U}$ on a number $x$, with

$\displaystyle\frac{\mathcal{U}-\mathcal{L}}{2\min(\left|\mathcal{U}\right|,\left|\mathcal{L}\right|)}\leq
r\text{ and }$ (34)
$\displaystyle\operatorname{sign}{\mathcal{U}}=\operatorname{sign}{L},$ (35)

then the relative error of the estimate $\frac{1}{2}(\mathcal{U}+\mathcal{L})$
is less than $r$, that is
$|\frac{\frac{1}{2}(\mathcal{U}+\mathcal{L})-x}{x}|\leq r$.

######  Remark 4.

In our experiments, we do not use $\frac{1}{2}(\mathcal{U}+\mathcal{L})$ as
estimator. When computing gradients for kernel parameters, the points from $s$
to $\tau$ enter the equation with a different weighting than the first $s$
datapoints for which there is no justification. In exploratory experiments
this was not problematic. Nevertheless, we instead use the biased estimator
$(N-\tau)\frac{1}{\tau}\log p(\boldsymbol{y}_{:\tau})$, since we use our own,
faster auto-differentiation implementation for this term.

######  Remark 5.

Another ingredient where we deviate from the theory before is the estimation
of the average correlation. In theory, the estimator is allowed to take only
every second entry of the off-diagonal. Yet, we could use the same estimator
using entries with an offset of $1$ and still would get a valid estimator.
Hence, we the average of the two should also be a decent estimator. We believe
that it should be possible to prove that the average of the two estimators is
also a decent estimator. Therefore, in practice, take the average over all
indices.

The question remains how to use bounds and stopping strategy to derive an
approximation algorithm. We transform the exact Cholesky decomposition for
that purpose. For brevity denote
$\boldsymbol{L}_{s}\colonequals\operatorname{chol}[k(\boldsymbol{X}_{:s},\boldsymbol{X}_{:s})+\sigma^{2}(\boldsymbol{X}_{:s})]$
and $\boldsymbol{T}_{s}\colonequals
k(\boldsymbol{X}_{s+1:},\boldsymbol{X}){\boldsymbol{L}_{s}}^{-\top}$. For any
$s\in\\{1,\dots N\\}$:

$\displaystyle\boldsymbol{L}_{N}=\begin{bmatrix}\boldsymbol{L}_{s}&\boldsymbol{0}\\\
\boldsymbol{T}&\operatorname{chol}\left[k(\boldsymbol{X}_{s+1:,s+1:})-\boldsymbol{T}\boldsymbol{T}^{\top}\right]\end{bmatrix}$
(36)

One can verify that $\boldsymbol{L}_{N}$ is indeed the Cholesky of
$\boldsymbol{K}_{N}$ by evaluating
$\boldsymbol{L}_{N}\boldsymbol{L}_{N}^{\top}$. Observe that
$k(\boldsymbol{X}_{s+1:,s+1:})-\boldsymbol{T}\boldsymbol{T}^{\top}$ is the
posterior covariance matrix of the $\boldsymbol{y}_{s+1:}$ conditioned on
$\boldsymbol{y}_{:s}$. Hence, in the step before the Cholesky of the posterior
covariance matrix is computed, we can estimate our log-determinant bounds.

Similar reasoning applies for solving the linear equation system. We can write

$\displaystyle\boldsymbol{\alpha}_{N}=\begin{bmatrix}\boldsymbol{\alpha}_{s}\\\
\operatorname{chol}\left[k(\boldsymbol{X}_{s+1:,s+1:})-\boldsymbol{T}\boldsymbol{T}^{\top}\right]^{-\\!1}\left(\boldsymbol{y}_{s+1:}-\boldsymbol{T}_{s}\boldsymbol{\alpha}_{s}\right)\end{bmatrix}$
(37)

Now observe that
$\boldsymbol{T}_{s}\boldsymbol{\alpha}_{s}=m_{s}(\boldsymbol{X}_{s+1:})$.
Hence, before the solving the lower equation system (and before computing the
posterior Cholesky), we can compute our bounds for the quadratic form. There
are different options to implement the Cholesky decomposition. We use a
blocked, row-wise implementation (George et al., 1986). For a practical
implementation see Algorithm 1 and Algorithm 2.

Algorithm 1 blocked and recursive formulation of Cholesky decomposition and
Gaussian elimination, augmented with our stopping conditions.

1procedure ACGP($k(\cdot,\cdot)$, $\mu(\cdot)$, $\sigma^{2}(\cdot)$,
$\boldsymbol{X}$, $\boldsymbol{y}$, $m$, $N_{\max}$)

2 $\boldsymbol{A}\leftarrow\boldsymbol{0}^{N_{\max}\times
N_{\max}},\boldsymbol{\alpha}\leftarrow\boldsymbol{0}^{N_{\max}}$ $\sslash$
allocate memory

3 $\boldsymbol{A}_{1:m,1:m}\leftarrow
k(\boldsymbol{X}_{1:m})+\sigma^{2}(\boldsymbol{X}_{1:m})$ $\sslash$ initialize
kernel matrix

4
$\boldsymbol{\alpha}_{1:m}\leftarrow\boldsymbol{y}_{1:m}-\mu(\boldsymbol{X}_{1:m})$
$\sslash$ evaluate mean function for the same datapoints

5
$\boldsymbol{A}_{1:m,1:m}\leftarrow\operatorname{chol}(\boldsymbol{A}_{1:m,1:m})$
$\sslash$ call to low-level Cholesky

6
$\boldsymbol{\alpha}_{1:m}\leftarrow\boldsymbol{A}_{1:m,1:m}^{-\\!1}\boldsymbol{\alpha}_{1:m}$
$\sslash$ second back-substitution step

7 $i\leftarrow m+1$, $j\leftarrow\min(i+m,N)$

8 while $i<N_{\max}$ do

9 $\boldsymbol{A}_{i:j,1:i}\leftarrow
k(\boldsymbol{X}_{i:j},\boldsymbol{X}_{1:i})$ $\sslash$ evaluate required
block-off-diagonal part of the kernel matrix

10
$\boldsymbol{A}_{i:j,1:i}\leftarrow\boldsymbol{A}_{i:j,1:i}\boldsymbol{A}_{1:i,1:i}^{-\top}$
$\sslash$ solve triangular linear equation system

11 $\boldsymbol{A}_{i:j,i:j}\leftarrow
k(\boldsymbol{X}_{i:j})+\sigma^{2}(\boldsymbol{X}_{i:j})$ $\sslash$ evaluate
required block-diagonal part of the kernel matrix

12
$\boldsymbol{\alpha}_{i:j}\leftarrow\boldsymbol{y}_{i:j}-\mu(\boldsymbol{X}_{i:j})$
$\sslash$ evaluate mean function for the same datapoints

13
$\boldsymbol{A}_{i:j,i:j}\leftarrow\boldsymbol{A}_{i:j,i:j}-\boldsymbol{A}_{i:j,1:i}\boldsymbol{A}_{i:j,1:i}^{\top}$$\sslash$
down-date

14 $\sslash$ now $\boldsymbol{A}_{i:j,i:j}=\boldsymbol{Q}_{s+1:j}$

15
$\boldsymbol{\alpha}_{i:j}\leftarrow\boldsymbol{\alpha}_{i:j}-\boldsymbol{A}_{i:j,1:i}\boldsymbol{\alpha}_{1:i}$$\sslash$
now $\boldsymbol{\alpha}_{i:j}$ contains
$\boldsymbol{y}_{i:j}-m_{i}(\boldsymbol{X}_{i:j})$

16 $\mathcal{L},\mathcal{U}\leftarrow$EvaluateBounds($i$, $j$) $\sslash$ costs
$\mathcal{O}(j-i)$

17 if Equations 34 and 35 fulfilled then

18 return estimator

19 end if

20
$\boldsymbol{A}_{i:j,i:j}\leftarrow\operatorname{chol}(\boldsymbol{A}_{i:j,i:j})$
$\sslash$ finish computing Cholesky for data-points up to index $j$

21
$\boldsymbol{\alpha}_{i:j}\leftarrow\boldsymbol{A}_{i:j,i:j}^{-\\!1}\boldsymbol{\alpha}_{i:j}$
$\sslash$ finish solving linear equation system for index up to $j$

22 $i\leftarrow i+m$, $j\leftarrow\min(i+m,N_{\max})$

23 end while$\sslash$ now $\boldsymbol{A}=\boldsymbol{L}$ and
$\boldsymbol{\alpha}=\boldsymbol{L}^{-\\!1}(\boldsymbol{y}-\mu(\boldsymbol{X}))$

24 return estimator

25end procedure

Algorithm 2 bound algorithm as used in our experiments. The algorithm deviates
slightly from our theory. We use Equation 64 for the upper bound in the
quadratic form, and we use all off-diagonal entries (instead of only every
second). inline]the upper bound on the quadratic is not up-to-date.

1procedure EvaluateBounds($s$, $t$)

2 $D\leftarrow\sum_{j=1}^{s}2\log\bm{A}_{jj}$$\sslash$ in practice we reuse
the sum from the last iteration

3 $Q\leftarrow\sum_{j=1}^{s}\boldsymbol{\alpha}_{j}^{2}$

4 $\mu\leftarrow\frac{1}{t-s}\sum_{j=s+1}^{t}2\log\bm{A}_{jj}$ $\sslash$
average variance of the new points conditioned on all points processed until
$s$

5 $\mathcal{U}_{D}\leftarrow D+(N-s)\mu$

6 $\rho\leftarrow\frac{1}{t-s-1}\sum_{j=s+1}^{t-1}\bm{A}_{j,j+1}^{2}$
$\sslash$ average square correlation (deviating from theory!)

7 $\psi\leftarrow\lfloor
s-1+\frac{2}{\rho}\left(\mu-\log\sigma^{2}\right)\rfloor$$\sslash$ number of
steps the probabilistic bound is better than the deterministic

8 $\mathcal{L}_{D}\leftarrow
D+(\psi-s)\left(\mu-\frac{\psi-s-1}{2}\rho\right)+(N-\psi)\log\sigma^{2}$

9 $\mathcal{U}_{Q}\leftarrow
Q+\frac{N-s}{t-s-1}\sum_{j=s+1}^{t-1}\frac{\boldsymbol{\alpha}^{2}_{j}}{\bm{A}_{j,j}\sigma^{2}(\boldsymbol{x}_{j})}\left(\sigma^{2}(\boldsymbol{x}_{j})+\frac{\bm{A}_{j,j+1}^{2}}{\sigma^{2}(\boldsymbol{x}_{j+1})}\right)$

10 $m\leftarrow\frac{1}{t-s}\sum_{j=s+1}^{t}\boldsymbol{\alpha}_{j}^{2}$
$\sslash$ mean square error

11
$\mu\leftarrow\frac{1}{t-s}\sum_{j=s+1}^{t}\boldsymbol{\alpha}_{j}^{2}\bm{A}_{j,j}$

12
$\rho\leftarrow\frac{1}{t-s-1}\sum_{j=s+1}^{t-1}\boldsymbol{\alpha}_{j}\boldsymbol{\alpha}_{j+1}\bm{A}_{j,j+1}$

13 $\alpha\leftarrow\frac{m}{\mu+(N-s-1)\rho}$

14 $\mathcal{L}_{Q}\leftarrow Q+\alpha(N-s)(2m-\alpha(\mu+\rho(N-s-1)))$

15 return $\mathcal{L}_{D}+\mathcal{L}_{Q},\mathcal{U}_{D}+\mathcal{U}_{Q}$

16end procedure

##  Appendix D Assumptions

######  Assumption 6.

Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space and let
$(\boldsymbol{x}_{j},y_{j})_{j=1}^{N}$ be a sequence of independent and
identically distributed random vectors with
$\boldsymbol{x}:\Omega\rightarrow\mathbb{R}^{D}$ and
$y:\Omega\rightarrow\mathbb{R}$.

######  Assumption 7.

For all $s,i,j,t$ with $s<i\leq j\leq N$ and functions
$f(\boldsymbol{x}_{j},\boldsymbol{x}_{i};\boldsymbol{x}_{1},\dots\boldsymbol{x}_{s})\geq
0$

$\displaystyle\mathbb{E}\left[f(\boldsymbol{x}_{j},\boldsymbol{x}_{i})\left(y_{j}-m_{j-1}(\boldsymbol{x}_{j})\right)^{2}\mid\mathcal{F}_{s}\right]\leq\mathbb{E}\left[f(\boldsymbol{x}_{j},\boldsymbol{x}_{i})\left(y_{j}-m_{s}(\boldsymbol{x}_{j})\right)^{2}\mid\mathcal{F}_{s}\right]$
(38)

where
$f(\boldsymbol{x}_{j},\boldsymbol{x}_{i})\in\left\\{\frac{1}{k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j})},\frac{k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{i})^{2}}{(k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j}))\sigma^{2}(\boldsymbol{x}_{j})\sigma^{2}(\boldsymbol{x}_{i})}\right\\}$.

That is, we assume that in expectation the estimator improves with more data.
Note that, $f$ can not depend on any entries of $\boldsymbol{y}$.

##  Appendix E Main Theorem

inline]adapt notation

######  Theorem 8.

Assume that 6 and 7 hold. For any even $m\in\\{2,4,\dots,N-2\\}$ and any
$s\in\\{1,\dots,N-m\\}$, the bounds defined in Equations 9, 8, 12 and 11 hold
in expectation:

$\displaystyle\mathbb{E}[\mathcal{L}_{D}\mid\mathcal{F}_{s}]\leq\mathbb{E}[\log(\det[\boldsymbol{K}])\mid\mathcal{F}_{s}]\leq\mathbb{E}[\mathcal{U}_{D}\mid\mathcal{F}_{s}]\text{
and}$
$\displaystyle\mathbb{E}[\mathcal{L}_{Q}\mid\mathcal{F}_{s}]\leq\mathbb{E}[\boldsymbol{y}^{\top}\boldsymbol{K}^{-\\!1}\boldsymbol{y}\mid\mathcal{F}_{s}]\leq\mathbb{E}[\mathcal{U}_{Q}\mid\mathcal{F}_{s}]\
.$

###### Proof.

Follows from Theorems 9, 16 and 12, and Theorem 28 in Bartels (2020).

∎

##  Appendix F Proof for the Lower Bound on the Determinant

######  Theorem 9.

Assume that 6 holds, and that $m\in\\{2,4,\dots,N\\}$ is an even number. Set
$t\colonequals s+m$, then, for all $s\in\\{1,\dots,N-m\\}$

$\mathbb{E}\left[\mathcal{L}_{D}\mid\mathcal{F}_{s}\right]\leq\mathbb{E}\left[D_{N}\mid\mathcal{F}_{s}\right].$
$\displaystyle\mathcal{L}_{D}$
$\displaystyle\colonequals\log\det\boldsymbol{K}_{:s,:s}+\left(\psi_{t}-s\right)\left(\log\underline{\mu}_{t}-\frac{\psi_{t}-s-1}{2}\tilde{\rho}_{t}\right)+(N-\psi_{t})\log\sigma^{2}$
(39) $\displaystyle\sslash\text{{\emph{the lower bound}}}$
$\displaystyle\log\underline{\mu}_{t}$
$\displaystyle\colonequals\frac{1}{m}\sum_{j=s+1}^{t}\log\left(\sigma^{2}(\boldsymbol{x}_{j})+k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})\right)$
(40) $\displaystyle\sslash\text{{\emph{(under-)estimate of the posterior
variance conditioned on $s$ points}}}$ $\displaystyle\tilde{\rho}_{t}$
$\displaystyle\colonequals\frac{2}{m}\sum_{j=\frac{s+1}{2}}^{\frac{t-1}{2}}\frac{k_{s}(\boldsymbol{x}_{2j+1},\boldsymbol{x}_{2j})^{2}}{\sigma^{2}(\boldsymbol{x}_{2j+1})\sigma^{2}(\boldsymbol{x}_{2j})}$
(41) $\displaystyle\sslash\text{{\emph{(over-)estimate of the correlation
conditioned on $s$ points}}}$ $\displaystyle\psi_{t}$
$\displaystyle\colonequals s+\max p\text{ where $p$ is such that}$ (42)
$\displaystyle\qquad
p\left(\log\underline{\mu}_{t}-\frac{p-1}{2}\tilde{\rho}_{t}\right)\geq
p\log\sigma^{2}$ (43) $\displaystyle\sslash\text{{\emph{number of steps that
we suspect the decrease in variance to be controllable}}}$
$\displaystyle=\max(N,\lfloor
s-1+\frac{2}{\rho_{D}}\left(\mu_{D}-\log\sigma^{2}\right)\rfloor)$ (44)

###### Proof.

$\displaystyle\mathbb{E}\left[\mathcal{L}_{D}\mid\mathcal{F}_{s}\right]-\mathbb{E}\left[\log\det\boldsymbol{K}\mid\mathcal{F}_{s}\right]=\mathbb{E}\left[\mathcal{L}_{D}-\log\det\boldsymbol{K}\mid\mathcal{F}_{s}\right]$
$\displaystyle=\mathbb{E}\left[\left(\psi_{t}-s\right)\left(\log\underline{\mu}_{t}-\frac{\psi_{t}-s-1}{2}\tilde{\rho}_{t}\right)+(N-\psi_{t})\log\sigma^{2}-\sum_{j=s+1}^{N}\log\left(k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j})\right)\mid\mathcal{F}_{s}\right]$
$\displaystyle\sslash\text{{\emph{using the definition of $\mathcal{L}_{t}$
and slightly simplifying using
\lx@cref{creftype~refnum}{lemma:log_det_as_variance}}}}$
$\displaystyle\leq\mathbb{E}\left[\left(\psi_{t}-s\right)\left(\log\underline{\mu}_{t}-\frac{\psi_{t}-s-1}{2}\tilde{\rho}_{t}\right)-\sum_{j=s+1}^{\psi_{t}}\log\left(k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j})\right)\mid\mathcal{F}_{s}\right]$
$\displaystyle\sslash\text{{\emph{using that
$\log\sigma^{2}\leq\log\left(k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j})\right)$
for all $j$}}}$
$\displaystyle=\left(\psi_{t}-s\right)\left(\mathbb{E}\left[\log\left(\sigma^{2}+k_{s}(\boldsymbol{x}_{t+1},\boldsymbol{x}_{t+1})\mid\mathcal{F}_{s}\right)\right]-\frac{\psi_{t}-s-1}{2}\mathbb{E}\left[\frac{k_{s}(\boldsymbol{x}_{t+1},\boldsymbol{x}_{t+2})^{2}}{\sigma^{2}(\boldsymbol{x}_{t+1})\sigma^{2}(\boldsymbol{x}_{t+2})}\mid\mathcal{F}_{s}\right]\right)$
$\displaystyle\quad-\mathbb{E}\left[\sum_{j=s+1}^{\psi_{t}}\log\left(k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j})\right)\mid\mathcal{F}_{s}\right]$
$\displaystyle\sslash\text{{\emph{using
\lx@cref{creftype~refnum}{assume:exchangeability}}}}$ $\displaystyle\leq 0$
$\displaystyle\sslash\text{{\emph{\lx@cref{creftype~refnum}{lemma:overestimating_expected_logdet}}}}$

∎

######  Lemma 10.

For $c>0$ and $b\geq a\geq 0:$

$\log(c+b-a)\geq\log(c+b)-\frac{a}{c}$

###### Proof.

For $a=0$, the statement is true with equality. We rewrite the inequality as

$\frac{a}{c}\geq\log\left(\frac{c+b}{c+b-a}\right)=\log\left(1+\frac{a}{c+b-a}\right).$

For the case $a=b$, apply the exponential function on both sides, and the
statement follows from $e^{x}\geq x+1$ for all $x$. For $a\in(0,b)$, consider
$f(a)\colonequals\log(c+b-a)+\frac{a}{c}-\log(c+b)$. The first derivative of
this function is $f^{\prime}(a)=-\frac{1}{c+b-a}+\frac{1}{c}$, which is always
positive for $a\in(0,b)$. Since $f(0)=0$, we must have $f(a)\geq 0$ for all
$a\in(a,b)$. ∎

######  Lemma 11.

For all $n\geq t\geq s$:

$\displaystyle\mathbb{E}\left[\sum_{j=t+1}^{n}\log\left(k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j})\right)\mid\mathcal{F}_{s}\right]\geq(n-t)\Bigg{(}$
$\displaystyle\mathbb{E}\left[\log\left(\sigma^{2}(\boldsymbol{x}_{t+1})+k_{s}(\boldsymbol{x}_{t+1},\boldsymbol{x}_{t+1})\right)\mid\mathcal{F}_{s}\right]$
$\displaystyle-\frac{n-t-1}{2\sigma^{4}}\mathbb{E}\left[k_{s}(\boldsymbol{x}_{t+1},\boldsymbol{x}_{t+2})^{2}\mid\mathcal{F}_{s}\right]\Bigg{)}$

###### Proof.

Introduce
$\overline{\boldsymbol{X}}_{j}\colonequals[\boldsymbol{x}_{s+1},\dots,\boldsymbol{x}_{j-1}]$
with the convention
$k_{s}(\boldsymbol{x}_{s+1},\overline{\boldsymbol{X}}_{s+1})\colonequals 0$.

$\displaystyle\mathbb{E}\left[\sum_{j=t+1}^{n}\log
k_{j-1}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j})\mid\mathcal{F}_{s}\right]$
(45)
$\displaystyle=\mathbb{E}\left[\sum_{j=t+1}^{n}\log\left(\sigma^{2}(\boldsymbol{x}_{j})+k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})-k_{s}(\boldsymbol{x}_{j},\overline{\boldsymbol{X}}_{j})\left(k_{s}(\overline{\boldsymbol{X}}_{j},\overline{\boldsymbol{X}}_{j})+\sigma^{2}(\overline{\boldsymbol{X}})\right)^{-\\!1}k_{s}(\overline{\boldsymbol{X}}_{j},\boldsymbol{x}_{j})\right)\mid\mathcal{F}_{s}\right]$
(46)
$\displaystyle\sslash\text{{\emph{\lx@cref{creftype~refnum}{lemma:post_kernel}}}}$
$\displaystyle\geq\mathbb{E}\left[\sum_{j=t+1}^{n}\left(\log\left(\sigma^{2}(\boldsymbol{x}_{j})+k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})\right)-\frac{1}{\sigma^{2}(\boldsymbol{x}_{j})}k_{s}(\boldsymbol{x}_{j},\overline{\boldsymbol{X}}_{j})\left(k_{s}(\overline{\boldsymbol{X}}_{j},\overline{\boldsymbol{X}}_{j})+\sigma^{2}(\overline{\boldsymbol{X}}_{j})\right)^{-\\!1}k_{s}(\overline{\boldsymbol{X}}_{j},\boldsymbol{x}_{j})\right)\mid\mathcal{F}_{s}\right]$
(47)
$\displaystyle\sslash\text{{\emph{\lx@cref{creftype~refnum}{lemma:log_trick}}}}$
$\displaystyle\geq\mathbb{E}\left[\sum_{j=t+1}^{n}\left(\log\left(\sigma^{2}(\boldsymbol{x}_{j})+k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})\right)-\frac{1}{\sigma^{2}(\boldsymbol{x}_{j})}k_{s}(\boldsymbol{x}_{j},\overline{\boldsymbol{X}}_{j})\left(\sigma^{2}(\overline{\boldsymbol{X}}_{j})\right)^{-\\!1}k_{s}(\overline{\boldsymbol{X}}_{j},\boldsymbol{x}_{j})\right)\mid\mathcal{F}_{s}\right]$
(48) $\displaystyle\sslash\text{{\emph{underestimating
$k_{s}(\overline{\boldsymbol{X}}_{j},\overline{\boldsymbol{X}}_{j})$ by
$\boldsymbol{0}$}}}$
$\displaystyle\geq\mathbb{E}\left[\sum_{j=t+1}^{n}\left(\log\left(\sigma^{2}(\boldsymbol{x}_{j})+k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})\right)-\frac{1}{\sigma^{2}(\boldsymbol{x}_{j})}\sum_{i=t+1}^{j-1}\frac{k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})^{2}}{\sigma^{2}(\boldsymbol{x}_{i})}\right)\mid\mathcal{F}_{s}\right]$
(49) $\displaystyle\sslash\text{{\emph{writing the vector multiplication as
sum}}}$
$\displaystyle=(n-t)\mathbb{E}\left[\log\left(\sigma^{2}+k_{s}(\boldsymbol{x}_{t+1},\boldsymbol{x}_{t+1})\right)\mid\mathcal{F}_{s}\right]$
(50)
$\displaystyle\quad+\frac{(n-t)(n-t-1)}{2}\mathbb{E}\left[\frac{k_{s}(\boldsymbol{x}_{t+1},\boldsymbol{x}_{t+2})^{2}}{\sigma^{2}(\boldsymbol{x}_{t+1})\sigma^{2}(\boldsymbol{x}_{t+2})}\mid\mathcal{F}_{s}\right]$
$\displaystyle\sslash\text{{\emph{using
\lx@cref{creftype~refnum}{assume:exchangeability} and then applying
\lx@cref{creftype~refnum}{lemma:little_gauss}}}}$

∎

##  Appendix G Proof for the Upper Bound on the Quadratic Form

######  Theorem 12.

Assume that 6 and 7 hold. Let $m\in\mathbb{N}$ be even, then for all
$s\in\\{1,\dots,N-m\\}$

$\mathbb{E}[\boldsymbol{y}^{\top}\boldsymbol{K}^{-\\!1}\boldsymbol{y}\mid\mathcal{F}_{s}]\leq\mathbb{E}[\mathcal{U}_{Q}\mid\mathcal{F}_{s}]\
,$

where

$\displaystyle\mathcal{U}_{Q}$
$\displaystyle\colonequals\boldsymbol{y}_{:s}^{\top}\boldsymbol{K}_{:s,:s}^{-\\!1}\boldsymbol{y}_{:s}+\left(N-s\right)\left(\mu_{t}+\frac{N-s-1}{2}\tilde{\rho}_{t}\right)$
(51) $\displaystyle\sslash\text{{\emph{the upper bound}}}$
$\displaystyle\mu_{t}$
$\displaystyle\colonequals\frac{1}{t-s}\sum_{j=s+1}^{t}\frac{(y_{j}-m_{s}(\boldsymbol{x}_{j}))^{2}}{k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j})}$
(52) $\displaystyle\tilde{\rho}_{t}$
$\displaystyle\colonequals\frac{2}{t-s}\sum_{j=\frac{s+2}{2}}^{\frac{t}{2}}\frac{(y_{2j}-m_{s}(\boldsymbol{x}_{2j}))^{2}k_{s}(\boldsymbol{x}_{2j},\boldsymbol{x}_{2j-1})^{2}}{\left(k_{s}(\boldsymbol{x}_{2j},\boldsymbol{x}_{2j})+\sigma^{2}(\boldsymbol{x}_{2j})\right)\sigma^{2}(\boldsymbol{x}_{2j})\sigma^{2}(\boldsymbol{x}_{2j-1})}\
.$ (53)

###### Proof.

$\displaystyle\mathbb{E}\left[\boldsymbol{y}^{\top}\boldsymbol{K}^{-\\!1}\boldsymbol{y}\mid\mathcal{F}_{s}\right]-\mathbb{E}\left[\mathcal{U}_{Q}\mid\mathcal{F}_{s}\right]$
$\displaystyle=\mathbb{E}\left[\sum_{j=s+1}^{N}\frac{\left(y_{j}-m_{j-1}(\boldsymbol{x}_{j})\right)^{2}}{k_{j-1}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j})}-\left(N-s\right)\left(\mu_{t}+\frac{N-s-1}{2}\tilde{\rho}_{t}\right)\mid\mathcal{F}_{s}\right]$
$\displaystyle\sslash\text{{\emph{using the definition of $\mathcal{U}_{Q}$
and slightly simplifying with \lx@cref{creftype~refnum}{thm:quad_form_sum}}}}$
$\displaystyle=\mathbb{E}\left[\sum_{j=s+1}^{\psi_{t}}\frac{\left(y_{j}-m_{j-1}(\boldsymbol{x}_{j})\right)^{2}}{k_{j-1}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j})}\mid\mathcal{F}_{s}\right]$
$\displaystyle-(N-s)\left(\mathbb{E}\left[\frac{(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1}))^{2}}{k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+1})+\sigma^{2}(\boldsymbol{x}_{s+1})}\mid\mathcal{F}_{s}\right]\right)$
$\displaystyle-(N-s)\frac{N-s-1}{2}\mathbb{E}\left[\frac{(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1}))^{2}k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+2})^{2}}{\left(k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+1})+\sigma^{2}(\boldsymbol{x}_{s+1})\right)\sigma^{2}(\boldsymbol{x}_{s+2})}\mid\mathcal{F}_{s}\right]$
$\displaystyle\sslash\text{{\emph{using
\lx@cref{creftype~refnum}{assume:exchangeability}}}}$ $\displaystyle\leq 0$
$\displaystyle\sslash\text{{\emph{\lx@cref{creftype~refnum}{lemma:overestimating_expected_qform}}}}$

∎

###  G.1 Tools for this proof

######  Lemma 13.

For $c>0$, $b,x\geq 0$ and $a\geq b$:

$\frac{x}{c+a-b}\leq\frac{x}{c}\left(1-\frac{a-b}{c+a}\right)=\frac{x(c+b)}{c(c+a)}$

###### Proof.

$\displaystyle\frac{x}{c+a-b}$
$\displaystyle=\frac{x}{c}\left(\frac{c}{c+a-b}\right)$
$\displaystyle=\frac{x}{c}\left(1-\frac{a-b}{c+a-b}\right)$
$\displaystyle\leq\frac{x}{c}\left(1-\frac{c+a-b}{c+a}\frac{a-b}{c+a-b}\right)$
$\displaystyle\sslash\text{{\emph{since $\frac{c+a-b}{c+a}\leq 1$}}}$
$\displaystyle=\frac{x}{c}\left(1-\frac{a-b}{c+a}\right)$
$\displaystyle\sslash\text{{\emph{cancelling terms}}}$

∎

######  Lemma 14.

For all $s,t$ and $n$ with $n\geq t\geq s$:

$\displaystyle\mathbb{E}\left[\sum_{j=t+1}^{n}\frac{(y_{j}-m_{j-1}(\boldsymbol{x}_{j}))^{2}}{k_{j-1}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j})}\mid\mathcal{F}_{s}\right]$
$\displaystyle\leq(n-t)\left(\mathbb{E}\left[\frac{(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1}))^{2}}{k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+1})+\sigma^{2}(\boldsymbol{x}_{s+1})}\mid\mathcal{F}_{s}\right]\right)$
$\displaystyle+(n-t)\left(\left(\frac{n+t+1}{2}-s\right)\mathbb{E}\left[\frac{(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1}))^{2}k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+2})^{2}}{\left(k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+1})+\sigma^{2}(\boldsymbol{x}_{s+1})\right)\sigma^{2}(\boldsymbol{x}_{s+2})}\mid\mathcal{F}_{s}\right]\right)$

###### Proof.

Introduce
$\overline{\boldsymbol{X}}_{j}\colonequals[\boldsymbol{x}_{s+1},\dots,\boldsymbol{x}_{j-1}]$
with the convention
$k_{s}(\boldsymbol{x}_{s+1},\overline{\boldsymbol{X}}_{s+1})\colonequals 0$.

$\displaystyle\mathbb{E}\left[\sum_{j=t+1}^{n}\frac{(y_{j}-m_{j-1}(\boldsymbol{x}_{j}))^{2}}{k_{j-1}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j})}\mid\mathcal{F}_{s}\right]$
(54)
$\displaystyle=\mathbb{E}\left[\sum_{j=t+1}^{n}\frac{(y_{j}-m_{j-1}(\boldsymbol{x}_{j}))^{2}}{\sigma^{2}(\boldsymbol{x}_{j})+k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})-k_{s}(\boldsymbol{x}_{j},\overline{\boldsymbol{X}}_{j})\left(k_{s}(\overline{\boldsymbol{X}}_{j},\overline{\boldsymbol{X}}_{j})+\sigma^{2}(\overline{\boldsymbol{X}})\right)^{-\\!1}k_{s}(\overline{\boldsymbol{X}}_{j},\boldsymbol{x}_{j})}\mid\mathcal{F}_{s}\right]$
(55)
$\displaystyle\sslash\text{{\emph{\lx@cref{creftype~refnum}{lemma:post_kernel}}}}$
$\displaystyle\leq\mathbb{E}\left[\sum_{j=t+1}^{n}\frac{(y_{j}-m_{j-1}(\boldsymbol{x}_{j}))^{2}\left(\sigma^{2}(\boldsymbol{x}_{j})+k_{s}(\boldsymbol{x}_{j},\overline{\boldsymbol{X}}_{j})\left(k_{s}(\overline{\boldsymbol{X}}_{j},\overline{\boldsymbol{X}}_{j})+\sigma^{2}(\overline{\boldsymbol{X}})\right)^{-\\!1}k_{s}(\overline{\boldsymbol{X}}_{j},\boldsymbol{x}_{j})\right)}{\sigma^{2}(\boldsymbol{x}_{j})\left(\sigma^{2}(\boldsymbol{x}_{j})+k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})\right)}\mid\mathcal{F}_{s}\right]$
(56)
$\displaystyle\sslash\text{{\emph{\lx@cref{creftype~refnum}{lemma:fraction_trick2}}}}$
$\displaystyle\leq\mathbb{E}\left[\sum_{j=t+1}^{n}\frac{(y_{j}-m_{j-1}(\boldsymbol{x}_{j}))^{2}\left(\sigma^{2}(\boldsymbol{x}_{j})+k_{s}(\boldsymbol{x}_{j},\overline{\boldsymbol{X}}_{j})\left(\sigma^{2}(\overline{\boldsymbol{X}})\right)^{-\\!1}k_{s}(\overline{\boldsymbol{X}}_{j},\boldsymbol{x}_{j})\right)}{\sigma^{2}(\boldsymbol{x}_{j})\left(\sigma^{2}(\boldsymbol{x}_{j})+k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})\right)}\mid\mathcal{F}_{s}\right]$
(57) $\displaystyle\sslash\text{{\emph{underestimating the eigenvalues of
$k_{s}(\overline{\boldsymbol{X}},\overline{\boldsymbol{X}})$ by $0$}}}$
$\displaystyle=\mathbb{E}\left[\sum_{j=t+1}^{n}\frac{(y_{j}-m_{j-1}(\boldsymbol{x}_{j}))^{2}\left(\sigma^{2}(\boldsymbol{x}_{j})+\sum_{i=s+1}^{j-1}\frac{k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{i})^{2}}{\sigma^{2}(\boldsymbol{x}_{i})}\right)}{\sigma^{2}(\boldsymbol{x}_{j})\left(\sigma^{2}(\boldsymbol{x}_{j})+k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})\right)}\mid\mathcal{F}_{s}\right]$
(58) $\displaystyle\sslash\text{{\emph{writing the vector-product explicitly
as a sum}}}$
$\displaystyle=\sum_{j=t+1}^{n}\left(\mathbb{E}\left[\frac{(y_{j}-m_{j-1}(\boldsymbol{x}_{j}))^{2}}{\sigma^{2}(\boldsymbol{x}_{j})+k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})}\mid\mathcal{F}_{s}\right]+\sum_{i=s+1}^{j-1}\mathbb{E}\left[\frac{(y_{j}-m_{j-1}(\boldsymbol{x}_{j}))^{2}k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{i})^{2}}{\left(k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j})\right)\sigma^{2}(\boldsymbol{x}_{i})\sigma^{2}(\boldsymbol{x}_{j})}\mid\mathcal{F}_{s}\right]\right)$
(59) $\displaystyle\sslash\text{{\emph{linearity of expectation}}}$
$\displaystyle=\sum_{j=t+1}^{n}\left(\mathbb{E}\left[\frac{(y_{j}-m_{s}(\boldsymbol{x}_{j}))^{2}}{\sigma^{2}(\boldsymbol{x}_{j})+k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})}\mid\mathcal{F}_{s}\right]+\sum_{i=s+1}^{j-1}\mathbb{E}\left[\frac{(y_{j}-m_{s}(\boldsymbol{x}_{j}))^{2}k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{i})^{2}}{\left(k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j})\right)\sigma^{2}(\boldsymbol{x}_{i})\sigma^{2}(\boldsymbol{x}_{j})}\mid\mathcal{F}_{s}\right]\right)$
(60) $\displaystyle\sslash\text{{\emph{by assumption
\lx@cref{creftype~refnum}{eq:expected_quadratic_form_assumptions}}}}$
$\displaystyle=\sum_{j=t+1}^{n}\left(\mathbb{E}\left[\frac{(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1}))^{2}}{\sigma^{2}(\boldsymbol{x}_{s+1})+k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+1})}\mid\mathcal{F}_{s}\right]+\sum_{i=s+1}^{j-1}\mathbb{E}\left[\frac{(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1}))^{2}k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+2})^{2}}{\left(k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+1})+\sigma^{2}(\boldsymbol{x}_{s+1})\right)\sigma^{2}(\boldsymbol{x}_{s+2})\sigma^{2}(\boldsymbol{x}_{s+1})}\mid\mathcal{F}_{s}\right]\right)$
(61) $\displaystyle\sslash\text{{\emph{using
\lx@cref{creftype~refnum}{assume:exchangeability}}}}$
$\displaystyle=(n-t)\left(\mathbb{E}\left[\frac{(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1}))^{2}}{\sigma^{2}(\boldsymbol{x}_{s+1})+k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+1})}\mid\mathcal{F}_{s}\right]\right)$
(62)
$\displaystyle+(n-t)\left(\left(\frac{n+t+1}{2}-s\right)\mathbb{E}\left[\frac{(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1}))^{2}k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+2})^{2}}{\left(k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+1})+\sigma^{2}(\boldsymbol{x}_{s+1})\right)\sigma^{2}(\boldsymbol{x}_{s+2})\sigma^{2}(\boldsymbol{x}_{s+1})}\mid\mathcal{F}_{s}\right]\right)$
(63) $\displaystyle\sslash\text{{\emph{by
\lx@cref{creftype~refnum}{lemma:little_gauss}}}}$

∎

######  Remark 15.

Similar to the proof of Theorem 9, we can improve the bound by monitoring how
many steps the sum of average correlations is below the average variance. More
precisely, we solve for the largest $\psi\leq N$ such that

$\mu_{t}+\frac{N-\psi-1}{2}\tilde{\rho}_{t}\leq\frac{1}{t-s}\sum_{j=s+1}^{t}\frac{(y_{j}-m_{s}(\boldsymbol{x}_{\boldsymbol{x}_{j}}))^{2}}{\sigma^{2}(\boldsymbol{x}_{j})}\text{,}$

and replace the upper bound by

$\displaystyle\mathcal{U}_{t}\colonequals\boldsymbol{y}_{:s}^{\top}\boldsymbol{K}_{:s,:s}^{-\\!1}\boldsymbol{y}_{:s}+\left(\psi-s\right)\left(\mu_{t}+\frac{\psi-s-1}{2}\tilde{\rho}_{t}\right)+\frac{N-\psi}{t-s}\sum_{j=s+1}^{t}\frac{(y_{j}-m_{s}(\boldsymbol{x}_{\boldsymbol{x}_{j}}))^{2}}{\sigma^{2}(\boldsymbol{x}_{j})}\text{.}$
(64)

##  Appendix H Proof for the Lower Bound on the Quadratic Form

inline] There seems to be a mistake here, pulling out $\alpha$. But it might
be possible to repair this (to a certain extent). I have to remember, I only
need to show
$\mathbb{E}[\mathcal{L}_{Q}]\leq\mathbb{E}[\boldsymbol{y}^{\top}\boldsymbol{K}^{-\\!1}\boldsymbol{y}]$
and NOT
$\mathcal{L}_{Q}\leq\mathbb{E}[\boldsymbol{y}^{\top}\boldsymbol{K}^{-\\!1}\boldsymbol{y}]$.
Idea: for $j\geq t$ $\mathbb{E}[\alpha
f_{j}\mid\mathcal{F}_{s}]=\mathbb{E}[\alpha\mid\mathcal{F}_{s}]\mathbb{E}[f_{j}\mid\mathcal{F}_{s}]=\mathbb{E}[\alpha\mid\mathcal{F}_{s}]\mathbb{E}[f_{s+k}\mid\mathcal{F}_{s}]$
However! Then
$\mathbb{E}[\alpha\mid\mathcal{F}_{s}]\mathbb{E}[f_{s+k}\mid\mathcal{F}_{s}]\neq\mathbb{E}[\alpha
f_{s+k}\mid\mathcal{F}_{s}]$ if $k<t-s$!!! This could be fixed by choosing
$\boldsymbol{b}\colonequals[\boldsymbol{0}_{m},\boldsymbol{1}_{N-t}]$, then
splitting the elements from $s$ to $t$ in half, using one half to estimate
$\alpha$ and the other half for $\mu$ and $\rho$ and such. inline] Looking
more closely, $\mu^{2}$ is actually a very problematic term. A fast and not
terrible solution is to estimate $\alpha$ from the previous block…

######  Theorem 16.

Assume that 6 and 7 hold. Let $m\in\\{2,\dots,N-2\\}$ be an even number less
than $N$. For $s\in\\{1,\dots,N-m\\}$, let $\alpha$ be
$\mathcal{F}_{s}$-measurable, then

$\mathbb{E}[\mathcal{L}_{Q}\mid\mathcal{F}_{s}]\leq\mathbb{E}[\boldsymbol{y}^{\top}\boldsymbol{K}^{-\\!1}\boldsymbol{y}\mid\mathcal{F}_{s}]$

where

$\displaystyle\mathcal{L}_{Q}$
$\displaystyle\colonequals\boldsymbol{y}_{:s}^{\top}\boldsymbol{K}_{:s,:s}^{-\\!1}\boldsymbol{y}_{:s}+\alpha(N-s)\left(2\mu_{t}-\alpha\mu_{t}^{\prime}-\alpha\frac{N-s}{2}\rho_{t}\right)$
(65) $\displaystyle\mu_{t}$
$\displaystyle\colonequals\frac{1}{t-s}\sum_{j=s+1}^{t}(y_{j}-m_{s}(\boldsymbol{x}_{j}))^{2}$
(66) $\displaystyle\mu^{\prime}_{t}$
$\displaystyle\colonequals\frac{1}{t-s}\sum_{j=s+1}^{t}(y_{j}-m_{s}(\boldsymbol{x}_{j}))^{2}(k_{s}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j}))$
(67) $\displaystyle\rho_{t}$
$\displaystyle\colonequals\frac{2}{t-s}\sum_{j=\frac{s+2}{2}}^{\frac{t}{2}}(y_{2j}-m_{s}(\boldsymbol{x}_{2j}))(y_{2j-1}-m_{s}(\boldsymbol{x}_{2j-1}))k_{s}(\boldsymbol{x}_{2j},\boldsymbol{x}_{2j-1})$
(68)

######  Remark 17.

In our implementation we choose
$\alpha\colonequals\frac{\mu_{t}}{\mu^{\prime}_{t}+(N-s)\rho_{t}}$ which
maximizes the lower bound, though violates the assumption that $\alpha$ is
$\mathcal{F}_{s}$-measurable.

###### Proof.

$\displaystyle\mathbb{E}\left[\mathcal{L}_{Q}\mid\mathcal{F}_{s}\right]-\mathbb{E}\left[\boldsymbol{y}^{\top}\boldsymbol{K}^{-\\!1}\boldsymbol{y}\mid\mathcal{F}_{s}\right]$
$\displaystyle=\mathbb{E}\left[\alpha(N-s)\left(2\mu_{t}-\alpha\mu_{t}^{\prime}-\alpha\frac{N-s}{2}\rho_{t}\right)-\sum_{j=s+1}^{N}\frac{\left(y_{j}-m_{j-1}(\boldsymbol{x}_{j})\right)^{2}}{k_{j-1}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j})}\mid\mathcal{F}_{s}\right]$
$\displaystyle\sslash\text{{\emph{using the definition of $\mathcal{L}_{Q}$
and slightly simplifying}}}$
$\displaystyle=2(N-s)\alpha\mathbb{E}\left[(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1}))^{2}\mid\mathcal{F}_{s}\right]$
$\displaystyle\quad-(N-s)\alpha^{2}\mathbb{E}\left[(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1}))^{2}(k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+1})+\sigma^{2}(\boldsymbol{x}_{s+1}))\mid\mathcal{F}_{s}\right]$
$\displaystyle\quad-\frac{(N-s)^{2}}{2}\alpha^{2}\mathbb{E}\left[(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1}))(y_{s+2}-m_{s}(\boldsymbol{x}_{s+2}))k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+2})\mid\mathcal{F}_{s}\right]$
$\displaystyle\quad-\sum_{j=s+1}^{N}\mathbb{E}\left[\frac{\left(y_{j}-m_{j-1}(\boldsymbol{x}_{j})\right)^{2}}{k_{j-1}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j})}\mid\mathcal{F}_{s}\right]$
$\displaystyle\sslash\text{{\emph{using
\lx@cref{creftype~refnum}{assume:exchangeability}}}}$ $\displaystyle\leq 0$
$\displaystyle\sslash\text{{\emph{using
\lx@cref{creftype~refnum}{thm:underestimating_expected_qform}}}}$

∎

######  Lemma 18.

For all $\mathcal{F}_{s}$-measurable $\alpha\in\mathbb{R}:$

$\displaystyle\mathbb{E}\left[\boldsymbol{y}^{\top}\boldsymbol{K}^{-\\!1}\boldsymbol{y}\mid\mathcal{F}_{s}\right]$
$\displaystyle\geq\boldsymbol{y}_{:s}^{\top}\boldsymbol{K}_{:s,:s}^{-\\!1}\boldsymbol{y}_{:s}+2\alpha(N-s)\mathbb{E}[(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1}))^{2}\mid\mathcal{F}_{s}]$
$\displaystyle\quad-\alpha^{2}(N-s)\mathbb{E}[(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1}))^{2}(k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+1})+\sigma^{2}(\boldsymbol{x}_{s+1}))\mid\mathcal{F}_{s}]$
$\displaystyle\quad-\alpha^{2}\frac{(N-s)^{2}}{2}\mathbb{E}[(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1}))(y_{s+2}-m_{s}(\boldsymbol{x}_{s+2}))k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+2})\mid\mathcal{F}_{s}]$
(69)

###### Proof.

Using Lemma 21, we can write the quadratic form as a sum of two quadratic
forms:

$\displaystyle\boldsymbol{y}^{\top}\boldsymbol{K}^{-\\!1}\boldsymbol{y}$
$\displaystyle=\boldsymbol{y}_{:s}^{\top}\boldsymbol{K}_{:s,:s}^{-\\!1}\boldsymbol{y}_{:s}+(\boldsymbol{y}_{s+1:}-m_{s}(\boldsymbol{X}_{s+1:}))^{\top}\left(k_{s}(\boldsymbol{X}_{s+1:},\boldsymbol{X}_{s+1:})+\sigma^{2}(\boldsymbol{X}_{s+1:})\right)^{-\\!1}(\boldsymbol{y}_{s+1:}-m_{s}(\boldsymbol{X}_{s+1:})).$
(70)

Define
$\boldsymbol{e}\colonequals(\boldsymbol{y}_{s+1:}-m_{s}(\boldsymbol{X}_{s+1:}))$
to rewrite the second addend as

$\displaystyle\boldsymbol{e}^{\top}\operatorname{Diag}(\boldsymbol{e})\left(\operatorname{Diag}(\boldsymbol{e})\left(k_{s}(\boldsymbol{X}_{s+1:},\boldsymbol{X}_{s+1:})+\sigma^{2}(\boldsymbol{X}_{s+1:})\right)\operatorname{Diag}(\boldsymbol{e})\right)^{-\\!1}\operatorname{Diag}(\boldsymbol{e})\boldsymbol{e}$
(71)

We use a trick we first encountered in Kim & Teh (2018):
$\boldsymbol{y}^{\top}\boldsymbol{A}^{-\\!1}\boldsymbol{y}\geq
2\boldsymbol{y}^{\top}\boldsymbol{b}-\boldsymbol{b}^{\top}\boldsymbol{A}\boldsymbol{b}$,
for any $\boldsymbol{b}$. Choose
$\boldsymbol{b}\colonequals\alpha\boldsymbol{1}$ and observe that
$\boldsymbol{b}$ is $\mathcal{F}_{s}$-measurable.

$\displaystyle\mathbb{E}\left[\boldsymbol{y}^{\top}\boldsymbol{K}^{-\\!1}\boldsymbol{y}\mid\mathcal{F}_{s}\right]$
(72)
$\displaystyle=\boldsymbol{y}_{:s}^{\top}\boldsymbol{K}_{:s,:s}^{-\\!1}\boldsymbol{y}_{:s}+\mathbb{E}\left[(\boldsymbol{y}_{s+1:}-m_{s}(\boldsymbol{X}_{s+1:}))^{\top}\left(k_{s}(\boldsymbol{X}_{s+1:},\boldsymbol{X}_{s+1:})+\sigma^{2}(\boldsymbol{X}_{s+1:})\right)^{-\\!1}(\boldsymbol{y}_{s+1:}-m_{s}(\boldsymbol{X}_{s+1:}))\mid\mathcal{F}_{s}\right]$
(73) $\displaystyle\sslash\text{{\emph{since
$\boldsymbol{y}_{:s}^{\top}\boldsymbol{K}_{:s,:s}^{-\\!1}\boldsymbol{y}_{:s}$
is $\mathcal{F}_{s}$-measurable}}}$
$\displaystyle\geq\boldsymbol{y}_{:s}^{\top}\boldsymbol{K}_{:s,:s}^{-\\!1}\boldsymbol{y}_{:s}+2\alpha\boldsymbol{1}^{\top}\mathbb{E}\left[\operatorname{Diag}(\boldsymbol{e})\boldsymbol{e}\right]-\alpha^{2}\boldsymbol{1}^{\top}\mathbb{E}\left[\operatorname{Diag}(\boldsymbol{e})\left(k_{s}(\boldsymbol{X}_{s+1:},\boldsymbol{X}_{s+1:})+\sigma^{2}(\boldsymbol{X}_{s+1:})\right)\operatorname{Diag}(\boldsymbol{e})\mid\mathcal{F}_{s}\right]\boldsymbol{1}$
(74) $\displaystyle\sslash\text{{\emph{applying the inequality for quadratic
forms and using the $\mathcal{F}_{s}$-measurabilituy of $\alpha$}}}$
$\displaystyle=\boldsymbol{y}_{:s}^{\top}\boldsymbol{K}_{:s,:s}^{-\\!1}\boldsymbol{y}_{:s}+2\alpha(N-s)\mathbb{E}[(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1}))^{2}\mid\mathcal{F}_{s}]$
$\displaystyle\quad-\alpha^{2}(N-s)\mathbb{E}[(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1}))^{2}(k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+1})+\sigma^{2}(\boldsymbol{x}_{s+1}))\mid\mathcal{F}_{s}]$
$\displaystyle\quad-\alpha^{2}\left((N-s)\frac{N-s+1}{2}-(N-s)\right)\mathbb{E}[(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1}))(y_{s+2}-m_{s}(\boldsymbol{x}_{s+2}))k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+2})\mid\mathcal{F}_{s}]$
(75) $\displaystyle\sslash\text{{\emph{using
\lx@cref{creftype~refnum}{assume:exchangeability}, grouping variance and
covariance terms separately}}}$
$\displaystyle=\boldsymbol{y}_{:s}^{\top}\boldsymbol{K}_{:s,:s}^{-\\!1}\boldsymbol{y}_{:s}+2\alpha(N-s)\mathbb{E}[(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1}))^{2}\mid\mathcal{F}_{s}]$
$\displaystyle\quad-\alpha^{2}(N-s)\mathbb{E}[(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1}))^{2}(k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+1})+\sigma^{2}(\boldsymbol{x}_{s+1}))\mid\mathcal{F}_{s}]$
$\displaystyle\quad-\alpha^{2}(N-s)\frac{N-s-1}{2}\mathbb{E}[(y_{s+1}-m_{s}(\boldsymbol{x}_{s+1}))(y_{s+2}-m_{s}(\boldsymbol{x}_{s+2}))k_{s}(\boldsymbol{x}_{s+1},\boldsymbol{x}_{s+2})\mid\mathcal{F}_{s}]$
(76) $\displaystyle\sslash\text{{\emph{simplifying}}}$

∎

##  Appendix I Utility Proofs

######  Lemma 19 (Bounding the relative error).

Let $D,\hat{D}\in[\mathcal{L},\mathcal{U}]$, and assume
$\operatorname{sign}(\mathcal{L})=\operatorname{sign}(\mathcal{U})\neq 0$.
Then the relative error of the estimator $\hat{D}$ can be bounded as

$\frac{|D-\hat{D}|}{|D|}\leq\frac{\max(\mathcal{U}-\hat{D},\hat{D}-\mathcal{L})}{\min(|\mathcal{L}|,|\mathcal{U}|)}\,.$

###### Proof.

First observe that if $D_{N}>\hat{D}$ then
$|D_{N}-\hat{D}|=D_{N}-\hat{D}\leq\mathcal{U}-\hat{D}$. If $D_{N}\leq\hat{D}$,
then $|D_{N}-\hat{D}|=\hat{D}-D_{N}\leq\hat{D}-\mathcal{L}$. Hence,

$|D_{N}-\hat{D}|\leq\max(\mathcal{U}-\hat{D},\hat{D}-\mathcal{L}).$

Case $\mathcal{L}>0$: In this case
$|D_{N}|=D_{N}\geq\mathcal{L}=|\mathcal{L}|$, and we obtain for the relative
error:

$\displaystyle\frac{\max(\mathcal{U}-\hat{D},\hat{D}-\mathcal{L})}{|D_{N}|}$
$\displaystyle\leq\frac{\max(\mathcal{U}-\hat{D},\hat{D}-\mathcal{L})}{|\mathcal{L}|}\,.$

Case $\mathcal{U}<0$: In that case
$|\mathcal{L}|\geq|D_{N}|\geq|\mathcal{U}|$, and the relative error can be
bounded as follows.

$\displaystyle\frac{\max(\mathcal{U}-\hat{D},\hat{D}-\mathcal{L})}{|D_{N}|}$
$\displaystyle\leq\frac{\max(\mathcal{U}-\hat{D},\hat{D}-\mathcal{L})}{|\mathcal{U}|}$

Since we assumed
$\operatorname{sign}(\mathcal{L})=\operatorname{sign}(\mathcal{U})$ these were
all cases that required consideration. Combining all observations yields

$\displaystyle\frac{|D_{N}-\hat{D}|}{|D_{N}|}$
$\displaystyle\leq\max(\mathcal{U}-\hat{D},\hat{D}-\mathcal{L})\max\left(\frac{1}{|\mathcal{U}|},\frac{1}{|\mathcal{L}|}\right)$
$\displaystyle=\frac{\max(\mathcal{U}-\hat{D},\hat{D}-\mathcal{L})}{\min(|\mathcal{U}|,|\mathcal{L}|)}$

∎

inline]introduce Corollary for $\frac{1}{2}\mathcal{L}+\frac{1}{2}\mathcal{U}$

######  Lemma 20.

The $\log$-determinant of a kernel matrix can be written as a sum of
conditional variances.

$\displaystyle\log\det\boldsymbol{K}=\sum_{j=1}^{N}\log(k_{j-1}(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}(\boldsymbol{x}_{j}))$
(77)

###### Proof.

Denote with $\boldsymbol{L}$ the Cholesky decomposition of $\boldsymbol{K}$.
Then we obtain

$\displaystyle\log\det\boldsymbol{K}$
$\displaystyle=\log\det(\boldsymbol{L}\boldsymbol{L}^{\top})$ (78)
$\displaystyle\sslash\text{{\emph{using
$\boldsymbol{A}=\boldsymbol{L}\boldsymbol{L}^{\top}$}}}$
$\displaystyle=\log\left(\det(\boldsymbol{L})\det(\boldsymbol{L}^{\top})\right)$
(79) $\displaystyle\sslash\text{{\emph{for square matrices
$\boldsymbol{B},\boldsymbol{C}$:
$\det(\boldsymbol{B}\boldsymbol{C})=\det(\boldsymbol{B})\det(\boldsymbol{C})$}}}$
$\displaystyle=\log\left(\prod_{j=1}^{N}\boldsymbol{L}_{jj}^{2}\right)$ (80)
$\displaystyle\sslash\text{{\emph{for triangular matrices the determinant is
the product of the diagonal elements}}}$
$\displaystyle=\sum_{j=1}^{N}2\log\boldsymbol{L}_{jj}$ (81)
$\displaystyle\sslash\text{{\emph{property of $\log$}}}$

With Lemma 24 the result follows. ∎

######  Lemma 21.

The term $\boldsymbol{y}^{\top}\boldsymbol{K}^{-\\!1}\boldsymbol{y}$ can be
written as

$\displaystyle\boldsymbol{y}^{\top}\boldsymbol{K}^{-\\!1}\boldsymbol{y}=\sum_{n=1}^{N}\frac{\left(y_{n}-m_{n-1}(\boldsymbol{x}_{n})\right)}{k_{n-1}(\boldsymbol{x}_{n},\boldsymbol{x}_{n})+\sigma^{2}(\boldsymbol{x}_{n})}.$
(82)

###### Proof.

Define

$\displaystyle\boldsymbol{k}_{j}(\boldsymbol{x})$
$\displaystyle\colonequals[k(\boldsymbol{x},\boldsymbol{x}_{1}),...,k(\boldsymbol{x},\boldsymbol{x}_{j})]^{\top}\in\mathbb{R}^{j}$
$\displaystyle\boldsymbol{k}_{j+1}$
$\displaystyle\colonequals\boldsymbol{k}_{j}(\boldsymbol{x}_{j+1})\in\mathbb{R}^{j}$
$\displaystyle p_{j}$ $\displaystyle\colonequals
k(\boldsymbol{x}_{j},\boldsymbol{x}_{j})+\sigma^{2}-\boldsymbol{k}_{j}^{\top}(\boldsymbol{K}_{j-1}+\sigma^{2}\boldsymbol{I})^{-\\!1}\boldsymbol{k}_{j}$
$\displaystyle\boldsymbol{\alpha}$
$\displaystyle\colonequals(\boldsymbol{K}_{n}+\sigma^{2}\boldsymbol{I})^{-1}\boldsymbol{k}_{n+1}$

First note, that using block-matrix inversion we can write

$\displaystyle(\boldsymbol{K}_{n+1}+\sigma^{2}\boldsymbol{I})^{-1}=\begin{bmatrix}(\boldsymbol{K}_{n}+\sigma^{2}\boldsymbol{I})^{-1}+\boldsymbol{\alpha}p_{n+1}^{-1}\boldsymbol{\alpha}^{\top}&-\boldsymbol{\alpha}p_{n+1}^{-1}\\\
-\boldsymbol{\alpha}^{\top}p_{n+1}^{-1}&p_{n+1}^{-1}\end{bmatrix}.$

This allows to write

$\displaystyle\boldsymbol{y}_{n+1}^{\top}(\boldsymbol{K}_{n+1}+\sigma^{2}\boldsymbol{I})^{-\\!1}\boldsymbol{y}_{n+1}$
$\displaystyle=\begin{bmatrix}\boldsymbol{y}_{n}^{\top}&y_{n+1}\end{bmatrix}\begin{bmatrix}(\boldsymbol{K}_{n}+\sigma^{2}\boldsymbol{I})^{-1}+\boldsymbol{\alpha}p_{n+1}^{-1}\boldsymbol{\alpha}^{\top}&-\boldsymbol{\alpha}p_{n+1}^{-1}\\\
-\boldsymbol{\alpha}^{\top}p_{n+1}^{-1}&p_{n+1}^{-1}\end{bmatrix}\begin{bmatrix}\boldsymbol{y}_{n}\\\
y_{n+1}\end{bmatrix}$ $\displaystyle\sslash\text{{\emph{using above
observation}}}$
$\displaystyle=\begin{bmatrix}\boldsymbol{y}_{n}^{\top}&y_{n+1}\end{bmatrix}\begin{bmatrix}(\boldsymbol{K}_{n}+\sigma^{2}\boldsymbol{I})^{-1}\boldsymbol{y}_{n}+\boldsymbol{\alpha}p_{n+1}^{-1}\boldsymbol{\alpha}^{\top}\boldsymbol{y}_{n}-\boldsymbol{\alpha}p_{n+1}^{-1}y_{n+1}\\\
-\boldsymbol{\alpha}^{\top}p_{n+1}^{-1}\boldsymbol{y}_{n}+p_{n+1}^{-1}y_{n+1}\end{bmatrix}$
$\displaystyle\sslash\text{{\emph{simplifying from the right}}}$
$\displaystyle=\boldsymbol{y}_{n}^{\top}(\boldsymbol{K}_{n}+\sigma^{2}\boldsymbol{I})^{-1}\boldsymbol{y}_{n}+\boldsymbol{y}_{n}^{\top}\boldsymbol{\alpha}p_{n+1}^{-1}\boldsymbol{\alpha}^{\top}\boldsymbol{y}_{n}-\boldsymbol{y}_{n}^{\top}\boldsymbol{\alpha}p_{n+1}^{-1}y_{n+1}-y_{n+1}\boldsymbol{\alpha}^{\top}p_{n+1}^{-1}\boldsymbol{y}_{n}+y_{n+1}p_{n+1}^{-1}y_{n+1}$
$\displaystyle\sslash\text{{\emph{simplifying from the left}}}$
$\displaystyle=\boldsymbol{y}_{n}^{\top}(\boldsymbol{K}_{n}+\sigma^{2}\boldsymbol{I})^{-1}\boldsymbol{y}_{n}+p_{n+1}^{-1}(\boldsymbol{y}_{n}^{\top}\boldsymbol{\alpha}\boldsymbol{\alpha}^{\top}\boldsymbol{y}_{n}-\boldsymbol{y}_{n}^{\top}\boldsymbol{\alpha}y_{n+1}-y_{n+1}\boldsymbol{\alpha}^{\top}\boldsymbol{y}_{n}+y_{n+1}y_{n+1})$
$\displaystyle\sslash\text{{\emph{pulling out $p_{n+1}^{-1}$}}}$
$\displaystyle=\boldsymbol{y}_{n}^{\top}(\boldsymbol{K}_{n}+\sigma^{2}\boldsymbol{I})^{-1}\boldsymbol{y}_{n}+p_{n+1}^{-1}((\boldsymbol{y}_{n}^{\top}\boldsymbol{\alpha})^{2}-2\boldsymbol{y}_{n}^{\top}\boldsymbol{\alpha}y_{n+1}+y_{n+1}^{2})$
$\displaystyle\sslash\text{{\emph{simplifying}}}$
$\displaystyle=\boldsymbol{y}_{n}^{\top}(\boldsymbol{K}_{n}+\sigma^{2}\boldsymbol{I})^{-1}\boldsymbol{y}_{n}+p_{n+1}^{-1}(\boldsymbol{y}_{n}^{\top}\boldsymbol{\alpha}-y_{n+1})^{2}$
$\displaystyle\sslash\text{{\emph{simplifying}}}$

Now observe that the last addend is indeed the mean square error divided by
the posterior variance. By induction the result follows. ∎

######  Lemma 22.

For all $t,m\in\mathbb{N}$ with $1\leq t+m\leq N$

$\displaystyle
k_{t+m}(\boldsymbol{x}_{a},\boldsymbol{x}_{b})=k_{t}(\boldsymbol{x}_{a},\boldsymbol{x}_{b})-k_{t}(\boldsymbol{x}_{a},\overline{\boldsymbol{X}})\left(k_{t}(\overline{\boldsymbol{X}})+\sigma^{2}\boldsymbol{I}_{m}\right)^{-\\!1}k_{t}(\overline{\boldsymbol{X}},\boldsymbol{x}_{b})$

where $k_{t}(\boldsymbol{x}_{a},\boldsymbol{x}_{b})\colonequals
k(\boldsymbol{x}_{a},\boldsymbol{x}_{b})-k(\boldsymbol{x}_{a},\boldsymbol{X}_{t})\left(k(\boldsymbol{X}_{:t},\boldsymbol{X}_{:t})+\sigma^{2}(\boldsymbol{X}_{:t})\right)^{-\\!1}k(\boldsymbol{X}_{t},\boldsymbol{x}_{b})$
and $\overline{\boldsymbol{X}}\colonequals\boldsymbol{X}_{t:t+m}$.

###### Proof.

$\displaystyle k_{t+m}(\boldsymbol{x}_{a},\boldsymbol{x}_{b})$
$\displaystyle=k(\boldsymbol{x}_{a},\boldsymbol{x}_{b})-k(\boldsymbol{x}_{a},\boldsymbol{X}_{t+m})\boldsymbol{A}_{t+m}^{-\\!1}k(\boldsymbol{X}_{t+m},\boldsymbol{x}_{b})$
$\displaystyle\sslash\text{{\emph{by definition}}}$
$\displaystyle=k(\boldsymbol{x}_{a},\boldsymbol{x}_{b})-\begin{bmatrix}k(\boldsymbol{x}_{a},\boldsymbol{X}_{t})&k(\boldsymbol{x}_{a},\overline{\boldsymbol{X}})\end{bmatrix}\begin{bmatrix}k(\boldsymbol{X}_{t})+\sigma^{2}\boldsymbol{I}_{t}&k(\boldsymbol{X}_{t},\overline{\boldsymbol{X}})\\\
k(\overline{\boldsymbol{X}},\boldsymbol{X}_{t})&k(\overline{\boldsymbol{X}})+\sigma^{2}\boldsymbol{I}_{t}\end{bmatrix}^{-\\!1}\begin{bmatrix}k(\boldsymbol{X}_{t},\boldsymbol{x}_{b})\\\
k(\overline{\boldsymbol{X}},\boldsymbol{x}_{b})\end{bmatrix}$
$\displaystyle\sslash\text{{\emph{in block notation}}}$
$\displaystyle=k(\boldsymbol{x}_{a},\boldsymbol{x}_{b})-\begin{bmatrix}k(\boldsymbol{x}_{a},\boldsymbol{X}_{t})&k(\boldsymbol{x}_{a},\overline{\boldsymbol{X}})\end{bmatrix}\begin{bmatrix}\boldsymbol{A}_{t}&k(\boldsymbol{X}_{t},\overline{\boldsymbol{X}})\\\
k(\overline{\boldsymbol{X}},\boldsymbol{X}_{t})&k(\overline{\boldsymbol{X}})+\sigma^{2}\boldsymbol{I}_{t}\end{bmatrix}^{-\\!1}\begin{bmatrix}k(\boldsymbol{X}_{t},\boldsymbol{x}_{b})\\\
k(\overline{\boldsymbol{X}},\boldsymbol{x}_{b})\end{bmatrix}$
$\displaystyle\sslash\text{{\emph{using the definition of
$\boldsymbol{A}_{t}$}}}$
$\displaystyle=k(\boldsymbol{x}_{a},\boldsymbol{x}_{b})-\begin{bmatrix}k(\boldsymbol{x}_{a},\boldsymbol{X}_{t})&k(\boldsymbol{x}_{a},\overline{\boldsymbol{X}})\end{bmatrix}\cdot$
$\displaystyle{\begin{bmatrix}\boldsymbol{A}_{t}^{-\\!1}+\boldsymbol{A}_{t}^{-\\!1}k(\boldsymbol{X}_{t},\overline{\boldsymbol{X}})\left(k(\overline{\boldsymbol{X}})+\sigma^{2}\boldsymbol{I}_{t}-k(\overline{\boldsymbol{X}},\boldsymbol{X}_{t})\boldsymbol{A}_{t}^{-\\!1}k(\boldsymbol{X}_{t},\overline{\boldsymbol{X}})\right)^{-\\!1}k(\overline{\boldsymbol{X}},\boldsymbol{X}_{t})\boldsymbol{A}_{t}^{-\\!1}&-\boldsymbol{A}_{t}^{-\\!1}k(\boldsymbol{X}_{t},\overline{\boldsymbol{X}})\left(k(\overline{\boldsymbol{X}})+\sigma^{2}\boldsymbol{I}_{t}-k(\overline{\boldsymbol{X}},\boldsymbol{X}_{t})\boldsymbol{A}_{t}^{-\\!1}k(\boldsymbol{X}_{t},\overline{\boldsymbol{X}})\right)^{-\\!1}\\\
-\left(k(\overline{\boldsymbol{X}})+\sigma^{2}\boldsymbol{I}_{t}-k(\overline{\boldsymbol{X}},\boldsymbol{X}_{t})\boldsymbol{A}_{t}^{-\\!1}k(\boldsymbol{X}_{t},\overline{\boldsymbol{X}})\right)^{-\\!1}k(\overline{\boldsymbol{X}},\boldsymbol{X}_{t})\boldsymbol{A}_{t}^{-\\!1}&\left(k(\overline{\boldsymbol{X}})+\sigma^{2}\boldsymbol{I}_{t}-k(\overline{\boldsymbol{X}},\boldsymbol{X}_{t})\boldsymbol{A}_{t}^{-\\!1}k(\boldsymbol{X}_{t},\overline{\boldsymbol{X}})\right)^{-\\!1}\end{bmatrix}\cdot}$
$\displaystyle\begin{bmatrix}k(\boldsymbol{X}_{t},\boldsymbol{x}_{b})\\\
k(\overline{\boldsymbol{X}},\boldsymbol{x}_{b})\end{bmatrix}$
$\displaystyle\sslash\text{{\emph{applying block-matrix inversion}}}$
$\displaystyle=k(\boldsymbol{x}_{a},\boldsymbol{x}_{b})-\begin{bmatrix}k(\boldsymbol{x}_{a},\boldsymbol{X}_{t})&k(\boldsymbol{x}_{a},\overline{\boldsymbol{X}})\end{bmatrix}\cdot$
$\displaystyle\begin{bmatrix}\boldsymbol{A}_{t}^{-\\!1}+\boldsymbol{A}_{t}^{-\\!1}k(\boldsymbol{X}_{t},\overline{\boldsymbol{X}})\left(k_{t}(\overline{\boldsymbol{X}})+\sigma^{2}\boldsymbol{I}_{t}\right)^{-\\!1}k(\overline{\boldsymbol{X}},\boldsymbol{X}_{t})\boldsymbol{A}_{t}^{-\\!1}&-\boldsymbol{A}_{t}^{-\\!1}k(\boldsymbol{X}_{t},\overline{\boldsymbol{X}})\left(k_{t}(\overline{\boldsymbol{X}})+\sigma^{2}\boldsymbol{I}_{t}\right)^{-\\!1}\\\
-\left(k_{t}(\overline{\boldsymbol{X}})+\sigma^{2}\boldsymbol{I}_{t}\right)^{-\\!1}k(\overline{\boldsymbol{X}},\boldsymbol{X}_{t})\boldsymbol{A}_{t}^{-\\!1}&\left(k_{t}(\overline{\boldsymbol{X}})+\sigma^{2}\boldsymbol{I}_{t}\right)^{-\\!1}\end{bmatrix}\cdot$
$\displaystyle\begin{bmatrix}k(\boldsymbol{X}_{t},\boldsymbol{x}_{b})\\\
k(\overline{\boldsymbol{X}},\boldsymbol{x}_{b})\end{bmatrix}$
$\displaystyle\sslash\text{{\emph{applying the definition of $k_{t}$}}}$
$\displaystyle=k(\boldsymbol{x}_{a},\boldsymbol{x}_{b})-\begin{bmatrix}k(\boldsymbol{x}_{a},\boldsymbol{X}_{t})&k(\boldsymbol{x}_{a},\overline{\boldsymbol{X}})\end{bmatrix}\cdot$
$\displaystyle\begin{bmatrix}\boldsymbol{A}_{t}^{-\\!1}k(\boldsymbol{X}_{t},\boldsymbol{x}_{b})+\boldsymbol{A}_{t}^{-\\!1}k(\boldsymbol{X}_{t},\overline{\boldsymbol{X}})\left(k_{t}(\overline{\boldsymbol{X}})+\sigma^{2}\boldsymbol{I}_{t}\right)^{-\\!1}k(\overline{\boldsymbol{X}},\boldsymbol{X}_{t})\boldsymbol{A}_{t}^{-\\!1}k(\boldsymbol{X}_{t},\boldsymbol{x}_{b})-\boldsymbol{A}_{t}^{-\\!1}k(\boldsymbol{X}_{t},\overline{\boldsymbol{X}})\left(k_{t}(\overline{\boldsymbol{X}})+\sigma^{2}\boldsymbol{I}_{t}\right)^{-\\!1}k(\overline{\boldsymbol{X}},\boldsymbol{x}_{b})\\\
-\left(k_{t}(\overline{\boldsymbol{X}})+\sigma^{2}\boldsymbol{I}_{t}\right)^{-\\!1}k(\overline{\boldsymbol{X}},\boldsymbol{X}_{t})\boldsymbol{A}_{t}^{-\\!1}k(\boldsymbol{X}_{t},\boldsymbol{x}_{b})+\left(k_{t}(\overline{\boldsymbol{X}})+\sigma^{2}\boldsymbol{I}_{t}\right)^{-\\!1}k(\overline{\boldsymbol{X}},\boldsymbol{x}_{b})\end{bmatrix}$
$\displaystyle\sslash\text{{\emph{evaluating multiplication with right-most
vector}}}$
$\displaystyle=k(\boldsymbol{x}_{a},\boldsymbol{x}_{b})-\begin{bmatrix}k(\boldsymbol{x}_{a},\boldsymbol{X}_{t})&k(\boldsymbol{x}_{a},\overline{\boldsymbol{X}})\end{bmatrix}\cdot$
$\displaystyle\begin{bmatrix}\boldsymbol{A}_{t}^{-\\!1}k(\boldsymbol{X}_{t},\boldsymbol{x}_{b})-\boldsymbol{A}_{t}^{-\\!1}k(\boldsymbol{X}_{t},\overline{\boldsymbol{X}})\left(k_{t}(\overline{\boldsymbol{X}})+\sigma^{2}\boldsymbol{I}_{t}\right)^{-\\!1}\left(k(\overline{\boldsymbol{X}},\boldsymbol{x}_{b})-k(\overline{\boldsymbol{X}},\boldsymbol{X}_{t})\boldsymbol{A}_{t}^{-\\!1}k(\boldsymbol{X}_{t},\boldsymbol{x}_{b})\right)\\\
\left(k_{t}(\overline{\boldsymbol{X}})+\sigma^{2}\boldsymbol{I}_{t}\right)^{-\\!1}\left(k(\overline{\boldsymbol{X}},\boldsymbol{x}_{b})-k(\overline{\boldsymbol{X}},\boldsymbol{X}_{t})\boldsymbol{A}_{t}^{-\\!1}k(\boldsymbol{X}_{t},\boldsymbol{x}_{b})\right)\end{bmatrix}$
$\displaystyle\sslash\text{{\emph{rearranging}}}$
$\displaystyle=k(\boldsymbol{x}_{a},\boldsymbol{x}_{b})-\begin{bmatrix}k(\boldsymbol{x}_{a},\boldsymbol{X}_{t})&k(\boldsymbol{x}_{a},\overline{\boldsymbol{X}})\end{bmatrix}\cdot$
$\displaystyle\begin{bmatrix}\boldsymbol{A}_{t}^{-\\!1}k(\boldsymbol{X}_{t},\boldsymbol{x}_{b})-\boldsymbol{A}_{t}^{-\\!1}k(\boldsymbol{X}_{t},\overline{\boldsymbol{X}})\left(k_{t}(\overline{\boldsymbol{X}})+\sigma^{2}\boldsymbol{I}_{t}\right)^{-\\!1}k_{t}(\overline{\boldsymbol{X}},\boldsymbol{x}_{b})\\\
\left(k_{t}(\overline{\boldsymbol{X}})+\sigma^{2}\boldsymbol{I}_{t}\right)^{-\\!1}k_{t}(\overline{\boldsymbol{X}},\boldsymbol{x}_{b})\end{bmatrix}$