# Numerically Stable Sparse Gaussian Processes
via Minimum Separation using Cover Trees

exander Terenin$*1,2$ David R. Burt$*1,3$ Artem Artemev$*2,5$ Seth Flaxman$4$
Mark van der Wilk$2$ Carl Edward Rasmussen$1,5$ Hong Ge$1$
$1$University of Cambridge $2$Imperial College London $3$MIT $4$University of
Oxford $5$Secondmind

###### Abstract

Gaussian processes are frequently deployed as part of larger machine learning
and decision-making systems, for instance in geospatial modeling, Bayesian
optimization, or in latent Gaussian models. Within a system, the Gaussian
process model needs to perform in a stable and reliable manner to ensure it
interacts correctly with other parts of the system. In this work, we study the
numerical stability of scalable sparse approximations based on inducing
points. To do so, we first review numerical stability, and illustrate typical
situations in which Gaussian process models can be unstable. Building on
stability theory originally developed in the interpolation literature, we
derive sufficient and in certain cases necessary conditions on the inducing
points for the computations performed to be numerically stable. For low-
dimensional tasks such as geospatial modeling, we propose an automated method
for computing inducing points satisfying these conditions. This is done via a
modification of the cover tree data structure, which is of independent
interest. We additionally propose an alternative sparse approximation for
regression with a Gaussian likelihood which trades off a small amount of
performance to further improve stability. We provide illustrative examples
showing the relationship between stability of calculations and predictive
performance of inducing point methods on spatial tasks.

$*$Equal contribution. Preprint.
Code available at: https://github.com/awav/conjugate-gradient-sparse-gp.

##  1 Introduction

Gaussian processes are a flexible framework and model class for learning
unknown functions. By way of being constructed in the language of Bayesian
learning, Gaussian process models provide an ability to incorporate prior
information into the model, and assess and propagate uncertainty in a
principled manner. This makes them well-suited for a wide variety of areas
where these capabilities are important, including statistical applications
such as spatial modeling cressie92 (8), and decision-making applications such
as Bayesian optimization snoek12 (44), sensor placement krause08 (32), and
active learning krause07 (31).

In many settings, the increased availability of data and need to accurately
model higher-resolution phenomena has led to a strong interest in working with
Gaussian processes at a larger scale. Unfortunately, classical Gaussian
process models generally scale cubically with training data size due to the
need to solve large linear systems of equations. This mismatch has led to a
longstanding and fruitful line of work on scalable Gaussian processes. In the
era of GPUs and automatic differentiation, two classes of scalable
approximations have been deployed within major Gaussian process software
packages, including _GPflow_ gpflow (35) and _GPyTorch_ gardner2018gpytorch
(17): those based on _inducing point methods_ quinonero2005unifying (40, 45,
22), and iterative methods such as the _(preconditioned) conjugate gradient_
algorithm gibbs1997efficient (18, 17).

Motivated by these advances, in this work we study a complementary question:
how can we _guarantee_ that the above algorithms run successfully, no matter
what kind of data they are provided with? This question is particularly
important in areas such as Bayesian optimization where the data is not
available in advance, and in latent variable models, where the data going into
the Gaussian process is reconstructed from auxiliary observations. We focus
chiefly on inducing point methods applied to geospatial or other low-
dimensional data, though the techniques we develop also reveal insights on how
the performance of conjugate-gradient-based approaches depends on data.

Our work complements existing analyses of numerical stability in Gaussian
processes. foster09 (15) study the subset of regressors approximation, and
describe in detail which numerical linear algebra routines one should use to
implement it in the most stable manner. basak21 (3) study maximum likelihood
estimation, observe that stability of log-determinant computations is
effectively controlled by the same quantities that control stability of linear
solves, and describe practical strategies for initializing and stopping
numerical optimization in a stable manner. In contrast, we describe how
different Gaussian process approximations and their hyperparameters affect the
numerical stability of the linear-algebraic routines which need to be
implemented, regardless of the details of how this implementation is done.

Our contributions are as follows. We (i) identify that mathematical tools
developed in the interpolation literature narcowich92 (37, 36, 42, 10) can be
used to show, in a large class of Gaussian process priors, that datasets
satisfying a _minimum separation_ criterion lead to linear systems which are
provably stable independent of data size. Using this insight, we (ii) develop
a clustering analogy and propose a technique to guarantee stability in sparse
approximations by using a refinement of the _cover tree_ data structure
beygelzimer06 (4) to automatically select inducing points that satisfy the
needed criterion while staying as close to the true data as possible. This
approximation is parameterized by a hyperparameter called the _spatial
resolution_ , which directly controls the tradeoff between approximation
accuracy and computational cost. Following this, we (iii) extend the
clustering analogy further, and propose an alternative, non-variational
inducing point approximation which trades off a small amount of performance to
further improve stability. Finally, we (iv) study error and computational
complexity of the proposed methods, and evaluate them on a number of examples.

##  2 Gaussian Processes

Let $X$ be a set. We say that a random function $f:X\rightarrow\mathbb{R}$ is
a _Gaussian process_ if, for any finite set of points $\boldsymbol{x}\in
X^{N}$, the random vector $f(\boldsymbol{x})$ is multivariate Gaussian. We
write $f\sim\operatorname{GP}(\mu,k)$, where
$\mu(\boldsymbol{x})=\operatorname{\mathbb{E}}(f(\boldsymbol{x}))$ is the
_mean function_ , and
$k(\boldsymbol{x},\boldsymbol{x}^{\prime})=\operatorname{Cov}(f(\boldsymbol{x}),f(\boldsymbol{x}^{\prime}))$
is the _covariance kernel_ , which determine the mean and covariance of the
corresponding multivariate Gaussians.

Define the likelihood
$\boldsymbol{y}=f(\boldsymbol{x})+\boldsymbol{\varepsilon}$ where
$\boldsymbol{\varepsilon}\sim\operatorname{N}(\boldsymbol{0},\mathbf{\Sigma})$.
Given training data $\boldsymbol{x},\boldsymbol{y}$, if we place a Gaussian
process prior $f\sim\operatorname{GP}(0,k)$, then the posterior
$f\mid\boldsymbol{y}$ is a Gaussian process with

$\displaystyle\operatorname{\mathbb{E}}(f\mid\boldsymbol{y})$
$\displaystyle=\mathbf{K}_{(\cdot)\boldsymbol{x}}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}+\mathbf{\Sigma})^{-1}\boldsymbol{y}$
$\displaystyle\operatorname{Cov}(f\mid\boldsymbol{y})$
$\displaystyle=\mathbf{K}_{(\cdot,\cdot^{\prime})}-\mathbf{K}_{(\cdot)\boldsymbol{x}}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}+\mathbf{\Sigma})^{-1}\mathbf{K}_{\boldsymbol{x}(\cdot^{\prime})}$
(1)

where $(\cdot)$ and $(\cdot^{\prime})$ denote arbitrary sets of test locations
and for sets $a,b$ we use $\mathbf{K}_{ab}$ to denote the kernel matrix formed
by evaluating $k$ at points in $a\times b$. We assume throughout that
$\mathbf{\Sigma}$ is diagonal. It is also possible to express the posterior in
the form of _pathwise conditioning_ wilson20 (58, 59) as

$(f\mid\boldsymbol{y})(\cdot)=f(\cdot)+\mathbf{K}_{(\cdot)\boldsymbol{x}}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}+\mathbf{\Sigma})^{-1}(\boldsymbol{y}-f(\boldsymbol{x})-\boldsymbol{\varepsilon})$
(2)

where equality holds in distribution. These expressions describe the true
posterior, whose computational costs are, classically, $\mathcal{O}(N^{3})$
owing to the Cholesky decomposition used to handle the matrix inverse. To
alleviate this, we now consider approximations.

###  2.1 Inducing Points

Many of the most widely used Gaussian process approximations are based on the
idea of _inducing points_ csato2002sparse (9, 43, 40). We present the
variational formulation of titsias09 (45, 22), using a pathwise conditioning
construction wilson20 (58, 59). Let $\boldsymbol{z}\in X^{M}$, and define the
_variational posterior_

$\displaystyle(f\mid\boldsymbol{u})(\cdot)$
$\displaystyle=f(\cdot)+\mathbf{K}_{(\cdot)\boldsymbol{z}}\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}^{-1}(\boldsymbol{u}-f(\boldsymbol{z}))$
$\displaystyle\boldsymbol{u}$
$\displaystyle\sim\operatorname{N}(\boldsymbol{m},\mathbf{S})$ (3)

where $\boldsymbol{m}$, $\mathbf{S}$ and $\boldsymbol{z}$ are approximation
parameters. We find these by minimizing the Kullback–Leibler divergence of the
approximation from the true posterior. By the chain rule for Kullback–Leibler
divergences, this quantity reduces to an easy-to-evaluate divergence between
finite-dimensional Gaussians matthews16 (34). Given $\boldsymbol{z}$, which
one can interpret as a learned set of approximate $\boldsymbol{x}$-values, the
optimal values for $\boldsymbol{m}$ and $\mathbf{S}$ can be solved for
analytically in closed form: see titsias09 (45) for details.

This approximation only requires factorizing a smaller kernel matrix, leading
to a cost of $\mathcal{O}(M^{3})$, a significant improvement if $M\ll N$.
Following hensman13 (22), the approximation also inherits variational methods’
compatibility with stochastic optimization and mini-batch methods.

###  2.2 Gaussian Processes for Geospatial Modelling

Gaussian processes models are widely used in spatial statistics both for
interpolation, and to allow for spatially correlated noise in spatial
regression models cressie92 (8). In low-dimensional settings, such as
geospatial data, inducing point approximations can lead to provably near-
linear computational scalability, particularly when input points, $x_{i}$, are
sampled from within a bounded region, and the Gaussian process prior is smooth
enough—see burt19 (6, 7). This makes them a reasonable choice for many
applications.

In geospatial modeling, specialized methods including multi-resolution
approaches katzfuss2017multi (26) or sparse-precision approaches vecchia1988
(51, 33, 27) often result in accurate approximations, especially when the
domain is much larger than the scale over which the function being modelled
varies. See, for instance, heaton2019case (21) for a comparison of a large
number of approaches on problems in this class. In such settings, low-rank
calculations very similar to those used in inducing point methods remain
central to the construction of a number of more tailored approaches, including
in particular multi-resolution methods katzfuss2017multi (26) and combinations
of low-rank and compactly supported methods sang2012fsa (41). A better
understanding of numerical stability in inducing point methods is therefore an
important step to better understanding when these more sophisticated methods
designed specifically for spatial data will be numerically stable.

###  2.3 Numerical Stability

On a computer, the calculations needed to deploy a Gaussian process, whether
exact or approximate, must be performed in floating-point arithmetic. Due to
the introduction of roundoff error, algorithms which solve linear systems
$\mathbf{A}^{-1}\boldsymbol{b}$ can fail if the system’s solution is too
sensitive to the numerical values of $\mathbf{A}$ or $\boldsymbol{b}$. We
assume throughout that $\mathbf{A}$ is symmetric positive definite, as will be
the case for the kernel matrices of interest. The key quantity used to
understand how numerically stable a linear system is the associated _condition
number_

$\operatorname{cond}(\mathbf{A})=\lim_{\varepsilon\rightarrow
0}\sup_{\left\lVert\boldsymbol{\delta}\right\rVert\leq\varepsilon\left\lVert\boldsymbol{b}\right\rVert}\frac{\left\lVert\mathbf{A}^{-1}(\boldsymbol{b}+\boldsymbol{\delta})-\mathbf{A}^{-1}\boldsymbol{b}\right\rVert_{2}}{\varepsilon\left\lVert\mathbf{A}^{-1}\boldsymbol{b}\right\rVert_{2}}=\big{\lVert}\mathbf{A}\big{\rVert}_{2}\big{\lVert}\mathbf{A}^{-1}\big{\rVert}_{2}=\frac{\lambda_{\max}(\mathbf{A})}{\lambda_{\min}(\mathbf{A})}$
(4)

of the matrix defining the system. Here, $\lambda_{\max}$ and $\lambda_{\min}$
are the maximum and minimum eigenvalues, respectively, and
$\left\lVert\cdot\right\rVert_{2}$ denotes the Euclidean norm and the
corresponding induced operator norm. A linear system’s condition number
quantifies how difficult it is to solve numerically. For a given floating-
point precision, if $\operatorname{cond}(\mathbf{A})$ is small enough and the
size of $\mathbf{A}$ is not too large, then Cholesky factorization is
guaranteed to succeed and return an accurate matrix square root.

######  Result 1.

Let $\mathbf{A}$ be a symmetric positive definite matrix of size $N\times N$.
Assume that $N>10$, that

$\operatorname{cond}(\mathbf{A})\leq\frac{1}{2^{-t}\times 3.9N^{3/2}}$ (5)

where $t$ is the length of the floating point mantissa, and that
$3N2^{-t}<0.1$. Then floating point Cholesky factorization will succeed,
producing a matrix $\mathbf{L}$ satisfying

$\displaystyle\mathbf{L}\mathbf{L}^{T}$ $\displaystyle=\mathbf{A}+\mathbf{E}$
$\displaystyle\left\lVert\mathbf{E}\right\rVert_{2}$ $\displaystyle\leq
2^{-t}\times 1.38N^{3/2}\left\lVert\mathbf{A}\right\rVert_{2}.$ (6)

###### Proof.

This follows from (kielbasinski87, 29, Corollary 2), and the relationship
between the Euclidean and Frobenius matrix norms. See also (wilkinson66, 57,
Theorem 2). ∎

For single precision arithmetic, we have $2^{-t}\approx 10^{-7.2}$, and for
double-precision arithmetic, we have $2^{-t}\approx 10^{-16}$. We therefore
see that well-conditioned systems of equations lead to numerically stable
Cholesky factorizations. Moreover, when an iterative algorithm is used to
solve a well-conditioned system, the iteration is often guaranteed to converge
quickly, leading to computational benefits. In particular, if one solves the
linear system with the conjugate gradient algorithm golub96 (19), then the
algorithm is guaranteed to converge in logarithmically many steps.

######  Result 2.

Let $\mathbf{A}$ be a positive semi-definite matrix of size $N\times N$, and
let $\varepsilon>0$. Then, in exact arithmetic, the conjugate gradient
algorithm for solving $\mathbf{A}^{-1}\boldsymbol{b}$ converges to within
$\varepsilon$ of the true solution, with respect to the norm induced by
$\mathbf{A}$, in

$\mathcal{O}\left(\sqrt{\operatorname{cond}(\mathbf{A})}\log\frac{\operatorname{cond}(\mathbf{A})\|\boldsymbol{b}\|}{\varepsilon}\right)$
(7)

total steps.

###### Proof.

The claim follows from (golub96, 19, Theorem 10.2.6) by a direct calculation
given in Section A.2. ∎

Working with stable linear systems is therefore crucial for both inducing
point methods which rely on Cholesky factorizations, and for algorithms based
on conjugate gradients. In the latter case, _preconditioning_ is often used to
improve stability and accelerate convergence. Preconditioning can be viewed as
a way to partially solve the linear system, resolving its unstable components.
State-of-the-art preconditioners in many areas of applied mathematics, such as
multigrid methods for elliptic partial differential equations e11 (14, 60),
rely on detailed properties of their respective settings for construction and
analysis. A key step in designing provably effective preconditioners for
Gaussian processes is therefore to understand how the model and data influence
their resulting linear systems: we focus on this and defer ideas on
preconditioning to future work.

###  2.4 Instability in Gaussian Process Models

Figure 1: Here we illustrate the condition numbers and resulting linear system
error for the Kac–Murdock–Szegö matrix, which is the kernel matrix of an
exponential kernel on a regularly-spaced one-dimensional grid. We vary the
spacing $r$ and respective correlation $\rho=\exp(-r)$ of neighboring points,
and the size $n$ of the matrix, with $\rho=0.999$ used in cases where $n$
varies, and $n=256$ used in cases where $\rho$ varies. We plot the condition
number, computed numerically using an eigenvalue factorization, along with its
theoretical lower and upper bounds. Then, we generate random vectors
$\boldsymbol{v}\sim\operatorname{N}(\boldsymbol{0},\mathbf{I})$, compute
$\boldsymbol{u}=\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}\boldsymbol{v}$,
solve for
$\boldsymbol{v}=\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}^{-1}\boldsymbol{u}$
numerically, and plot the median error norm over $10^{4}$ samples, along with
25% and 75% quantiles. We observe that condition numbers asymptotically grow
as $\rho\to 1$, but not as $N\to\infty$, and increase hand-in-hand with
numerical error.

To begin, we first observe that, generically, condition numbers of kernel
matrices need not be well-behaved. To illustrate what can go wrong, consider
the one-dimensional exponential kernel

$k(x,x^{\prime})=\exp(-|x-x^{\prime}|).$ (8)

Suppose that $x$ is a time series which lies on a regular one-dimensional grid
with spacing $r$, so that $x_{i}=ri$. Then the resulting kernel matrix is a
_Kac–Murdock–Szegö matrix_ trench2001properties (47, 11), meaning it takes the
form

$\displaystyle\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}=\begin{pmatrix}1&\rho&\rho^{2}&\dots&\rho^{n-1}\\\
\rho&1&\rho&\dots&\rho^{n-2}\\\ \vdots&\vdots&\ddots&\ddots&\vdots\\\
\rho^{n-1}&\rho^{n-2}&\rho^{n-3}&\dots&1\end{pmatrix}$ (9)

for $\rho=\exp(-r)$. To compute a posterior Gaussian process under exact
observations, the condition number of the matrix we need to invert to compute
the posterior satisfies (trench2001properties, 47, p. 9) the inequality

$\frac{(1+\rho)^{2}}{(1-\rho)^{2}}\leq\operatorname{cond}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}})\leq\frac{(1+\rho)^{2}+2\rho\epsilon}{(1-\rho)^{2}-2\rho\epsilon}$
(10)

where $\varepsilon=\frac{\pi^{2}}{(N+1)^{2}}$ and for the upper bound we have
assumed $(1-\rho)^{2}>2\rho\epsilon$. For large $N$, the value $\varepsilon$
is small so the upper and lower bounds essentially match. As $r\rightarrow 0$,
we have $\rho\rightarrow 1$, and the condition number diverges to infinity for
any given value of $N$. On the other hand, if the spread of the data, $r$ is
fixed, then regardless of how much data is collected, $\rho$ and the condition
number of $\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}$ both remain bounded. In
cases where data is observed under a Gaussian likelihood, one instead needs to
invert $\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}+\mathbf{\Sigma}$. If
$\mathbf{\Sigma}$ is diagonal, the minimum eigenvalue of
$\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}+\mathbf{\Sigma}$ is bounded below
by the minimum diagonal entry of $\mathbf{\Sigma}$, which we will refer to as
$\sigma^{2}$ by analogy to the homoscedastic noise case when
$\mathbf{\Sigma}=\sigma^{2}\mathbf{I}$. Unfortunately, this is not enough to
ensure stability: if one samples the time series more densely, then as
$r\rightarrow 0$ and $N\rightarrow\infty$, the condition number diverges to
infinity for all values of $\sigma^{2}$ due to the growth of the largest
eigenvalue of $\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}$. We now observe that
the behavior illustrated by this example is neither specific to the
exponential kernel, nor to data with algebraic structure arising from grids,
nor limited to one dimension.

propositionPropCovOpLimit Let $k$ be a continuous stationary kernel, and let
$\boldsymbol{x}$ be randomly sampled from a uniform distribution on some
compact subset $X\subseteq\mathbb{R}^{d}$. Define the covariance operator

$\displaystyle\mathcal{K}:L^{2}(X)$ $\displaystyle\rightarrow L^{2}(X)$
$\displaystyle\mathcal{K}:\phi$
$\displaystyle\mapsto\int_{X}\phi(x)k(x,\cdot)\mathop{}\\!\mathrm{d}x$ (11)

acting on the Hilbert space of (equivalence classes of) square-integrable
functions. Then the eigenvalues of $\mathcal{K}$ are countable, non-negative,
admit a maximum, and can be ordered to form a non-increasing sequence whose
limit is zero. Fix and $n>0$, then as as $N\rightarrow\infty$, each the $n$th
largest eigenvalue of $\frac{1}{N}\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}$
converges almost surely to the corresponding eigenvalue of $\mathcal{K}$. As a
consequence,
$\operatorname{cond}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}})\rightarrow\infty$
almost surely.

###### Proof.

This follows from (koltchinskii2000, 30, Corollary 3.3), after a bit of
technical reformulation and some calculations, which we perform in Section A.
∎

Figure 1 implies that $\lambda_{\max}$ grows linearly with $N$, while Cauchy’s
interlacing theorem together with an elementary argument shows that
$\lambda_{\min}\rightarrow 0$, causing the condition number of
$\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}$ to grow without bound. From the
one-dimensional example and Figure 1, we therefore see that to control the
condition numbers of kernel matrices, we need to ensure they are sufficiently
far away from this ill-conditioned limiting case. We conclude that if the data
is sampled too closely relative to the kernel’s length scale, computing a
Cholesky factorization of
$\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}+\mathbf{\Sigma}$ numerically is not
possible for large matrices. We proceed to explore this in more general cases
in the sequel.

##  3 Numerical Stability in Scalable Gaussian Process Approximations

To build towards Gaussian process approximations which are numerically stable,
we study what properties of models and data lead to numerical stability of
kernel matrices. We focus in particular on how their condition numbers depend
on the number of input points under general regimes.

###  3.1 Numerical Stability in Sparse Approximations via Minimum Separation

In both exact and sparse Gaussian processes, the linear systems that need to
be solved arise from kernel matrices
$\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}$ and
$\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}$, which depend on the prior and
data. In sparse Gaussian processes, however, $\boldsymbol{z}$ is not fixed: it
is a variational parameter, and is usually selected to maximize approximation
accuracy. We now study whether we can use this freedom to select
$\boldsymbol{z}$ in a way that ensures stability, while maintaining
performance. To this end, we introduce the following notion.

######  Definition 3.

Let $\boldsymbol{z}$ be the inducing points. Define the _separation radius_

$\operatorname{sep}(\boldsymbol{z})=\min_{i\neq j}\left\lVert
z_{i}-z_{j}\right\rVert.$ (12)

We now show that separation is closely connected with numerical stability. For
this, we need a mild regularity condition on the kernel.

[Spatial decay]assumptionAsmSpatialDecay Let $k$ be a kernel on
$X\subseteq\mathbb{R}^{d}$. We say that $k$ has _spatial decay_ if there is a
decreasing function $\psi:[0,\infty)\rightarrow[0,\infty)$ such that for all
$x,x^{\prime}\in X$ we have

$\displaystyle|k(x,x^{\prime})|$ $\displaystyle\leq\psi(\left\lVert
x-x^{\prime}\right\rVert)$ $\displaystyle\psi(m)$
$\displaystyle=\mathcal{O}\left(\frac{1}{m^{d}\log(m)^{2}}\right).$ (13)

In words, spatial decay requires that if we consider points far apart, the
covariance between the function at those points must tend to zero sufficiently
quickly with the distance between them. This assumption is satisfied for a
very large set of kernels, including potentially non-stationary kernels, as
well as the squared exponential kernel and all Matérn kernels, whose decay is
controlled above and below by products of an exponential with polynomials. We
now state the main results.

######  Proposition 4.

Let $X\subseteq\mathbb{R}^{d}$, and let $k$ satisfy spatial decay. Then there
is a constant $C_{\max}^{k,\delta}$ such that for any $M$ and any
$\boldsymbol{z}$ of size $M$ with
$\operatorname{sep}(\boldsymbol{z})\geq\delta>0$, we have

$\lambda_{\max}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}})\leq
C_{\max}^{k,\delta}.$ (14)

###### Proof.

Section A, Proposition 19. ∎

From a technical perspective, the claim is proven by combining Gershgorin’s
Circle Theorem with a packing argument. Arguments of this kind have a long
history in the interpolation literature—see for instance narcowich92 (37, 36,
42, 10). In cases where the Gaussian process approximation requires inversion
of $\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda}$, from this we
can immediately conclude a condition number bound.

######  Corollary 5.

Under the conditions of Proposition 4, for diagonal $\mathbf{\Lambda}$,
letting $\Lambda_{\max}$ and $\Lambda_{\min}$ denote its respective maximum
and minimum entry, we have

$\operatorname{cond}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})\leq\frac{C_{\max}^{k,\delta}+\Lambda_{\max}}{\Lambda_{\min}}.$
(15)

###### Proof.

Combine Proposition 4 with the triangle inequality for operator norm, noting
that for positive semi-definite matrices the operator norm is equal to the
largest eigenvalue. ∎

This generalizes immediately to the case where $\mathbf{\Lambda}$ is not
diagonal, as long as $\Lambda_{\max}$ and $\Lambda_{\min}$ are replaced with
its maximum and minimum eigenvalues. Since in many practical cases
$\mathbf{\Lambda}$ is diagonal, we retain this presentation for simplicity. In
the noiseless case, we can also conclude a condition number bound, but now
additionally require stationarity to do so. We will need the following result.

######  Result 6.

Let $X\subseteq\mathbb{R}^{d}$, and let $k$ be stationary and continuous. Then
there is a constant $C_{\min}^{k,\delta}>0$ such that for any $M$ and any
$\boldsymbol{z}$ of size $M$ with
$\operatorname{sep}(\boldsymbol{z})\geq\delta>0$, we have

$\lambda_{\min}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}})\geq
C_{\min}^{k,\delta}.$ (16)

###### Proof.

(wendland04, 56, Theorem 12.3). See Section A for a sketch of the main
argument. ∎

We can combine this the eigenvalue upper bound to bound the condition number,
which we present for completeness.

######  Corollary 7.

Under the conditions of Proposition 4 and Result 6, there is a constant
$C_{\operatorname{cond}}^{k,\delta}>1$ such that for any $M$ and any
$\boldsymbol{z}$ of size $M$ with $\operatorname{sep}(\boldsymbol{z})<\delta$,
we have

$\operatorname{cond}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}})\leq
C_{\operatorname{cond}}^{k,\delta}.$ (17)

###### Proof.

Combine Proposition 4 and Result 6 with
$\operatorname{cond}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}})=\frac{\lambda_{\max}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}})}{\lambda_{\min}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}})}$,
taking
$C_{\operatorname{cond}}^{k,\delta}=\frac{C_{\max}^{k,\delta}}{C_{\min}^{k,\delta}}$.
∎

These results therefore reveal what properties the set of inducing points
$\boldsymbol{z}$ needs to have in order to yield numerically stable linear
systems. Separation distance and spatial decay rule out the bad behavior shown
in Figure 1, by preventing the limiting covariance operator from controlling
properties of the corresponding kernel matrices. The analysis shown is based
on the techniques presented by (wendland04, 56, Ch. 12) for proving Result 6,
which originate in the polynomial interpolation literature narcowich92 (37,
42). For non-stationary kernels, mirroring Proposition 4, we conjecture that a
similar bound to Result 6 for $\lambda_{\min}$ holds for a wide class of
kernels, but requires a different argument which avoids the use of Fourier
analysis. To conclude, we observe that for most kernels, if one allows for an
arbitrarily large number of inducing points and no noise is added to the
kernel matrix, then the minimum separation condition is necessary to ensure
bounded condition numbers.

propositionPropSeparationNecessary Let $z_{m}\in X$, $m=1,..,\infty$, be a
sequence, and let $k$ be a Lipschitz continuous kernel on $X\times X$. Define
$\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}^{(M)}\in\mathbb{R}^{M\times M}$ to
be the kernel matrix formed by evaluating $k$ at $z_{1},..,z_{M}$. If there
exists a $C>0$ such that for all $M$,
$\operatorname{cond}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}^{(M)})\leq C$,
then $z_{m}$ satisfy minimum separation.

###### Proof.

Section A.5. ∎

###  3.2 Inducing Points for Geospatial Data via Cover Trees

To obtain an algorithmic solution that guarantees numerical stability, we
study how to construct inducing points which satisfy a user-specified minimum
separation radius. Since separation is defined purely in terms of the inducing
points, we need a second criterion to quantify how well the inducing points
summarize a given dataset. A clear choice is to place the inducing points as
close to the data as possible.

######  Definition 8.

Let $\boldsymbol{x}$ be the data, and $\boldsymbol{z}$ be the inducing points.
Define the _spatial resolution_

$\displaystyle\operatorname{res}_{\boldsymbol{x}}(\boldsymbol{z})$
$\displaystyle=\max_{i=1,..,N}\min_{j=1,..,M}\left\lVert
x_{i}-z_{j}\right\rVert.$ (18)

In simple terms, the spatial resolution is the maximum distance from a data
point to the closest inducing point. This notion is closely related to _fill
distance_ which appears in the the analysis of posterior contraction rates
kanagawa18 (25), but instead involves both a dataset $\boldsymbol{x}$ and
inducing points $\boldsymbol{z}$. Among all sets of inducing points with a
given spatial resolution, finding ones which satisfy separation is always
possible. More precisely, if we take the inducing points $\boldsymbol{z}$ to
be the centers of balls which make up a maximal separated subset of radius
$\operatorname{sep}(\boldsymbol{z})=\varepsilon$, the resulting choice
automatically has spatial resolution
$\operatorname{res}_{\boldsymbol{x}}(\boldsymbol{z})=\varepsilon$. A similar
observation applies if one takes $\boldsymbol{z}$ to be an optimal covering of
$\boldsymbol{x}$—see Section A for details. This shows that the notions of
separation distance and spatial resolution are closely connected with optimal
covering and packing numbers—see (vershynin18, 52, Section 4.2).

We now study numerical techniques for computing inducing points that
simultaneously satisfy separation and covering guarantees. Our goal is to
obtain a practical near-linear-time method producing a set of inducing points
with small spatial resolution, to ensure accurate approximations of the
posterior, and large separation to ensure stability. To achieve this, we
propose an approach we term the _$R$ -neighbor cover tree_, which is a
refinement of the _cover tree_ data structure kibriya07 (28, 24) originally
proposed by beygelzimer06 (4) for efficient nearest neighbor search. This
approach is well-suited to cases where the inputs are low-dimensional such as
geospatial data. We note that the proposed algorithm will also perform well in
cases when the dataset has low _intrinsic dimension_ —see beygelzimer06
(4)—even if the dataset is embedded in a high dimensional space.

(a) Place covering nodes

(b) Compute $R$-neighbors

(c) Repeat steps (a) and (b) recursively

Figure 2: Two iterations of Algorithm 9 on a tree with $L=3$ total levels.
Given a region, the algorithm first (a) computes a covering of the region by
picking data points one-by-one from those not yet covered. Then, it (b)
computes the $R$-neighbors of each node and uses this to efficiently compute a
Voronoi tessellation of the region. Finally, the algorithm (c) repeats the
process recursively.

A cover tree is a tree where nodes are points associated with metric balls.
The root node is a ball which covers all of the data $\boldsymbol{x}$. Each
parent node contains the centers of its child nodes, which consist of balls of
smaller radius than the parent. At each level, the tree enforces separation
and covering properties similar to those in Definition 3 and Definition 8,
with the separation distance and spatial resolution proportional to the radius
of balls at that depth. beygelzimer06 (4) present a depth-first algorithm for
building such a tree.

Since we are interested in inducing point selection with a fixed spatial
resolution $\varepsilon$, rather than nearest neighbors search, we seek a
different set of guarantees than those provided by beygelzimer06 (4), or by
improved versions such as that of izbicki15 (24). The most important
difference is that, in the inducing point case, we want separation to hold
_globally_ in every level of the tree, and not only for children with a common
parent. To do this, we present a modified breadth-first construction, which is
of independent interest, and conceptually works as follows.

  1. 1.

Initialize the root node $z_{1,1}$ at the mean of the data, and assign all
data to it.

  2. 2.

Loop over tree depth $\ell$, starting at $\ell=0$.

    1. (a)

Loop over parent nodes $z_{\ell,p}$.

      1. i.

Select a point $z^{\prime}$ from the parent’s assigned data.

      2. ii.

_Optional:_ compute the local average of the parent node’s assigned data
around $z^{\prime}$, and set $z^{\prime}$ to be this average, as long as it is
not too close to another node.

      3. iii.

Create a child node $z_{\ell+1,c}$ centered at $z^{\prime}$.

      4. iv.

Reassign all points near $z_{\ell+1,c}$ in nodes at level $\ell$ to
$z_{\ell+1,c}$.

    2. (b)

_Optional:_ reassign all points at level $\ell+1$ according to the Voronoi
partition of $z_{\ell+1,c}$.

A formal description is given in Algorithm 9. The key insight behind this top-
down, breadth-first construction is that, for a dataset with a given intrinsic
dimensionality, it can be performed in near-linear time. The reason for this
is that every step is _local_ , and can be computed by searching only the
children of nodes sufficiently close to the current parent node. More
precisely, we show that if the parent node’s radius is at level $\ell$ with
radius $R$ and the total depth of the tree is $L$, then only nodes within a
distance of $4(1-1/2^{(L-\ell)})R\leq 4R$ need to be searched. The number of
such nodes is bounded, can be tracked recursively, and is controlled by the
intrinsic dimensionality of the data.

The two optional steps do not alter the algorithm’s complexity or its
guarantees, but improve practical performance. The local averaging step we
propose is similar to the Lloyd’s iteration of $K$-means, and allows the tree
to use fewer nodes by placing them in-between data points, rather than exactly
on top of data points. The Voronoi repartitioning step, originally proposed by
izbicki15 (24), makes the tree better-balanced. We now prove the algorithm
works as intended.

Algorithm 9 Cover Tree Inducing Points

1: input spatial resolution $\varepsilon>0$ and dataset $\boldsymbol{x}$.
Define notation $B_{z}(R)=\\{x\in X:\left\lVert z-x\right\rVert\leq R\\}$.

2: Initialize root node $z_{0,1}=\frac{1}{N}\sum_{i=1}^{N}x_{i}$, assigned
data $\mathcal{A}_{1,1}=\boldsymbol{x}$, $R$-neighbors
$\mathcal{R}_{1,1}=\\{1\\}$, as well as constants
$d_{\max}=\max_{i=1,..,N}\left\lVert z_{1,1}-x_{i}\right\rVert$,
$L=\left\lceil\log_{2}\frac{d_{\max}}{\varepsilon}\right\rceil$, $M_{0}=1$,
and $R_{0}=2^{L}\varepsilon$.

3: for tree depth level $\ell=1,..,L$ do

4: Initialize number of nodes $M_{\ell}=0$, child node index $c=1$, and radius
$R_{\ell}=\frac{1}{2}R_{\ell-1}$.

5: for parent node index $p=1,..,M_{\ell-1}$ do

6: Initialize child nodes $\mathcal{C}_{\ell-1,p}=\varnothing$ and,
_optionally_ , assigned data copy
$\mathcal{A}^{\prime}_{\ell-1,p}=\mathcal{A}_{\ell-1,p}$.

7: while assigned data $\mathcal{A}_{\ell-1,p}\neq\varnothing$ do

8: Choose an arbitrary data point $\zeta\in\mathcal{A}_{\ell-1,p}$.

9:  _Optional:_ compute local average
$\zeta^{\prime}=\frac{1}{|\mathcal{Z}|}\sum_{z\in\mathcal{Z}}z$, where
$\mathcal{Z}=\mathcal{A}_{\ell-1,p}\cap B_{\zeta}(R_{\ell})$.

10: if
$\min_{z\in\mathcal{Z}^{\prime}}\left\lVert\zeta-z\right\rVert>R_{\ell}$ where
$\mathcal{Z}^{\prime}=\bigcup_{r\in\mathcal{R}_{\ell-1,p_{\ell,c}}}\mathcal{C}_{\ell-1,p}$
then set $\zeta=\zeta^{\prime}$.

11: Create child node $z_{\ell,c}=\zeta$.

12: for parent $R$-neighbor index $r$ in $\mathcal{R}_{\ell-1,p}$ do

13: Update assigned data
$\mathcal{A}_{\ell,c}=\mathcal{A}_{\ell,c}\cup(\mathcal{A}_{\ell-1,r}\cap
B_{z_{\ell,c}}(R_{\ell}))$,
$\mathcal{A}_{\ell-1,r}=\mathcal{A}_{\ell-1,r}\setminus
B_{z_{\ell,c}}(R_{\ell})$.

14: Update children
$\mathcal{C}_{\ell-1,p}=\mathcal{C}_{\ell-1,p}\cup\\{c\\}$, parent
$p_{\ell,c}=p$, index $c=c+1$, and $M_{\ell}=M_{\ell}+1$.

15: for child node index $\varsigma=1,..,c$ do

16: Compute child $R$-neighbors
$\mathcal{R}_{\ell,\varsigma}=\bigcup_{r\in\mathcal{R}_{\ell-1,p_{\ell,c}}}\mathcal{C}_{\ell-1,r}\cap
B_{z_{\ell,c}}(4(1-1/2^{(L-\ell)})R_{\ell})$.

17:  _Optional:_ recompute
$\mathcal{A}_{\ell,\varsigma}=\\{x\in\bigcup_{r\in\mathcal{R}_{\ell-1,p_{\ell-1,c}}}\mathcal{A}^{\prime}_{\ell,r}:\left\lVert
x-z_{\varsigma}\right\rVert\leq\left\lVert
x-z_{\varsigma^{\prime}}\right\rVert,\forall\varsigma^{\prime}\neq\varsigma\\}$.

18: return $\boldsymbol{z}=\\{z_{\ell,m}:\ell=L\\}$.

theoremThmCoverTree For a given target spatial resolution $\varepsilon>0$, and
dataset with $\bar{x}=\frac{1}{N}\sum_{i=1}^{N}x_{i}$ and
$d_{\max}=\max_{i}\left\lVert x_{i}-\bar{x}\right\rVert$, the cover tree
inducing point algorithm, with or without the optional steps, produces a tree
with $L=\left\lceil\log_{2}\frac{d_{\max}}{\varepsilon}\right\rceil+1$ levels,
terminates in

$\mathcal{O}\left(P_{B(4)}^{\operatorname{ext},1/2}\left(P_{B(1)}^{\operatorname{ext},1/2}\right)^{2}N\log\frac{d_{\max}}{\varepsilon}\right)$
(19)

steps, where $B(r)$ denotes a ball of radius $r$, and
$P_{X}^{\operatorname{ext},\delta}$ is the external packing number, namely the
maximum number of disjoint balls of radius $\delta$ that one can choose so
that their centers lie in $X$. Moreover, the algorithm guarantees that the
nodes at each level $\ell$ satisfy

$\displaystyle\operatorname{res}_{\boldsymbol{x}}(\boldsymbol{z}_{\ell})$
$\displaystyle\leq 2^{L-\ell}\varepsilon$
$\displaystyle\operatorname{sep}(\boldsymbol{z}_{\ell})$ $\displaystyle\geq
2^{L-\ell}\varepsilon.$ (20)

###### Proof.

Section C. ∎

This gives an efficient and practical algorithm for choosing inducing points
which guarantees the properties needed in Corollary 7 hold. The algorithm’s
output is fairly-sharp: compared to an optimal covering, the cover tree
construction may require a slightly larger number of inducing points, but
achieves the same separation guarantee, irrespective of the number of levels.
Using this, we obtain the following numerical stability guarantee.

######  Corollary 10.

Let $\boldsymbol{x}$ be a dataset of size $N$, and let $\varepsilon>0$ Let
$\boldsymbol{z}$ be a set of inducing points computed from $\boldsymbol{x}$ by
the cover tree algorithm with target spatial resolution $\varepsilon$. Let $k$
be stationary, continuous, and satisfy spatial decay. Then
$\operatorname{cond}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}})$ is bounded
independently of $N$.

###### Proof.

Combine Algorithm 9 with Corollary 7. ∎

For kernels which possess a length scale parameter, the desired spatial
resolution in a given problem might depend on the value of the length scale.
This parameter is often set by training the model via maximum marginal
likelihood, and can therefore change as training progresses. In such settings,
one can use the hierarchical structure of the tree to dynamically adjust the
number of inducing points as the length scale changes. This can be done by
switching to inducing points obtained from an intermediate level of the tree
if the length scale becomes coarser, or by computing an extra level of the
tree if it becomes finer. The simplest way to do this is to adjust the model
after initial training is complete, since this avoids non-smooth changes in
the inducing point approximation during training, and fine-tune the model
after the update if necessary.

###  3.3 Additional Stability via the Clustered-data Inducing Point
Approximation

In the preceding sections, we studied properties of data needed to obtain a
stability guarantee for inducing point methods. We now study a complementary
question: if selecting $\boldsymbol{z}$ using a cover tree, can the inducing
point approximation be modified to improve stability further? To proceed, we
begin by defining a low-rank analog bui2017unifying (5, 39, 1) of the opper09
(38) approximation

$\displaystyle(f\mid\boldsymbol{u})(\cdot)$
$\displaystyle=f(\cdot)+\mathbf{K}_{(\cdot)\boldsymbol{z}}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})^{-1}(\boldsymbol{u}-f(\boldsymbol{z})-\boldsymbol{\epsilon})$
$\displaystyle\boldsymbol{\epsilon}\sim\operatorname{N}(\boldsymbol{0},\mathbf{\Lambda})$
(21)

where $\boldsymbol{z}$, diagonal noise matrix $\mathbf{\Lambda}$, and the
inducing mean value $\boldsymbol{u}$ are now variational parameters. Since its
optimization objective is non-convex in the variational parameters, this
approximation is generally used with natural gradients adam2021dual (1) or
stochastic optimization panos18 (39, 49). Similar to the earlier inducing
point approximation of titsias09 (45, 22), it recovers the true posterior when
$M=N$, $\boldsymbol{z}=\boldsymbol{x}$, $\mathbf{\Lambda}=\mathbf{\Sigma}$,
and $\boldsymbol{u}=\boldsymbol{y}$.

This approximation is more numerically stable thanks to the diagonal matrix
$\mathbf{\Lambda}$ panos18 (39, 49), which lower bounds the minimum eigenvalue
of $\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda}$ by
$\min_{i=1,..,M}\Lambda_{ii}+\lambda_{\min}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}})>0$.
For finite datasets, this is a strict improvement in stability over
approximations involving only $\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}$.
Since $\mathbf{\Lambda}$ effectively controls the width of the posterior error
bars, the degree to which this benefit manifests itself depends on the data,
with more uncertainty leading to more numerically stable linear systems.

(a) Original Gaussian process

(b) Clustered-data approx.

(c) Equivalent sparse process

Figure 3: The idea behind the clustered-data approximation of Figure 3. To
construct this approximation, we move each data point in $\boldsymbol{x}$ to
its nearest cluster, transforming (a) to (b). We can then reinterpret and
represent (b) in a sparse manner by merging repeated data points, replacing
$\boldsymbol{y}$ with the cluster means $\boldsymbol{u}$, and re-weighting
each data point according to cluster size—obtaining (c), which is equal in
distribution to (b).

Before studying convergence further, we observe that this approximation’s
relationship with the true posterior is deeper than might first meets the eye.
For a given set of inducing points $\boldsymbol{z}$ and dataset
$\boldsymbol{x}$, define the _nearest inducing point clustering_ and _cluster
size_ maps

$\displaystyle\operatorname{cl}(x)$
$\displaystyle=\operatorname*{\arg\min}_{z_{j},\,j=1,..,M}\left\lVert
x-z_{j}\right\rVert$ $\displaystyle N_{\operatorname{cl}}(x)$
$\displaystyle=|\\{x^{\prime}\in\boldsymbol{x}:\operatorname{cl}(x)=\operatorname{cl}(x^{\prime})\\}|.$
(22)

With these notions, we show that, under certain choices of the variational
parameters, the Opper–Archambeau approximation is exactly the correct
posterior under a perturbed dataset, where $\boldsymbol{x}$ is replaced with
$\operatorname{cl}(\boldsymbol{x})$ and all other data and model parameters
are unchanged.

propositionPropInducingPointRepr Let
$u_{j}=\frac{1}{N_{\operatorname{cl}}(z_{j})}\sum_{\operatorname{cl}(x_{i})=z_{j}}y_{i}$.
The Bayesian models defined by $f\sim\operatorname{GP}(0,k)$ and

$\displaystyle y_{i}\mid f$
$\displaystyle\sim\operatorname{N}(f(\operatorname{cl}(x_{i})),\sigma^{2})$
$\displaystyle u_{i}\mid f$
$\displaystyle\sim\operatorname{N}\left(f(z_{i}),\frac{\sigma^{2}}{N_{\operatorname{cl}}(z_{i})}\right).$
(23)

admit respective posterior distributions $f\mid\boldsymbol{y}$ and
$f\mid\boldsymbol{u}$ which are equal in distribution.

###### Proof.

Section B. ∎

In cases where the spatial resolution of Definition 8 is small—which the cover
tree guarantees—we can use Figure 3 to accurately estimate the correct
variational parameters without needing to perform optimization. This is done
by choosing $\boldsymbol{z}$ to be the clusters, choosing $\boldsymbol{u}$ to
be the cluster mean, and choosing the entries of the heteroskedastic noise
matrix $\mathbf{\Lambda}$ to be equal to the original noise variance divided
by the cluster size. We call this the _clustered data inducing point
approximation_ , and illustrate it in Figure 3.

The explicit description of the approximation made by use of inducing points
in Figure 3 enables us to better understand how to select them
algorithmically. To develop this line of thought in more detail, we proceed to
study how to choose inducing points to optimize the tradeoff between
approximation error and computational costs using conjugate gradients.

In many applications, employing this approximation involves computing its
marginal likelihood objective, or related quantities such as prior density
evaluations or Kullback–Leibler divergences. One of the noteworthy properties
of this approximation is that the Kullback–Leibler divergence between it and
the true prior simplifies into

$D_{\operatorname{KL}}(q\mid\mid
p)=\frac{1}{2}\ln\frac{|\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda}|}{|\mathbf{\Lambda}|}-\frac{1}{2}\operatorname{tr}((\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})^{-1}\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}})+\frac{1}{2}\boldsymbol{v}^{T}\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}\boldsymbol{v}$
(24)

where
$\boldsymbol{v}=(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})^{-1}\boldsymbol{u}$,
and the minus sign in the trace term appears due to cancellation with the
usual dimension term. All of the linear systems appearing in this term are
numerically stable. This shows, for instance, that if one uses variational
inference with the clustered-data approximation, _all computations_ needed to
train the Gaussian process are numerically stable, not just the ones for
obtaining predictive means and variances.

##  4 Experiments

We investigate the performance of the clustered-data inducing point
approximation on a series of tests designed to illustrate both the behavior of
its components and the overall picture. Each experiment focuses on a different
aspect of the method, including approximation error, scalability, and behavior
under varying quantities of data and inducing points. Full details regarding
the experimental setup are given in Section D.

###  4.1 Spatial Resolution and Empirical Approximation Error

To better understand how different spatial resolutions resolve the tradeoff
between approximation error and performance, we performed a number of
numerical experiments designed to highlight this tradeoff. To quantify
approximation error, we chose the Wasserstein distance, because it is finite
in the setting of interest, readily estimatable numerically, and controls
expectations of a large class of posterior functionals villani08 (53).

To assess these quantities, we generated synthetic data by sampling $N=1000$
input points uniformly on a hypercube of dimension $1,2,4$ and $8$, and then
sampling from the prior process, which uses a squared exponential kernel. We
then varied the spatial resolution, computed inducing points using the cover
tree algorithm, and empirically estimated the Wasserstein distance between the
approximate and true Gaussian processes. To better contextualize results,
since the cover tree construction produces a variable number of inducing
points, we also recorded the number of inducing points produced for spaces of
different dimension.

Results can be seen in Figure 4. Immediately, we see that behavior of the
cover tree depends on the dimension of the space the data lives on. In higher
dimension, more inducing points are needed to construct a separated covering
with the same spatial resolution. On the other hand, we see spatial resolution
directly controls approximation accuracy in Wasserstein distance, with finer
spatial resolutions more accurate than coarser ones, and the effect being
similar in different dimensions.

Next, to assess the consequence of different spatial resolutions on stability
and corresponding difficulty of solving the resulting linear systems, we
assess empirical performance of conjugate gradients in the settings of
interest. We adopt the same data-generation process used in the preceding
variant, but this time focus on the linear systems. We again vary spatial
resolution, but now compute the numerical condition numbers of the resulting
kernel matrices, and run the conjugate gradient algorithm to convergence,
obtaining the number of iterations needed for convergence.

In Figure 4, we see that finer spatial resolutions generally result in more
expensive-to-solve and less-stable linear systems. The effect varies according
to dimension, and is least pronounced in $d=8$: in high dimensions, when
sampling uniformly on a hypercube, the inputs tend to end up far away from
each other, leading to less extreme eigenvalue s. In lower dimensions, where
the posterior is more concentrated, the effect of spatial resolution on
condition numbers is more pronounced. In all cases, higher condition numbers
lead to more iterations conjugate gradient iterations needed. In total, we see
that spatial resolution is a directly interpretable quantity which facilitates
the tradeoff between approximation accuracy with stability and computational
cost.

Figure 4: We illustrate how spatial resolution affects approximation accuracy
and computational cost. We see that increasing the spatial resolution leads to
using fewer inducing points, dropping rapidly in low dimension, and slower in
higher dimensions. On the other hand, the error in approximation, as measured
by Wasserstein distance between the approximate posterior and exact posterior,
increases as the spatial resolution increases. Computational costs follow the
opposite relationship: finer spatial resolutions are less stable and more
expensive to compute, as evidenced by condition numbers and the number of
iterations needed for conjugate gradients to converge.

###  4.2 Comparison to Alternative Inducing Point Selection Schemes

Next, to understand how different inducing point selection methods perform, we
performed an experiment to assess the tradeoffs in terms of performance,
stability, and computational cost. We compare five methods for selecting
inducing points: the cover tree approach described in Algorithm 9, the online
inducing point selection algorithm proposed by galyfajou21 (16), the partial
pivoted Cholesky approach studied by foster09 (15), along with baselines
consisting of ordinary $K$-means, $K$-means++ arthur2006k (2), optimizing the
inducing points directly by minimizing the respective Kullback–Leibler
divergence, and uniform random sampling from the training data.

To assess performance and stability, we evaluated the methods in terms of
their produced kernel matrix condition numbers, and three sparse Gaussian
process regression performance metrics, namely root mean squared error,
negative log predictive density, and the evidence lower bound. Together, these
give a view of both training and test-set performance, as well as variational
approximation error. We chose the sparse approximation of titsias09 (45) owing
to its widespread use, working with the _GPFlow_ , whose linear algebra is
implemented in a stable manner using the $V$-method suggested by foster09
(15). We use the squared exponential kernel and a jitter value of $10^{-6}$.
We evaluated the different methods on two datasets: the two-dimensional _East
Africa_ geospatial dataset studied by ton2018 (46), and the four-dimensional
_Combined Cycle Power Plant_ dataset from the UCI Machine Learning Repository
tufekci2014power (48, 12). These give a view of performance in both geospatial
and non-geospatial settings, while ensuring all methods run successfully and
can therefore be meaningfully compared.

Figure 5: We compare how different inducing point selection algorithms vary
according to performance and stability, using the geospatial _East Africa_
dataset, and higher-dimensional _Power_ dataset. The number of inducing points
produced by algorithms whose output is variable is shown rounded to the
nearest increment. Note that the computational complexity of these methods
differs: cover tree is $\mathcal{O}(N\log N)$, online and the $K$-means
variants are $\mathcal{O}(NM)$, partial pivoted Cholesky is
$\mathcal{O}(NM^{2})$, and optimizing the inducing points is
$\mathcal{O}(NM^{2})$ per optimization iteration. We see that using more
inducing points results in higher performance, but reduced stability, with the
precise effect varying according to dataset, dimensionality, and inducing
point selection method.

Results can be seen in Figure 5. Immediately, we see that different methods
produce different tradeoffs between performance and stability. These tradeoffs
need not be achieved in an optimal manner: one might intuitively expect the
ordinary $K$-means baseline to produce inducing points that are well-separated
due to its clustering nature, and therefore result in stable linear systems.
This doesn’t happen: instead, $K$-means tends to simultaneously achieve the
worst performance and worst numerical stability among all methods except
possibly uniform subsampling. In comparison, $K$-means++, which uses a
sophisticated initialization scheme that promotes separation, produces much
better performance and stability than $K$-means, though still not as well as
methods with specified separation guarantees. This shows that one must
explicitly encourage separation to obtain good performance-stability tradeoffs
via such approaches.

The strongest predictive performance is achieved by optimizing the inducing
points directly by minimizing the respective Kullback–Leibler divergence using
gradient descent. The resulting inducing points, however, can be as
numerically unstable as those produced by $K$-means or uniform subsampling.
Since optimizing inducing point locations requires one to specify an
initialization, this suggests it may be advantageous to choose this
initialization carefully, using one of the stable inducing point selection
algorithms, to minimize potential numerical failures during the optimization
process.

On the geospatial example, again somewhat-surprisingly, all methods with
stability guarantees produced near-identical condition numbers and performance
for a given number of inducing points. The main difference with these methods
is therefore their computational cost, which in low dimension are
$\mathcal{O}(N\log N)$ for the cover tree, $\mathcal{O}(NM)$ for the online
approach and $K$-means++, and $\mathcal{O}(NM^{2})$ for partial pivoted
Cholesky. The online method’s complexity can, for some kernels, be improved to
$\mathcal{O}(N\log N)$ by for instance precomputing a tree-based data
structure that enables fast nearest-neighbor lookups, in which case the
resulting variant starts to resemble the final loop of the cover tree
construction algorithm. Overall, we conclude, in this setting, that
significant numerical stability improvements can be obtained without reducing
performance—and, using the cover tree, without paying excessive computational
costs.

On the higher-dimensional non-geospatial example, no method of inducing point
selection achieves a clear improvement on top of uniform subsampling of the
data in terms of performance. On the other hand, the cover tree, online
method, $K$-means++, and partial pivoted Cholesky all achieve better numerical
stability compared to $K$-means, uniform sampling, and optimization, with
partial pivoted Cholesky being the most stable. Due to the curse of
dimensionality, the $\mathcal{O}(NM^{2})$ computational costs of partial
pivoted Cholesky are similar or favorable to the other approaches, making this
method competitive. In this setting, we therefore also conclude that
significant numerical stability improvements can be obtained—however, unlike
in the lower-dimensional setting, the computational costs needed to do so can
become expensive.

In summary, we conclude from the comparisons that different approaches’
strengths are most effective in different settings. In low-dimensional
geospatial problems, among all non-iterative methods, the cover tree algorithm
provides excellent performance with guaranteed stability, while running in
near-linear time. The online approach performs comparably, and can be more
appropriate for settings where data arrives sequentially. $K$-means++, which
promotes but does not enforce separation, performs slightly worse. In high-
dimensional problems, the partial pivoted Cholesky algorithm results in strong
performance, provided the desired number of inducing points is small enough
that running it is practical. One can obtain the best performance by
optimizing the inducing points directly, at cost of having to run an iterative
procedure, and producing worse stability. In both classes of problems, our
experiments suggest that methods with stability guarantees can offer a better
performance-stability tradeoff compared to naïve methods.

###  4.3 Reliability in Floating Point Precision: Geospatial Illustrative
Example

$\varepsilon=0.09$
$M=902$

$\varepsilon=0.06$
$M=1934$

$\varepsilon=0.03$
$M=6851$

Figure 6: Here we illustrate the effect of spatial resolution $\varepsilon$
and resulting number of inducing points $M$ on the stability and performance
of the Gaussian process model for the clustered-data approximation (blue), and
sparse variational Gaussian process baseline (orange). We see that coarser
spatial resolutions result in blurrier predictions and, in turn, higher RMSE.
While the sparse variational Gaussian process baseline achieves better
performance on coarser resolutions, its Cholesky factorization eventually
fails, requiring us to increase jitter—this is displayed by the dashed line.
The clustered-data approximation with conjugate gradients, on the other hand,
runs reliably in all cases, achieving slightly worse but overall comparable
performance.

One of the appealing properties of using numerically stable approaches in
variational inference is that, in principle, they can enable one to use a
larger set of inducing points than customarily used. This is for two reasons:
(i) the quadratic cost of conjugate gradients and (ii) the ability to save
memory by effectively running in floating point precision. To illustrate that
using more inducing points can improve performance on applied examples, we
train the variational approximation of titsias09 (45) and with the clustered-
data approximation of Section 3.3 on the East Africa land surface temperature
dataset wan02 (54, 55, 46). In both cases, we use floating-point precision,
and the same inducing points obtained using a cover tree.

One of the defining characteristics of this dataset is that it has a
relatively small number of data points per unit length scale ball, due to the
fine-scale variation in land surface temperature. This renders approximations
based on a sparse set of inducing points less accurate, unless the number of
inducing points approaches the order of the size of the dataset. By varying
the spatial resolution from relatively long ($\varepsilon=0.09$) to very short
($\varepsilon=0.03$) we see both qualitative and quantitative improvement,
shown in Figure 6.

When using larger sets of inducing points and running in floating point
precision, numerical stability becomes increasingly crucial for running
successfully. For the variational approximation of titsias09 (45), this
manifests itself as Cholesky factorization failure. A common heuristic for
alleviating this in the Gaussian process literature is to add _jitter_ to the
kernel matrix: this is done by replacing
$\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}$ with
$\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\varepsilon\mathbf{I}$ for a small
value $\varepsilon>0$. This increases all of the eigenvalues of the matrix by
$\varepsilon$, and ensures that the smallest eigenvalue cannot be arbitrarily
close to $0$. We found that using a larger jitter value compared to the
_GPflow_ recommended default was necessary when working with the finest
spatial resolution of $\varepsilon=0.03$, which produced $M=6851$ inducing
points, indicating that this separation distance alone was too small to ensure
stability in this example.

In comparison, the clustered-data approximation does not require jitter, since
it only involves inversion of the matrices
$\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda}$ which include the
diagonal matrix $\mathbf{\Lambda}$. This makes it a more reliable approach,
since it does not involve risk of failure, but results in slightly worse
predictive performance. In settings where the Gaussian process is trained in
an automated manner without human oversight, or where data is collected online
and not available in advance, reliability can be a bigger concern than
performance. In such cases, our results indicate that the clustered-data
approximation can be an alternative worth considering.

In this example, numerical stability is first and foremost a consequence of
the minimum kernel matrix eigenvalue
$\lambda_{\min}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}})$ being too close to
zero. Both introducing jitter and using the clustered-data approximation can
alleviate this, but cannot help with issues arising from the maximum
eigenvalue being too large, which must eventually occur if many inducing
points accumulate in a bounded region. Our experiment indicates that this does
not occur on the scales considered, but might instead occur when working with
substantially larger sets of inducing points, for instance in the presence of
sparse or otherwise structured kernel matrices durrande19 (13), which can
often be inverted at substantially larger scale.

###  4.4 Numerical Stability and Dataset Size

Figure 7: Here we illustrate how numerical stability and performance change
when data size and approximation accuracy are varied simultaneously. Viewing a
single plot top-to-bottom, we see that for a fixed data size, increasing the
number of inducing points produces better performance and very slightly
improves stability. Viewing the same plot left-to-right, we see that for a
fixed level of variational approximation expressiveness, using more inducing
inducing points produces better performance but decreases stability. Viewing
the plot diagonally, we see that increasing data size and variational
approximation expressiveness simultaneously results in better performance but
worse stability, with different methods choosing tradeoffs in favor of either
performance or stability, depending on the details of how they are designed.

The preceding sections show that, in the examples considered, numerical
stability tends to deteriorate as approximation accuracy increases, as
measured by the number of inducing points. On the other hand, in many
settings, the more data one has, the more inducing points one needs to
accurately represent the true posterior. One can therefore also think of the
number of inducing points as measuring the expressive capacity of the
variational family. This naturally raises the question: what happens to
numerical stability when the number of data points grows, and expressive
capacity grows with it?

To understand this, we decided to examine how stability and predictive
performance change depending on data size and the number of inducing points
simultaneously. For this, we randomly sampled subsets from the East Africa
dataset used in the preceding illustrative example, and computed both the
respective kernel matrix condition number, and predictive root-mean-squared
error.

Results can be seen in Figure 7. For a fixed data size (Figure 7, top-to-
bottom), we see that using more inducing points improves performance, but
makes stability worse, with the tradeoff determined by the individual method:
numerically stable methods favor stability, whereas optimizing the inducing
points favors performance but quickly results in unstable linear algebra. On
the other hand, for a fixed number of inducing points (Figure 7, left-to-
right), using more data always improves performance, and can actually slightly
improve stability. In light of Corollary 7, this makes sense, because using
more data provides a larger region within which the algorithms can place a
fixed number of inducing points, allowing them to be spaced further apart from
one another.

If we grow both the number of data points and the number of inducing points at
the same time while keeping their ratio approximately fixed (Figure 7,
diagonal, bottom-left to top-right), we see that numerical stability
deteriorates as data size and the number of inducing points increase
simultaneously. The severity of the tradeoff depends on the algorithm used:
stable methods provide a tradeoff which favors numerical stability, while
optimized inducing points favor performance and become unstable faster.
$K$-means achieves the worst of both worlds, having a level of performance
comparable to stable methods, and a level of numerical stability comparable to
optimizing the inducing points. This highlights the value of methods designed
to produce a favorable performance-stability tradeoff.

##  5 Conclusion

We study numerical stability of variational approximations used for Gaussian
process models. We survey and introduce a number of results which, in in
total, enable one to conclude that the _minimum separation distance_ between
inducing points controls stability properties of the numerical linear algebra
needed to train variational approximations. We introduce two practical
techniques for ensuring stability: inducing points constructed by a _cover
tree_ , and the _clustered-data approximation_. We provide examples
illustrating the tradeoff between performance, stability, and dataset size.
Compared to standard inducing point placement methods, we find that approaches
that promote minimum separation produce a favorable performance-stability
tradeoff.

## Acknowledgments

DRB was supported by the Qualcomm Innovation Fellowship. SF was supported by
the EPSRC (EP/V002910/2). We further acknowledge support from Huawei Research
and Development, the Imperial College COVID-19 Research Fund, and UK Research
and Innovation.

## References

  * (1) Vincent Adam, Paul Chang, Mohammad Emtiyaz E Khan and Arno Solin  “Dual Parameterization of Sparse Variational Gaussian Processes”  In _Advances in Neural Information Processing Systems_ , 2021
  * (2) David Arthur and Sergei Vassilvitskii  “$K$-means++: The advantages of careful seeding”  In _ACM Symposium on Discrete Algorithms_ , 2007
  * (3) Subhasish Basak, Sébastien Petit, Julien Bect and Emmanuel Vazquez  “Numerical issues in maximum likelihood parameter estimation for Gaussian process interpolation”  In _International Conference on Machine Learning, Optimization, and Data Science_ , 2021, pp. 116–131  Springer
  * (4) Alina Beygelzimer, Sham Kakade and John Langford  “Cover trees for nearest neighbor”  In _International Conference on Machine Learning_ , 2006, pp. 97–104
  * (5) Thang D Bui, Josiah Yan and Richard E Turner  “A unifying framework for Gaussian process pseudo-point approximations using power expectation propagation”  In _The Journal of Machine Learning Research_ , 2017
  * (6) David R Burt, Carl E Rasmussen and Mark Wilk  “Rates of Convergence for Sparse Variational Gaussian Process Regression”  In _International Conference on Machine Learning_ , 2019
  * (7) David R. Burt, Carl Edward Rasmussen and Mark Wilk  “Convergence of Sparse Variational Inference in Gaussian Processes Regression”  In _Journal of Machine Learning Research_ , 2020
  * (8) Noel Cressie  “Statistics for Spatial Data”  Wiley, 1992
  * (9) Lehel Csató and Manfred Opper  “Sparse on-line Gaussian processes”  In _Neural Computation_ , 2002
  * (10) Benedikt Diederichs and Armin Iske  “Improved estimates for condition numbers of radial basis function interpolation matrices”  In _Journal of Approximation Theory_ , 2019
  * (11) Murray Dow  “Explicit inverses of Toeplitz and associated matrices”  In _Australian and New Zealand Industrial and Applied Mathematics Journal_ , 2002
  * (12) Dheeru Dua and Casey Graff  “UCI Machine Learning Repository”, 2017  URL: http://archive.ics.uci.edu/ml
  * (13) Nicolas Durrande et al.  “Banded Matrix Operators for Gaussian Markov Models in the Automatic Differentiation Era”  In _Artificial Intelligence and Statistics_ , 2019
  * (14) Weinan E  “Principles of Multiscale Modeling”  Cambridge University Press, 2011
  * (15) Leslie Foster et al.  “Stable and Efficient Gaussian Process Calculations.”  In _Journal of Machine Learning Research_ , 2009
  * (16) Théo Galy-Fajou and Manfred Opper  “Adaptive Inducing Points Selection For Gaussian Processes”  In _Workshop on Continual Learning at the International Conference on Machine Learning_ , 2021
  * (17) Jacob Gardner et al.  “Gpytorch: Blackbox matrix-matrix Gaussian process inference with GPU acceleration”  In _Advances in Neural Information Processing Systems_ , 2018
  * (18) Mark Gibbs and David Mackay  “Efficient Implementation of Gaussian Processes”, 1997
  * (19) Gene H Golub and Charles F Van Loan  “Matrix Computations”  Johns Hopkins University Press, 1996
  * (20) Daniel H Greene and Donald Ervin Knuth  “Mathematics for the Analysis of Algorithms”  Springer, 1990
  * (21) Matthew J Heaton et al.  “A case study competition among methods for analyzing large spatial data”  In _Journal of Agricultural, Biological and Environmental Statistics_ 24 Springer, 2019, pp. 398–425
  * (22) James Hensman, Nicolò Fusi and Neil D Lawrence  “Gaussian processes for big data”  In _Uncertainty in Artificial Intelligence_ , 2013
  * (23) Michael F. Hutchinson  “A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines”  In _Communications in Statistics – Simulation and Computation_ , 1989
  * (24) Mike Izbicki and Christian Shelton  “Faster cover trees”  In _International Conference on Machine Learning_ , 2015
  * (25) Motonobu Kanagawa, Philipp Hennig, Dino Sejdinovic and Bharath K Sriperumbudur  “Gaussian processes and kernel methods: A review on connections and equivalences”  In _arXiv:1807.02582_ , 2018
  * (26) Matthias Katzfuss  “A multi-resolution approximation for massive spatial datasets”  In _Journal of the American Statistical Association_ , 2017
  * (27) Matthias Katzfuss and Joseph Guinness  “A General Framework for Vecchia Approximations of Gaussian Processes”  In _Statistical Science_ , 2021
  * (28) Ashraf M Kibriya and Eibe Frank  “An empirical comparison of exact nearest neighbour algorithms”  In _European Conference on Principles of Data Mining and Knowledge Discovery_ , 2007
  * (29) Andrzej Kiełbasiński  “A note on rounding-error analysis of Cholesky factorization”  In _Linear Algebra and its applications_ 88, 1987
  * (30) Vladimir Koltchinskii and Evarist Giné  “Random matrix approximation of spectra of integral operators”  In _Bernoulli_ , 2000
  * (31) Andreas Krause and Carlos Guestrin  “Nonmyopic active learning of Gaussian processes: an exploration-exploitation approach”  In _International Conference on Machine learning_ , 2007
  * (32) Andreas Krause, Ajit Singh and Carlos Guestrin  “Near-optimal sensor placements in Gaussian processes: Theory, efficient algorithms and empirical studies.”  In _Journal of Machine Learning Research_ , 2008
  * (33) Finn Lindgren, Håvard Rue and Johan Lindström  “An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 73.4, 2011, pp. 423–498
  * (34) Alexander G. G. Matthews, James Hensman, Richard Turner and Zoubin Ghahramani  “On sparse variational methods and the Kullback–Leibler divergence between stochastic processes”  In _Artificial Intelligence and Statistics_ , 2016
  * (35) Alexander G. G. Matthews et al.  “GPflow: A Gaussian Process Library using TensorFlow”  In _Journal of Machine Learning Research_ , 2017
  * (36) F.J. Narcowich, N. Sivakumar and J.D. Ward  “On Condition Numbers Associated with Radial-Function Interpolation”  In _Journal of Mathematical Analysis and Applications_ , 1994
  * (37) Francis J Narcowich and Joseph D Ward  “Norm estimates for the inverses of a general class of scattered-data radial-function interpolation matrices”  In _Journal of Approximation Theory_ , 1992
  * (38) Manfred Opper and Cédric Archambeau  “The variational Gaussian approximation revisited”  In _Neural Computation_ , 2009
  * (39) Aristeidis Panos, Petros Dellaportas and Michalis K Titsias  “Fully scalable Gaussian processes using subspace inducing inputs”  In _arXiv:1807.02537_ , 2018
  * (40) Joaquin Quiñonero-Candela and Carl Edward Rasmussen  “A unifying view of sparse approximate Gaussian process regression”  In _Journal of Machine Learning Research_ , 2005
  * (41) Huiyan Sang and Jianhua Z. Huang  “A full scale approximation of covariance functions for large spatial data sets”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 2012
  * (42) Robert Schaback  “Error estimates and condition numbers for radial basis function interpolation”  In _Advances in Computational Mathematics_ , 1995
  * (43) Matthias W Seeger, Christopher KI Williams and Neil D Lawrence  “Fast forward selection to speed up sparse Gaussian process regression”  In _Artificial Intelligence and Statistics_ , 2003
  * (44) Jasper Snoek, Hugo Larochelle and Ryan P Adams  “Practical Bayesian optimization of machine learning algorithms”  In _Advances in Neural Information Processing Systems_ , 2012
  * (45) Michalis Titsias  “Variational learning of inducing variables in sparse Gaussian processes”  In _Artificial Intelligence and Statistics_ , 2009
  * (46) Jean-Francois Ton, Seth Flaxman, Dino Sejdinovic and Samir Bhatt  “Spatial mapping with Gaussian processes and nonstationary Fourier features”  In _Spatial Statistics_ 28, 2018, pp. 59–78
  * (47) William F Trench  “Properties of some generalizations of Kac–Murdock–Szegö matrices”  In _Contemporary Mathematics_ , 2001
  * (48) Pınar Tüfekci  “Prediction of full load electrical power output of a base load operated combined cycle power plant using machine learning methods”  In _International Journal of Electrical Power and Energy Systems_ , 2014
  * (49) Mark Wilk, Artem Artemev and James Hensman  “Improved Inverse-Free Variational Bounds for Sparse Gaussian Processes”  In _Symposium on Advances in Approximate Bayesian Inference_ , 2022
  * (50) Ramon Handel  “Probability in high dimension”  Princeton University, 2014
  * (51) A.. Vecchia  “Estimation and Model Identification for Continuous Spatial Processes”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 1988
  * (52) Roman Vershynin  “High-dimensional probability”  Cambridge University Press, 2018
  * (53) Cédric Villani  “Optimal Transport: old and new”  Springer, 2008
  * (54) Zhengming Wan, Yulin Zhang, Qincheng Zhang and Zhao-liang Li  “Validation of the land-surface temperature products retrieved from Terra Moderate Resolution Imaging Spectroradiometer data”  In _Remote Sensing of Environment_ , 2002
  * (55) Daniel J Weiss et al.  “An effective approach for gap-filling continental scale remotely sensed time-series”  In _ISPRS Journal of Photogrammetry and Remote Sensing_ , 2014
  * (56) Holger Wendland  “Scattered Data Approximation”  Cambridge University Press, 2004
  * (57) J H Wilkinson  “A Priori Error Analysis of Algebraic Processes”  In _Proceedings of the International Congress of Mathematicians_ , 1966, pp. 119–129
  * (58) James T. Wilson et al.  “Efficiently Sampling Functions from Gaussian Process Posteriors”  In _International Conference on Machine Learning_ , 2020
  * (59) James T. Wilson et al.  “Pathwise Conditioning of Gaussian Processes”  In _Journal of Machine Learning Research_ , 2021
  * (60) Jinchao Xu and Ludmil Zikatanov  “Algebraic Multigrid Methods”  In _Acta Numerica_ , 2017

##  A Eigenvalues and Condition Numbers

###  A.1 Covariance Operator Limits

*

###### Proof.

Define the (random) matrix $\mathbf{\tilde{K}}_{\boldsymbol{x}\boldsymbol{x}}$
to be zero on the diagonal and equal to
$\mathbf{\tilde{K}}_{\boldsymbol{x}\boldsymbol{x}}$ on the off-diagonal. Let
$\Pi$ denotes the set of all permutations of the natural numbers. For an
operator $\mathcal{A}$ with discrete spectrum, we view the spectrum as an
element of $\ell^{p}/\Pi$, denoted by $\lambda(\mathcal{A})$. By convention,
we extend the spectrum of a symmetric, positive definite matrix $\mathbf{M}$
to an element of $\ell^{p}/\Pi$ by appending zeros, denoted by
$\lambda(\mathbf{M})$. For $p\in(1,\infty]$ define the the $\ell^{p}$
rearrangement distance by,

$d_{p}(x,y)=\inf_{\pi\in\Pi}\left(\sum_{i=1}^{\infty}|x_{i}-y_{\pi(i)}|^{p}\right)^{1/p},$
(25)

with the unusual convention that $p=\infty$ corresponds to the element-wise
supremum. One can show that $d_{p}$ defines a pseudo-metric on $\ell^{p}$ and
a metric on $\ell^{p}/\Pi$: the only property that is tricky to check is
triangle inequality, which is shown in (koltchinskii2000, 30, p. 116), in the
case $p=2$, which immediately generalizes to other $p$. Since $k$ is
stationary,
$\mathbf{\tilde{K}}_{\boldsymbol{x}\boldsymbol{x}}=\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}-c\mathbf{I}$
for some $c>0$. Therefore,

$d_{\infty}\left(\lambda(\tfrac{1}{N}\mathbf{\tilde{K}}_{\boldsymbol{x}\boldsymbol{x}}),\lambda(\tfrac{1}{N}\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}})\right)=\frac{c}{N}\rightarrow
0.$ (26)

Since $k$ is bounded, the conditions of (koltchinskii2000, 30, Corollary 3.3)
are satisfied, so

$d_{2}\left(\lambda(\tfrac{1}{N}\mathbf{\tilde{K}}_{\boldsymbol{x}\boldsymbol{x}}),\lambda(\mathcal{K})\right)\rightarrow
0.$ (27)

and therefore by the standard $(\ell^{2},\ell^{\infty})$ inequality

$d_{\infty}\left(\lambda(\tfrac{1}{N}\mathbf{\tilde{K}}_{\boldsymbol{x}\boldsymbol{x}}),\lambda(\mathcal{K})\right)\rightarrow
0$ (28)

which implies that by the triangle inequality

$d_{\infty}\Big{(}\lambda(\tfrac{1}{N}\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}),\lambda(\mathcal{K})\Big{)}\rightarrow
0.$ (29)

We have

$\Big{|}\lambda_{1}(\tfrac{1}{N}\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}})-\lambda_{1}(\mathcal{K})\Big{|}\leq
d_{\infty}\Big{(}\lambda(\tfrac{1}{N}\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}),\lambda(\mathcal{K})\Big{)}\rightarrow
0,$ (30)

so by the reverse triangle inequality

$\lambda_{1}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}})=N\lambda_{1}(\mathcal{K})+o(N).$
(31)

As $\mathcal{K}$ is not the zero operator and is positive definite, this is
$\Omega(N)$. On the other hand, there exists an $N_{0}$ such that for all
$n>N_{0}$, we have $\lambda_{n}(\mathcal{K})\leq\varepsilon/2$. Also, there
exists an $N_{0}^{\prime}$ such that for all $n\geq N_{0}^{\prime}$

$d_{\infty}(\lambda(\tfrac{1}{N}\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}),\lambda(\mathcal{K}))\leq\varepsilon/2.$
(32)

Choose $M_{0}=\max(N_{0},N_{0}^{\prime})$. Then for all $n\geq M_{0}$,

$\frac{1}{N}\lambda_{N}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}})\leq
d_{\infty}(\lambda(\tfrac{1}{N}\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}),\lambda(\mathcal{K}))+\lambda_{n}(\mathcal{K})\leq\varepsilon,$
(33)

and so $\lambda_{N}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}})=o(N)$. In fact,
we can say much more: as the data is compactly supported, for any
$\varepsilon>0$ and for all $N$ sufficient large there exists a $1\leq i,j\leq
N$ such that $k(x_{i},x_{j})\geq(1-\varepsilon)k(x_{i},x_{i})$. By Cauchy’s
Interlacing Theorem,

$\displaystyle\lambda_{N}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}})\leq\lambda_{2}\begin{pmatrix}k(x_{i},x_{i})&k(x_{i},x_{j})\\\
k(x_{j},x_{i})&k(x_{j},x_{j})\end{pmatrix}\leq k(x_{i},x_{i})\varepsilon.$
(34)

As $\varepsilon$ is arbitrary, we conclude
$\lambda_{N}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}})\rightarrow 0$. ∎

###  A.2 Convergence of Conjugate Gradients

######  Result 11.

Let $\boldsymbol{v^{*}}$ satisfy
$\mathbf{A}\boldsymbol{v^{*}}=\boldsymbol{b}$. For any $t\in\mathbb{N}$ and
$\boldsymbol{v_{0}}\in\mathbb{R}^{n}$ and $\boldsymbol{v_{t}}$ denote the
solution found by running conjugate gradients on this system of equations with
initial vector $\boldsymbol{v_{0}}$, then

$\|\boldsymbol{v^{*}}-\boldsymbol{v_{t}}\|_{\mathbf{A}}\leq
2\left(1-\frac{2}{\sqrt{\operatorname{cond}(\mathbf{A})}+1}\right)^{t}\|\boldsymbol{v^{*}}-\boldsymbol{v_{0}}\|_{\mathbf{A}}.$
(35)

###### Proof.

(golub96, 19, Theorem 10.2.6). ∎

######  Corollary 12.

A sufficient condition for
$\|\boldsymbol{v^{*}}-\boldsymbol{v_{t}}\|_{\mathbf{A}}\leq\varepsilon$ is

$t\geq\left(\log\Big{(}1+\tfrac{2}{\sqrt{\operatorname{cond}(\mathbf{A})}+1}\Big{)}\right)^{-1}\log\frac{2\|\boldsymbol{v^{*}}-\boldsymbol{v_{0}}\|_{\mathbf{A}}}{\varepsilon}=\left(\tfrac{\sqrt{\operatorname{cond}(\mathbf{A})}}{2}+\mathcal{O}(1)\right)\log\frac{2\|\boldsymbol{v^{*}}-\boldsymbol{v_{0}}\|_{\mathbf{A}}}{\varepsilon}$
(36)

where the implicit constant is between $\frac{1}{2}$ and $\frac{3}{2}$.

###### Proof.

The first inequality is immediate from inverting the bound in Result 11. The
asymptotic formula can be derived directly from applying the inequalities
$1-\frac{1}{x}\leq\log x\leq x-1$ to yield

$\frac{\sqrt{\operatorname{cond}(\mathbf{A})}+1}{2}\leq\left(\log\Big{(}1+\tfrac{2}{\sqrt{\operatorname{cond}(\mathbf{A})}+1}\Big{)}\right)^{-1}\leq\frac{\sqrt{\operatorname{cond}(\mathbf{A})}+3}{2}.$
(37)

∎

###  A.3 Upper Eigenvalue Estimates

Here we prove the eigenvalue upper bound under separation and spatial decay
via a packing argument. We begin with a basic eigenvalue estimate.

######  Result 13 (Gershgorin’s Circle Theorem).

Let $\mathbf{K}$ be a positive semi-definite symmetric matrix with eigenvalues
$\lambda(\mathbf{K})=\\{\lambda_{i}\\}$. Define the _Gershgorin radius_
$R_{i}=\sum_{j\neq i}|k_{ij}|$ to be the sum of the off-diagonal entries of
$\mathbf{K}$. Then
$\lambda(\mathbf{K})\in\bigcup_{i}[k_{ii}-R_{i},k_{ii}+R_{i}]$.

###### Proof.

(golub96, 19, Theorem 8.1.3). ∎

We can use this result to get upper eigenvalue estimates as

$\lambda_{\max}\leq\max_{i}k_{ii}+R_{i}$ (38)

given a set of points that are separated and a kernel satisfying spatial
decay. We now introduce these notions.

######  Definition 14.

We say that a set of points $\boldsymbol{x}\subset\mathbb{R}^{d}$ are
$\delta$-_separated_ if, for every $x_{i}\neq x_{j}$ we have
$\|x_{i}-x_{j}\|\geq\delta>0$.

######  Definition 15.

We define the ball of radius $\delta$-centered at $x\in\mathbb{R}^{d}$ to be
the open set,

$B_{x,\delta}=\\{x^{\prime}\in X:\|x-x^{\prime}\|<\delta\\}.$ (39)

We write $B_{\delta}=B_{0,\delta}$.

*

The argument for the maximum eigenvalue bound is essentially that if we have
spatial decay and uniformly bounded variance, then the kernel matrix has
bounded diagonal and bounded Gershgorin radius, both independent of $N$.

######  Definition 16 (Covering and packing numbers).

Let $X\subset\mathbb{R}^{d}$ and let $\delta>0$. Define the following.

  1. 1.

A subset $C\subseteq C^{\prime}\subseteq\mathbb{R}^{d}$ is an _external cover_
of $X$ in $C^{\prime}$ if for every $x\in X$ there is a $c\in C$ such that
$\left\lVert x-c\right\rVert\leq\delta$.

  2. 2.

A subset $C\subseteq X$ is an _internal cover_ of $X$ if for every $x\in X$
there is a $c\in C$ such that $\left\lVert x-c\right\rVert\leq\delta$.

  3. 3.

A subset $P\subseteq X$ is an _external packing_ of $X$ if any distinct
$x,x^{\prime}\in P$ satisfy $\left\lVert x-x^{\prime}\right\rVert>\delta$.

  4. 4.

A subset $P\subseteq X$ is an _internal packing_ of $X$ if it is an external
packing, and moreover satisfies $B_{x,\delta}\subseteq X$ for all $x\in P$.

Define the internal and external covering numbers
$C^{\operatorname{int},\delta}_{X}$ and
$C^{\operatorname{ext},\delta}_{X,C^{\prime}}$ to be the smallest cardinality
of such sets, and the internal and external packing numbers
$P^{\operatorname{int},\delta}_{X}$ and $P^{\operatorname{ext},\delta}_{X}$ to
be the largest cardinality of such sets.

######  Lemma 17.

Let $A_{m,\delta}=B_{0,(m+1)\delta}\setminus
B_{0,m\delta}=\\{x\in\mathbb{R}^{d}:m\delta\leq\left\lVert
x\right\rVert<(m+1)\delta\\}$ be an annulus with external radius $(m+1)\delta$
and internal radius $m\delta$. Then we have

$\displaystyle\operatorname{vol}(B_{\delta})$
$\displaystyle=\dfrac{\pi^{d/2}}{\Gamma(\frac{d}{2}+1)}\delta^{d}$
$\displaystyle\operatorname{vol}(A_{m,\delta})=\frac{\pi^{d/2}}{\Gamma(\frac{d}{2}+1)}((m+1)^{d}-m^{d})\delta^{d}.$
(40)

We also have for $m>0$

$(m+2)^{d}-m^{d}\leq 2dm^{d-1}.$ (41)

###### Proof.

For (41), expand the first term on the left hand side with the Binomial
Theorem. The lead term cancels with the second term on the left hand side,
while the right hand side is the next term. The remaining terms are non-
negative, which when dropped give the inequality. ∎

######  Lemma 18.

For $m\geq 1$,

$\left(\frac{1}{\delta}\right)^{d}\frac{\operatorname{vol}(A_{m,\delta})}{\operatorname{vol}(B_{\delta})}\leq
C^{\operatorname{ext},\delta}_{A_{m,\delta}}\leq
P^{\operatorname{ext},\delta}_{A_{m,\delta}}\leq
2d\left(\frac{2}{\delta}\right)^{d}\left(m-\frac{1}{2}\right)^{d-1}.$ (42)

and $P^{\operatorname{ext},\delta}_{A_{0,\delta}}=1$.

###### Proof.

The first part of the claim follows immediately from (vanhandel14, 50, Lemma
5.13). To bound $P^{\operatorname{ext},\delta}_{A_{m,\delta}}$ for $m\geq 1$,
use (vershynin18, 52, Proposition 4.2.12), along with the fact that the
Minkowski sum of an annulus with a ball is another annulus, to bound

$P^{\operatorname{ext},\delta}_{A_{m,\delta}}\leq\frac{\operatorname{vol}(B_{0,(m+\frac{3}{2})\delta}\setminus
B_{0,(m-\frac{1}{2})\delta})}{\operatorname{vol}(B_{\delta/2})}=\frac{\left(m+\frac{3}{2}\right)^{d}-\left(m-\frac{1}{2}\right)^{d}}{\left(\frac{\delta}{2}\right)^{d}}\leq
2d\left(\frac{2}{\delta}\right)^{d}\left(m-\frac{1}{2}\right)^{d-1}.$ (43)

For the case $m=0$, $A_{0,\delta}=B_{\delta}$ is non-empty, so
$P^{\operatorname{ext},\delta}_{A_{0,\delta}}\geq 1$. On the other hand, for
any $x,x^{\prime}\in B_{\delta},\left\lVert x-x^{\prime}\right\rVert<\delta$,
hence they cannot be part of the same $\delta$-packing. ∎

######  Proposition 19.

Let $X=\mathbb{R}^{d}$, and let $k$ satisfy spatial decay. Then there is a
constant $C_{\max}^{k,\delta}$ such that for any $N$ and any
$\delta$-separated $\boldsymbol{x}$ of size $N$ we have

$\lambda_{\max}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}})\leq
C_{\max}^{k,\delta}$ (44)

###### Proof.

Let $\psi$ be the upper-bounding function. For an arbitrary $x_{i}$, we have

$\lambda_{\max}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}})\leq
k_{ii}+R_{i}\leq\sum_{j=1}^{N}\psi(\left\lVert x_{i}-x_{j}\right\rVert)$ (45)

so it suffices to bound this sum. For that, we have

$\sum_{j=1}^{N}\psi(\left\lVert
x_{i}-x_{j}\right\rVert)=\sum_{m=0}^{\infty}\sum_{x_{j}\in
A_{m,\delta}}\psi(\left\lVert
x_{i}-x_{j}\right\rVert)\leq\sum_{m=0}^{\infty}P^{\operatorname{ext},\delta}_{A_{m,\delta}}\psi(m\delta).$
(46)

The first equality uses that $\\{A_{m,\delta}\\}_{m=0}^{\infty}$ forms a
partition of $\mathbb{R}^{d}$. The second uses that for $x_{j}\in
A_{m,\delta}$, we have $\psi(\left\lVert
x_{i}-x_{j}\right\rVert)\leq\psi(m\delta)$, and the maximum number of $x_{j}$
in each $A_{m,\delta}$ is bounded by the external packing number
$P^{\operatorname{ext},\delta}_{A_{m,\delta}}$ by the definition of a
$\delta$-separated set. By the preceding lemma, spatial decay, and the
integral test for infinite series, we conclude

$\sum_{m=0}^{\infty}P^{\operatorname{ext},\delta}_{A_{m,\delta}}\psi(m\delta)\leq\psi(0)+2d\left(\frac{2}{\delta}\right)^{d}\sum_{m=1}^{\infty}\left(m-\frac{1}{2}\right)^{d-1}\psi(m\delta)=C_{\max}^{k,\delta}<\infty.$
(47)

∎

###  A.4 Lower Eigenvalue Estimates

For lower bounds, we apply (wendland04, 56, Theorem 12.3), which in the
stationary case follows from our assumptions. This argument, which builds on
the work of narcowich92 (37, 42), is somewhat involved. To aid understanding,
rather than simply quote the result, below we present a sketch of its key
steps.

The idea is to invoke the lower estimate given by Gershgorin’s Circle
Theorem—Result 13—which was previously used for the upper bounds. However,
applying this estimate directly leads in general to the inequality
$\lambda_{\min}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}})\geq\widetilde{C}_{\min}^{k,\delta}$
where potentially $\widetilde{C}_{\min}^{k,\delta}\leq 0$, which is vacuous.
To overcome this, the idea is to carefully construct a different covariance
kernel $\phi$ whose eigenvalues lower-bound those of $k$, and which admit
better Gershgorin estimates. This is done via the following lemma.

######  Lemma 20.

Let $k$ and $\phi$ be stationary kernels with respective spectral densities
$\rho$ and $\varsigma$. Suppose that $\varsigma\leq\rho$. Then for all
$\boldsymbol{x}$ we have

$\lambda_{\min}(\mathbf{\Phi}_{\boldsymbol{x}\boldsymbol{x}})\leq\lambda_{\min}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}).$
(48)

###### Proof.

For any $\boldsymbol{v}\in\mathbb{R}^{N}$ and by Bochner’s Theorem—see
(wendland04, 56, Theorem 6.6)—we have

$\displaystyle\boldsymbol{v}^{T}\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}\boldsymbol{v}$
$\displaystyle=\sum_{i=1}^{N}\sum_{j=1}^{N}v_{i}k(0,x_{i}-x_{j})v_{j}$ (49)
$\displaystyle=\sum_{i=1}^{N}\sum_{j=1}^{N}v_{i}v_{j}\int_{\mathbb{R}^{d}}e^{2\pi
i\langle\omega,x_{i}-x_{j}\rangle}\mathop{}\\!\mathrm{d}\rho(\omega)$ (50)
$\displaystyle=\int_{\mathbb{R}^{d}}\sum_{i=1}^{N}v_{i}e^{2\pi
i\langle\omega,x_{i}\rangle}\sum_{j=1}^{N}v_{j}e^{2\pi
i\langle\omega,x_{j}\rangle}\mathop{}\\!\mathrm{d}\rho(\omega)$ (51)
$\displaystyle=\int_{\mathbb{R}^{d}}\Bigg{|}\sum_{j=1}^{N}v_{j}e^{2\pi i\omega
x_{j}}\Bigg{|}^{2}\mathop{}\\!\mathrm{d}\rho(\omega)$ (52)
$\displaystyle\geq\int_{\mathbb{R}^{d}}\Bigg{|}\sum_{j=1}^{N}v_{j}e^{2\pi
i\omega x_{j}}\Bigg{|}^{2}\mathop{}\\!\mathrm{d}\varsigma(\omega)$ (53)
$\displaystyle=\boldsymbol{v}^{T}\mathbf{\Phi}_{\boldsymbol{x}\boldsymbol{x}}\boldsymbol{v}$
(54)

where, in the third inequality, we use that since $v_{j}$ is real, it is equal
to its complex conjugate. Taking an infimum over $\boldsymbol{v}$ on the unit
sphere on the left hand side gives the result. ∎

With this established, we can obtain lower eigenvalue estimates for any kernel
$\phi$ whose Gershgorin radius is small enough, by invoking separation and
using a packing argument similar to the one used to prove Proposition 19. This
is done as follows.

######  Lemma 21.

Let $X=\mathbb{R}^{d}$, and let $\phi$ decay sufficiently fast. Then there is
a constant $C^{\phi,\delta}_{\min}>0$ such that for any $N$ any any
$\delta$-separated $\boldsymbol{x}$ of size $N$ we have

$\lambda_{\min}(\mathbf{\Phi}_{\boldsymbol{x}\boldsymbol{x}})\geq
C^{\phi,\delta}_{\min}.$ (55)

###### Proof.

For arbitrary $x_{i}$, letting $R_{i}$ be the Gershgorin radius, we have

$\lambda_{\min}(\mathbf{\Phi}_{\boldsymbol{x}\boldsymbol{x}})\geq\Phi_{ii}-R_{i}$
(56)

so it suffices to upper-bound $R_{i}$ sufficiently tightly. Applying the
packing argument of Proposition 4, we get

$\sum_{j\neq i}\phi(x_{i}-x_{j})\leq
d\left(\frac{3}{\delta}\right)^{d}\sum_{m=1}^{\infty}m^{d-1}\phi(m\delta)$
(57)

where we have substituted the lower-bounding kernel $\phi$ in to the argument,
modified slightly to exclude the diagonal part, which corresponds to $m=0$.
The claim follows, provided

$d\left(\frac{3}{\delta}\right)^{d}\sum_{m=1}^{\infty}m^{d-1}\phi(m\delta)<\phi(0)$
(58)

which holds provided that $\phi$ decays sufficiently fast. ∎

From here, all that remains is to show that for any stationary kernel $k$, it
is always possible to choose a $\phi$ that decays sufficiently fast.
wendland04 (56) shows that for every kernel $k$ there is a constant
$0<M_{k}<\infty$ for which taking the spectral density $\varsigma$ of $\phi$
to be the convolution of an indicator function of the ball $B(0,M_{k})$ with
itself, namely

$\varsigma\propto\mathbbold{1}_{B(0,M_{k})}*\mathbbold{1}_{B(0,M_{k})}$ (59)

with proportionality scaling chosen to ensure the spectral lower-bound needed
for Lemma 20 holds, suffices to yield a non-vacuous bound. The idea behind
this choice, which dates back at least to ideas of narcowich92 (37, 42) in the
polynomial interpolation literature, is to leverage the uncertainty principle
property of the Fourier transform. Specifically, any kernel concentrated
enough to decay sufficiently fast in the spatial domain must possess a
sufficiently diffuse spectral measure. The simplest choice is to take
$\varsigma\propto\mathbbold{1}_{B(0,M_{k})}$, but this fails due to issues
with smoothness, which are fixed by applying a convolution and rescaling
appropriately, yielding the presented choice.

At this point, the main technical work that remains is to calculate the form
of the resulting kernel $\phi$, and deduce how to choose the constant $M_{k}$
and proportionality scaling to ensure it both lower-bounds $k$ and gives a
non-vacuous estimate. For that, we refer to (wendland04, 56, Ch. 12.2).

######  Result 22.

Let $X=\mathbb{R}^{d}$, and let $k$ be stationary. Then there is a constant
$C_{\min}^{k,\delta}$ such that for any $N$ and any $\delta$-separated
$\boldsymbol{x}$ of size $N$ we have

$\lambda_{\min}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}})\geq
C_{\min}^{k,\delta}.$ (60)

###### Proof.

(wendland04, 56, Theorem 12.3). ∎

###  A.5 Settings where Separation is a Necessary Condition

Let $k$ be a kernel. We say that $k$ is a _Lipschitz continuous kernel_ if $k$
is Lipschitz in both arguments separately. We extend the notion of minimum
separation from finite vectors to infinite sequences in the obvious pairwise
manner.

*

###### Proof.

We argue by contradiction: suppose minimum separation does not hold. We will
show that there does not exist a $C$ such that, for all $M$,
$\operatorname{cond}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}^{(M)})\leq C$.
By Cauchy’s Interlacing Theorem, any submatrix of
$\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}^{(M)}$ has condition number no
larger than
$\operatorname{cond}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}^{(M)})$. Hence,
it suffices to find a sequence of submatrices $\mathbf{K}_{m}$, for
$m=1,..,\infty$, such that the rows and columns of $\mathbf{K}_{m}$ are
contained in the rows and columns of
$\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}^{(m)}$, and show that this sequence
has unbounded condition number. We therefore consider such sequences,
consisting of $2\times 2$ matrices, in the remainder of the proof.

We now split into two cases. The first case handles when the variance of the
kernel can differ by an arbitrary factor at different points in input space.
This is handled as a separate case, as the technique in the second, arguably
more intuitive case, will break down if the variance can be arbitrarily close
to zero.

_Case I: for any $C>0$, there exists an $i,j$ with
$\frac{k(z_{i},z_{i})}{k(z_{j},z_{j})}>C$_. Let $\mathbf{K}_{m}$ be the matrix
formed by

$\displaystyle z^{(m)}_{i}$
$\displaystyle=\operatorname*{\arg\max}_{i^{\prime}\leq
m}k(z_{i^{\prime}},z_{i^{\prime}})$ $\displaystyle z^{(m)}_{j}$
$\displaystyle=\operatorname*{\arg\min}_{j^{\prime}\leq
m}k(z_{j^{\prime}},z_{j^{\prime}})$ (61)

where we take take the smallest such $i,j$ if the maximum and minimum are not
unique. Then

$\displaystyle\operatorname{cond}(\mathbf{K}_{m})\geq\operatorname{cond}\begin{pmatrix}k(z^{(m)}_{i},z^{(m)}_{i})&0\\\
0&k(z^{(m)}_{j},z^{(m)}_{j})\end{pmatrix}=\frac{k(z^{(m)}_{i},z^{(m)}_{i})}{k(z^{(m)}_{j},z^{(m)}_{j})}.$
(62)

The inequality can be verified by inspection using the formula for eigenvalues
of a $2\times 2$ matrix. By the assumption in this case, this ratio tends to
infinity with $m$, contradicting that the condition number is bounded.

_Case II: there exists a $C>0$ such that for all pairs $i,j$ we have
$\frac{k(z_{i},z_{i})}{k(z_{j},z_{j})}\leq C$._ By the assumption in this
case, for all $i$, we have

$0<\frac{1}{C}k(z_{1},z_{1})\leq k(z_{i},z_{i})\leq Ck(z_{1},z_{1}).$ (63)

Define $a=\frac{1}{C}k(z_{1},z_{1})$. Take $z^{(m)}_{i},z^{(m)}_{j}$ to be any
pair such that $\lVert z^{(m)}_{i}-z^{(m)}_{j}\rVert\leq\left\lVert
z_{i^{\prime}}-z_{j^{\prime}}\right\rVert$ for any pair
$i^{\prime},j^{\prime}\leq m$, and $\mathbf{K}_{m}$ to be the submatrix
associated to this pair of points. By Lipschitzness of the kernel, there
exists an $L$ such that

$\displaystyle|k(z^{(m)}_{i},z^{(m)}_{i})-k(z^{(m)}_{i},z^{(m)}_{j})|$
$\displaystyle\leq L\operatorname{sep}_{m}$
$\displaystyle|k(z^{(m)}_{i},z^{(m)}_{i})-k(z^{(m)}_{k},z^{(m)}_{j})|$
$\displaystyle\leq 2L\operatorname{sep}_{m}.$ (64)

where $\operatorname{sep}_{m}$ denotes the separation distance of
$(z_{1},..,z_{m})$. Combining these estimates gives the Frobenius norm bound

$\left\lVert\frac{1}{k(z_{i},z_{i})}\mathbf{K}_{m}-\mathbf{1}\right\rVert_{F}\leq\frac{L\operatorname{sep}_{m}\sqrt{6}}{k(z_{i},z_{i})}\leq\frac{L\operatorname{sep}_{m}\sqrt{6}}{a},$
(65)

where $\mathbf{1}$ is the $2\times 2$ matrix of all ones. As the condition
number of $\frac{1}{k(z_{i},z_{i})}\mathbf{K}_{m}$ and $\mathbf{K}_{m}$ are
the same, we study the the condition number of
$\frac{1}{k(z_{i},z_{i})}\mathbf{K}_{m}$. The operator norm is less than the
Frobenius norm, so by reverse triangle inequality applied to the operator norm

$\lambda_{1}(\frac{1}{k(z_{i},z_{i})}\mathbf{K}_{m})\geq
2-\left\lVert\frac{1}{k(z_{i},z_{i})}\mathbf{K}_{m}-\mathbf{1}\right\rVert_{2}\geq
2-\frac{L\operatorname{sep}_{m}\sqrt{6}}{a}$ (66)

and

$\lambda_{2}(\frac{1}{k(z_{i},z_{i})}\mathbf{K}_{m})\leq\frac{L\operatorname{sep}_{m}\sqrt{6}}{a}.$
(67)

By our assumption that $\operatorname{sep}_{m}$ is not bounded below, this
ratio of these tends to infinity with $m$. ∎

###  A.6 Effect of Dimension on Condition Number

We now consider the simple case of independent inputs in $\mathbb{R}^{d}$ with
independent and identically distributed dimensions. The key observation is
that for large $d$, by concentration of measure, we should most pairwise
distances to concentrate around the average distance. This means that there
will not be many tightly clustered sets of data points. Instead, the data will
be automatically well-separated with high probability in this specific
setting.

Let $\rho$ be the distribution of each input. We work with uniformly generated
inputs for Figure 4, specifically $\rho=\operatorname{U}(-2\pi,2\pi)$. For any
$i\neq j$ we have

$\left\lVert
x_{i}-x_{j}\right\rVert^{2}_{2}=\sum_{\ell=1}^{d}\Delta_{\ell}^{2},$ (68)

with $\Delta_{\ell}=(\xi_{i,\ell}-\xi_{j,\ell})^{2}$ where
$\xi_{i,\ell}\sim\rho$ are all independent, hence $\Delta_{\ell}$ are also all
independent.

If we consider the isotropic squared exponential kernel with length scale
$\sqrt{d}$, which is used for generating Figure 4, the off-diagonal entries of
the kernel matrix are independently and identically distributed random
variables of the form

$k_{ij}=\exp\left(-\frac{1}{2d}\sum_{\ell=1}^{d}\Delta_{\ell}^{2}\right).$
(69)

By the Law of Large Numbers, as $d\rightarrow\infty$, we have

$\frac{1}{2d}\sum_{\ell=1}^{d}\Delta_{\ell}^{2}\rightarrow\frac{1}{2}\operatorname{\mathbb{E}}(\Delta_{\ell}^{2}),$
(70)

and so by the Continuous Mapping Theorem we get

$k_{ij}\rightarrow\exp\left(-\frac{1}{2}\operatorname{\mathbb{E}}(\Delta_{\ell}^{2})\right)$
(71)

for large $d$. On the other hand, the diagonal entries of the kernel matrix
are all surely $1$.

In the case of Figure 4, we have
$(x_{i})_{\ell}\sim\operatorname{U}(-2\pi,2\pi)$. Hence,

$\operatorname{\mathbb{E}}(\Delta_{\ell}^{2})=\frac{4}{3}\pi^{2},$ (72)

so for $d$ large, we have
$k_{ij}\rightarrow\exp\left(-\frac{2\pi^{2}}{3}\right)$. For fixed $N$ and
large $d$, we have the approximation

$\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}\rightarrow\left(1-\exp\left(-\frac{2\pi^{2}}{3}\right)\right)\mathbf{I}+\exp\left(-\frac{2\pi^{2}}{3}\right)\boldsymbol{1}\boldsymbol{1}^{T}.$
(73)

The spectrum of the matrix on the right hand side is

$\lambda(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}})=\left(1+(N-1)\exp(-\tfrac{2\pi^{2}}{3}),1-\exp(-\tfrac{2\pi^{2}}{3}),\dots,1-\exp(-\tfrac{2\pi^{2}}{3})\right).$
(74)

Since, in our experiments, $\exp(-\tfrac{2\pi^{2}}{3})\approx\frac{1}{N}$, the
condition number of $\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}$ will be small.
Since the entries of $\mathbf{\Lambda}$ are small, this implies the condition
number of $\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}+\mathbf{\Lambda}$ will be
small. While this argument only holds asymptotically one might expect the
central limit theorem to apply relatively quickly, since the inputs are light-
tailed, thus explaining observed behavior for, say, $d=8$ in Figure 4.

##  B Inducing Points

Here we develop the inducing-point-specific theory used in the manuscript.

###  B.1 Clustering-based Inducing Points

*

###### Proof.

The likelihood can be written

$p(y_{i}\mid
f(x_{i}))\propto\exp\left(-\sum_{i=1}^{N}\frac{(y_{i}-f(\operatorname{cl}(x_{i})))^{2}}{2\sigma^{2}}\right)=\exp\left(-\sum_{j=1}^{M}\sum_{\operatorname{cl}(x_{i})=z_{j}}\frac{(y_{i}-f(z_{j}))^{2}}{2\sigma^{2}}\right)$
(75)

and the claim follows by noting that by Gaussianity, each
$\left(u_{j},\frac{1}{N_{\operatorname{cl}}(z_{j})}\right)$ pair is a
sufficient statistic for all $y_{i}$ whose indices satisfy
$\operatorname{cl}(x_{i})=z_{j}$. ∎

###  B.2 Stochastic Maximum Marginal Likelihood

We train hyperparameters using maximum marginal likelihood via doubly
stochastic gradient descent. The minimization problem resembles the usual
variational inference problem in most inducing point approaches, and is

$\operatorname*{\arg\min}_{q_{f}\in\mathbb{Q}}D_{\operatorname{KL}}(q_{f}\mid\mid\pi_{f})-\frac{N}{2}\log(2\pi\sigma^{2})-\frac{1}{2\sigma^{2}}\operatorname*{\mathbb{E}}_{f\sim
q_{f}}(\boldsymbol{y}-f(\boldsymbol{x}))^{T}(\boldsymbol{y}-f(\boldsymbol{x}))$
(76)

where $\mathbb{Q}$ is the space of all Gaussian processes of the form

$(f\mid\boldsymbol{u})(\cdot)=f(\cdot)+\mathbf{K}_{(\cdot)\boldsymbol{z}}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})^{-1}(\boldsymbol{u}-f(\boldsymbol{z})-\boldsymbol{\varepsilon})$
(77)

and $\boldsymbol{z}=\operatorname{cl}(\boldsymbol{x})$,
$u_{j}=\frac{1}{N_{\operatorname{cl}}(z_{j})}\sum_{\operatorname{cl}(x_{i})=z_{j}}y_{i}$,
and $\Lambda_{ii}=\frac{\sigma^{2}}{N_{\operatorname{cl}}(z_{i})}$. The only
parameters which are minimized numerically are the kernel hyperparameters. We
simplify the variational expectation term using the identity

$\operatorname{\mathbb{E}}(\boldsymbol{y}-f(\boldsymbol{x}))^{T}(\boldsymbol{y}-f(\boldsymbol{x}))=\sum_{i=1}^{N}\operatorname{\mathbb{E}}(y_{i}-f(x_{i}))^{2}=\sum_{i=1}^{N}(y_{i}-\operatorname{\mathbb{E}}(f(x_{i})))^{2}+\operatorname{Var}(f(x_{i})).$
(78)

To train the process, we need to do two things: compute the prior
Kullback–Leibler divergence, and compute the likelihood term. For the latter,
we sample a mini-batch of data, and use it compute an unbiased estimator of
the sum by evaluating the inner terms on the mini-batch. For a mini-batch
$\boldsymbol{x}^{\prime}$, the mean and variance terms are

$\displaystyle\operatorname{\mathbb{E}}(f(\boldsymbol{x}^{\prime}))$
$\displaystyle=\mathbf{K}_{\boldsymbol{x}^{\prime}\boldsymbol{z}}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})^{-1}\boldsymbol{u}$
$\displaystyle\operatorname{Var}(f(\boldsymbol{x}^{\prime}))$
$\displaystyle=\mathbf{K}_{\boldsymbol{x}^{\prime}\boldsymbol{z}}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})^{-1}\mathbf{K}_{\boldsymbol{z}\boldsymbol{x}^{\prime}}$
(79)

which can be computed using matrix-vector products whose cost, in the case of
the variance, is quadratic for each data point in the mini-batch. For the
prior term, we have the identity

$\displaystyle 2D_{\operatorname{KL}}(q_{f}\mid\mid\pi_{f})$
$\displaystyle=\ln\frac{|\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}|}{|\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}-\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})^{-1}\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}|}-d$
(80)
$\displaystyle\quad+\operatorname{tr}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}^{-1}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}-\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})^{-1}\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}))$
(81)
$\displaystyle\quad+\boldsymbol{u}^{T}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})^{-1}\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}^{-1}\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})^{-1}\boldsymbol{u}$
(82)
$\displaystyle=\ln\frac{|\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda}|}{|\mathbf{\Lambda}|}-\operatorname{tr}((\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})^{-1}\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}})$
(83)
$\displaystyle\quad+\boldsymbol{u}^{T}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})^{-1}\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})^{-1}\boldsymbol{u}.$
(84)

This means that all linear systems that need to be solved involve the well-
conditioned matrix
$(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})^{-1}$.
Unfortunately,
$(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})^{-1}\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}$
has a matrix-valued rather than vector-valued right-hand-side, and thus
requires $M$ total matrix-vector products, leading to cubic costs, even with
conjugate gradients. To avoid this, we apply Hutchinson’s trace estimator,
which we now describe.

To illustrate Hutchinson’s trace estimator, let $\mathbf{A}$ an arbitrary
matrix, for instance
$\mathbf{A}=(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})^{-1}\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}$.
We are interested in evaluating $\operatorname{tr}(\mathbf{A})$ in cases where
we do not have direct access to $\mathbf{A}$, but can compute matrix-vector
products with $\mathbf{A}$. hutchinson1989stochastic (23) observed that for
any vector $\boldsymbol{v}$ such that
$\operatorname{\mathbb{E}}(\boldsymbol{v}\boldsymbol{v}^{T})=\mathbf{I}$ we
have

$\operatorname{tr}(\mathbf{A})=\operatorname{tr}(\operatorname{\mathbb{E}}(\boldsymbol{v}\boldsymbol{v}^{T})\mathbf{A})=\operatorname{\mathbb{E}}(\boldsymbol{v}^{T}\mathbf{A}\boldsymbol{v}).$
(85)

This motivates one to employ the unbiased Monte Carlo estimator

$\frac{1}{L}\sum_{\ell=1}^{L}\boldsymbol{v}_{\ell}^{T}\mathbf{A}\boldsymbol{v}_{\ell}$
(86)

where the $\boldsymbol{v}_{\ell}$ are independent and identically distributed.
Typically, $\boldsymbol{v}_{\ell}$ is chosen to be Rademacher distributed,
where each entry of $\boldsymbol{v}$ is independent and takes the values
$\\{-1,1\\}$ with equal probability. This choice of $\boldsymbol{v}$ is leads
to the minimum variance—see (hutchinson1989stochastic, 23, Proposition 1).

In the case
$\mathbf{A}=(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})^{-1}\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}$,
we sample the Rademacher random variables, solve the linear system
$(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})^{-1}\boldsymbol{v}$,
then multiply the result by
$\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}\boldsymbol{v}$. The complexity of
the matrix-vector products required for both operations is quadratic.

The log-determinant term is handled similarly. Rather than evaluate it
directly, we follow gibbs1997efficient (18) and apply stochastic estimation on
its gradient. Using that

$\frac{\partial{|\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda}|}}{\partial{\theta}}=\operatorname{tr}\left(\frac{\partial{(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})}}{\partial{\theta}}(\mathbf{K}_{\boldsymbol{z}\boldsymbol{z}}+\mathbf{\Lambda})^{-1}\right)$
(87)

we apply Hutchinson’s trace estimator as in the case of the trace discussed
above.

##  C Cover Trees

To prove the key claim, we formalize the cover tree as a labeled rooted tree
in the sense of a directed graph, defined inductively by the steps of the
algorithm. We distinguish individual nodes from their locations in space,
which we view as labels.

######  Definition 23.

Let $T=(V,E,v_{0})$ be a rooted tree, with root node $v_{0}\in V$. Define the
_tree order_ $\prec$ by $v\prec v^{\prime}$ if the unique path from $v_{0}$ to
$v^{\prime}$ passes through $v$. Define the _parent_ of a node to be

$\operatorname{pa}:V\setminus\\{v_{0}\\}\rightarrow V$ (88)

which maps each $v$ to the unique vertex $\operatorname{pa}(v)\prec v$ for
which $(v,\operatorname{pa}(v))\in E$, which is well-defined because $T$ is a
rooted tree. Define the _children_ of a node to be

$\displaystyle\operatorname{ch}$ $\displaystyle:V\rightarrow 2^{V}$
$\displaystyle\operatorname{ch}(v)$
$\displaystyle=\operatorname{pa}^{-1}(\\{v\\})$ (89)

which is the set-theoretic preimage of the parent function. Define the _level_

$\operatorname{lv}:V\rightarrow\mathbb{N}_{0}$ (90)

of a node to be function which maps each $v$ to the length
$\operatorname{lv}(v)$ of the path from $v_{0}$ to $v$, which is also well-
defined because $T$ is a rooted tree.

With these notions, we can define an intermediate state of the cover tree
algorithm.

######  Definition 24.

We say that a _cover tree_ is a tuple $(T,z,\mathcal{A})$, where $T$ is a
rooted tree and

$\displaystyle z$ $\displaystyle:V\rightarrow\boldsymbol{x}$
$\displaystyle\mathcal{A}$ $\displaystyle:V\rightarrow 2^{\boldsymbol{x}}$
(91)

are, respectively, the _spatial location_ and _assigned data_ of each node.

The cover tree algorithm, then, is a function which maps a cover tree into
another cover tree, by constructing the next level. Our goal is to prove that
this algorithm terminates in finite time, producing a rooted tree whose nodes
satisfy the desired properties. We analyze a more general form of the
algorithm parameterized by a constant $a>1$ and sequence $b_{\ell}>0$, where
$a$ determines the rate at which the radius of each ball decays, and
$b_{\ell}$ determines the size of an $R$-neighbor ball at level $\ell$, with
Algorithm 9 corresponding to $a=2$ and $b_{\ell}=4(1-1/2^{\ell})R_{0}$, which
will be derived as the optimal sequence for the given choice of $a$.

######  Definition 25.

Let $V_{0}=\\{v_{0}\\}$, $E_{0}=\varnothing$,
$z_{0}(v_{0})=\frac{1}{N}\sum_{i=1}^{n}x_{i}$, and
$\mathcal{A}_{0}=\boldsymbol{x}$. Let $(T_{\ell},z_{\ell},\mathcal{A}_{\ell})$
be the output of one step of Algorithm 9, defined inductively starting from
$\ell=0$.

From this, we can define the inducing points produced by the algorithm.

######  Definition 26.

For a given tree $T$, define the _level- $\ell$ nodes_ and _level- $\ell$
inducing points_

$\displaystyle V_{\ell}(T)$ $\displaystyle=\\{v\in
V:\operatorname{lv}(v)=\ell\\}$ $\displaystyle\boldsymbol{z}_{\ell}(T)$
$\displaystyle=\\{z(v):v\in V:\operatorname{lv}(v)=\ell\\}.$ (92)

We first verify that the nodes and hence inducing points at a fixed level
remain fixed as the algorithm proceeds.

######  Lemma 27.

For $\ell^{\prime}\geq\ell$, we have

$\displaystyle V_{\ell}(T_{\ell^{\prime}})$ $\displaystyle=V_{\ell}(T_{\ell})$
$\displaystyle\operatorname{sep}(\boldsymbol{z}_{\ell}(T_{\ell^{\prime}}))$
$\displaystyle=\operatorname{sep}(\boldsymbol{z}_{\ell}(T_{\ell}))$
$\displaystyle\operatorname{res}_{\boldsymbol{x}}(\boldsymbol{z}_{\ell}(T_{\ell^{\prime}}))$
$\displaystyle=\operatorname{res}_{\boldsymbol{x}}(\boldsymbol{z}_{\ell}(T_{\ell})).$
(93)

###### Proof.

The second and third claims follow from the first one, which is immediate by
the construction of Algorithm 9. ∎

Next, we check that the separation and spatial resolution properties hold at
the root node.

######  Lemma 28.

We have

$\displaystyle\operatorname{sep}(\boldsymbol{z}_{0}(T_{0}))$
$\displaystyle=\infty$
$\displaystyle\operatorname{res}_{\boldsymbol{x}}(\boldsymbol{z}_{0}(T_{0}))\leq\max_{x\in\boldsymbol{x}}\left\lVert
x-z(v_{0})\right\rVert.$ (94)

###### Proof.

Because $V_{0}=\\{v_{0}\\}$, the supremum in the definition of
$\operatorname{sep}$ is taken over an empty set, and the first claim follows.
For the second claim, recall that
$z(v_{0})=\frac{1}{\boldsymbol{x}}\sum_{x^{\prime}\in\boldsymbol{x}}x^{\prime}$.
For any $x\in\boldsymbol{x}$,

$\big{\lVert}x-z(v_{0})\big{\rVert}=\frac{1}{|\boldsymbol{x}|}\Big{\lVert}\sum_{x^{\prime}\in\boldsymbol{x}}x-x^{\prime}\Big{\rVert}\leq\frac{1}{|\boldsymbol{x}|}\sum_{x^{\prime}\in\boldsymbol{x}}\left\lVert
x-x^{\prime}\right\rVert\leq\max_{x^{\prime}\in\boldsymbol{x}}\left\lVert
x-x^{\prime}\right\rVert.$ (95)

Taking a maximum over $x\in\boldsymbol{x}$ completes the proof. ∎

We now define the notion of $R$-neighbors.

######  Definition 29.

Define $R_{\ell}=\frac{R_{\ell-1}}{a}$ inductively, where
$R_{0}=\max_{x\in\boldsymbol{x}}\left\lVert x-z(v_{0})\right\rVert$. For
$b>0$, define the _$R_{\ell}$ -neighbors_ of a node inductively by

$\displaystyle\mathcal{R}(v_{0})$ $\displaystyle=\\{v_{0}\\}$
$\displaystyle\mathcal{R}(v)$
$\displaystyle=\left\\{v^{\prime}\in\bigcup_{s\in\mathcal{R}(\operatorname{pa}(v))}\operatorname{ch}(s):\left\lVert
z(v)-z(v^{\prime})\right\rVert\leq b_{\ell}R_{\ell}\right\\}$ (96)

The next step is to deduce what sequences $b_{\ell}$ are valid for a given
value of $a$.

######  Lemma 30.

Let

$b_{\ell}=\frac{2(1+c(\frac{1}{a})^{\ell})}{1-\frac{1}{a}}.$ (97)

for $c\in\\{-1,-a^{L}\\}$. Then, letting

$\mathcal{R}^{\prime}(v)=\left\\{v^{\prime}\in V_{\ell}(T_{\ell}):\left\lVert
z(v)-z(v^{\prime})\right\rVert\leq b_{\ell}R_{\ell}\right\\}$ (98)

be the set of $R_{\ell}$-nearby nodes, we have

$\mathcal{R}^{\prime}(v)=\mathcal{R}(v).$ (99)

###### Proof.

First, note for all $v\in V_{\ell}(T_{\ell})$ that
$\mathcal{R}(v)\subseteq\mathcal{R}^{\prime}(v)$, since

$\bigcup_{s\in\mathcal{R}(\operatorname{pa}(v))}\operatorname{ch}(s)\subseteq
V_{\ell}(T_{\ell})$ (100)

since the set of children of the parents of a node are contained in the same
level as the original node. It therefore suffices to show
$\mathcal{R}^{\prime}(v)\subseteq\mathcal{R}(v)$. The proof of this proceeds
by induction on depth. In the base case $\ell=0$ there is a single node
$v_{0}$, and $v_{0}\in\mathcal{R}(v_{0})=\\{v_{0}\\}$, so the proposition is
true. We take the inductive hypothesis that

$\mathcal{R}^{\prime}(v^{\prime\prime})\subseteq\mathcal{R}(v^{\prime\prime})$
(101)

for all $v^{\prime\prime}\in V_{\ell-1}(T_{\ell-1})=V_{\ell-1}(T_{\ell})$,
where the equality follows from Lemma 27note . We now show the claim holds for
$v\in V_{\ell}(T_{\ell})$. Let

$v^{\prime}\in\mathcal{R}^{\prime}(v)$ (102)

for which we would like to show that $v^{\prime}\in\mathcal{R}(v)$. The
strategy will be to consider these nodes’ parents, and show that
$\operatorname{pa}(v^{\prime})\in\mathcal{R}(\operatorname{pa}(v))$, in which
case the result follows from the definition of $\mathcal{R}$. By the
combination of the inductive hypothesis and first part of the proof, we have
$\mathcal{R}^{\prime}(v^{\prime\prime})=\mathcal{R}(v^{\prime\prime})$ for all
$v^{\prime\prime}\in V_{\ell-1}$, so it suffices to show

$\left\lVert
z(\operatorname{pa}(v^{\prime}))-z(\operatorname{pa}(v))\right\rVert\leq
b_{\ell-1}R_{\ell-1}.$ (103)

We therefore upper bound the distance between the locations of the parents of
the nodes in question. We have

$\displaystyle\left\lVert
z(\operatorname{pa}(v^{\prime}))-z(\operatorname{pa}(v))\right\rVert$
$\displaystyle\leq\left\lVert
z(\operatorname{pa}(v^{\prime}))-z(v^{\prime})\right\rVert+\left\lVert
z(\operatorname{pa}(v))-z(v)\right\rVert+\left\lVert
z(v^{\prime})-z(v)\right\rVert$ (104)
$\displaystyle\leq\|z(\operatorname{pa}(v^{\prime}))-z(v^{\prime})\|+\|z(\operatorname{pa}(v))-z(v)\|+b_{\ell}R_{\ell}$
(105) $\displaystyle\leq 2R_{\ell-1}+b_{\ell}R_{\ell}.$ (106)
$\displaystyle=2R_{\ell-1}+\frac{b_{\ell}}{a}R_{\ell-1}.$ (107)

The first inequality by the triangle inequality. The second inequality follows
by definition of $v^{\prime}$. The third inequality uses that the locations of
the children of a node are contained in the assigned data of the node, and
that all points contained in the assigned data of a node in level $\ell-1$
must be within a distance of $R_{\ell-1}$ of the label, or in equations that
for all $v^{\prime\prime}\in V_{\ell-1}(T_{\ell})$ and
$x\in\mathcal{A}_{\ell-1}(v)$ we have

$\left\lVert z(v^{\prime\prime})-x\right\rVert\leq R_{\ell-1}.$ (108)

Combining (107) and (103) gives non-homogenous linear recurrence relation

$2+\frac{1}{a}b_{\ell}\leq b_{\ell-1}.$ (109)

The solution can be derived using standard theory of recurrence relations
(greene1990mathematics, 20, Chapter 2), which gives

$b_{\ell}=\frac{2(1+c(\frac{1}{a})^{\ell})}{1-\frac{1}{a}}$ (110)

for some constant $c$ to be determined by the initial conditions—one can
easily check by plugging in the formula that this satisfies the recurrence. We
require that $b_{\ell}\geq 0$ for $1\leq\ell\leq L$, where $L$ is the maximum
tree depth. The optimal choice of initial conditions given this constraint, in
the sense of minimizing the size of neighborhoods while retaining guarantees,
come from selecting $b_{0}=0$ or $b_{L}=0$. The former gives $c=-1$, while the
latter gives $c=-a^{L}$. This completes the proof. ∎

We now argue that the spatial resolution of the algorithm satisfies the
claimed inequality.

######  Lemma 31.

We have

$\operatorname{res}_{\boldsymbol{x}}(\boldsymbol{z}_{\ell}(T_{\ell}))\leq
R_{\ell}.$ (111)

###### Proof.

By induction on depth, one can show that for every level,
$\\{\mathcal{A}_{\ell}(v)\\}_{v\in V_{\ell}}$ is a partition of
$\bigcup_{i=1}^{N}\\{x_{i}\\}$, where $x_{i}$ are entries of $\boldsymbol{x}$.
For $x\in X$, let $v_{x,\ell}$ denote the unique $v\in V_{\ell}$ such that
$x\in\mathcal{A}_{\ell}(v_{x,\ell})$. Then

$\min_{v\in V_{\ell}}\|x-z(v)\|\leq\|x-z(v_{x,\ell})\|\leq R_{\ell},$ (112)

where the second inequality uses that the active set is contained in the ball
of radius $R_{\ell}$ centered at $z(v_{x,\ell})$, which holds by definition of
Algorithm 9, specifically 13. If optional Voronoi repartitioning is performed,
we further upper bound $\left\lVert x-z(v_{x,\ell})\right\rVert$ by $\lVert
x-z(v^{\prime}_{x,\ell})\rVert$, where $v^{\prime}_{x,\ell}$ denotes the
vertex to which $x$ was assigned prior to the Voronoi step. The results
follows from taking a maximum over $x\in\boldsymbol{x}$. ∎

Next, we argue the separation distance inequality holds.

######  Lemma 32.

We have

$\operatorname{sep}(\boldsymbol{z}_{\ell}(T_{\ell}))\geq R_{\ell}.$ (113)

###### Proof.

We prove the contrapositive: suppose that there exist a $v,v^{\prime}\in
V_{\ell}$ such that $\left\lVert z(v)-z(v^{\prime})\right\rVert\leq R_{\ell}$.
By the reverse triangle inequality, this yields

$\displaystyle R_{\ell}$ $\displaystyle\geq\left\lVert
z(v)-z(v^{\prime})\right\rVert$ (114) $\displaystyle\geq\left\lVert
z(\operatorname{pa}(v))-z(\operatorname{pa}(v^{\prime}))\right\rVert-\left\lVert
z(\operatorname{pa}(v))-z(v)\right\rVert-\left\lVert
z(\operatorname{pa}(v^{\prime}))-z(v^{\prime})\right\rVert.$ (115)

On the other hand, we have,

$\displaystyle\left\lVert z(\operatorname{pa}(v))-z(v)\right\rVert$
$\displaystyle\leq R_{\ell-1}$ $\displaystyle\left\lVert
z(\operatorname{pa}(v^{\prime}))-z(v^{\prime})\right\rVert$ $\displaystyle\leq
R_{\ell-1}$ (116)

since the assigned data of $\operatorname{pa}(v)$ is contained in the ball of
radius $R_{\ell+1}$ centered at $v$ by 13 of Algorithm 9, and $z(v)$ is in the
convex hull of the assigned data. Rearranging, we conclude that

$\|z(\operatorname{pa}(v))-z(\operatorname{pa}(v^{\prime}))\|\leq
2R_{\ell-1}+R_{\ell}=(2+1/a)R_{\ell-1}.$ (117)

From Lemma 30 and since for $a>1$, $2+1/a\leq 2a/(a-1)$, we conclude
$\operatorname{pa}(v^{\prime})\in\mathcal{R}(\operatorname{pa}(v))$. Without
loss of generality, assume that $v$ was added to the tree prior to
$v^{\prime}$, which is unambiguous because nodes are added sequentially. Then
by 13 of Algorithm 9, and since
$\operatorname{pa}(v^{\prime})\in\mathcal{R}(\operatorname{pa}(v))$, any node
selected in 8, and 10 if optional Lloyd’s averaging is performed, must satisfy

$\left\lVert z(v^{\prime})-z(v)\right\rVert>R_{\ell}$ (118)

which is a contradiction. ∎

######  Lemma 33.

We have

$|\operatorname{ch}(v)|\leq P^{\operatorname{ext},1/a}_{B_{1}}\leq(2a+1)^{d}.$
(119)

###### Proof.

Suppose $v\in V_{\ell}$. As $\mathcal{A}_{\ell}(v)$ is contained in the ball
of radius $R_{\ell}$ centered at $z(v)$, for any $w\in\operatorname{ch}(v)$,
$z(w)$ is contained in this ball. On the other hand, since $w\in V_{\ell+1}$
and since $\operatorname{sep}(z(V_{\ell+1}))\geq R_{\ell+1}$, balls of radius
$R_{\ell+1}$ centered at $\operatorname{ch}(v)$ form a $R_{\ell+1}$ (exterior)
packing of a sphere of radius $R_{\ell+1}$. Since packing numbers are
invariant to rescaling, this is the same as the packing number of the unit
sphere by spheres of radius $1/a$. The first inequality then follows from the
definition of an exterior packing, and the second from (vershynin18, 52,
Corollary 4.2.13). ∎

######  Lemma 34.

We have

$|\operatorname{ch}(v)|\leq
P^{\operatorname{ext},1/(4a)}_{B_{1}}\leq(8a+1)^{d}.$ (120)

###### Proof.

The proof is identical to Lemma 33, upon noting that $\mathcal{R}(v)$ is
contained in a ball of radius $4R_{\ell}$ centered at $z(v)$, as given in 16
of Algorithm 9. ∎

We are now ready to prove the main result.

*

###### Proof.

The claim follows by combining the above results. ∎

##  D Experimental Details

All experiments were run on a single Nvidia V100 GPU with 32GB RAM. For the
Wasserstein distance and condition number experiment of Figure 4, the dataset
of size $N=1000$ was generated by sampling from a GP prior with squared
exponential kernel, where inputs were uniformly sampled on the $d$ dimensional
cube $[-5,5]$ with $d$ equal to $1248$. The length scale of the kernel was set
to $0.5\sqrt{d}$. The parameters of the clustered-data approximation were set
via a cover tree with spatial resolutions $0.054.0$. Experiments for each $d$
were repeated 20 times to assess variability.

In the surface temperature geospatial data illustrative example shown of
Figure 6, we ran experiments with using the sparse Gaussian process regression
(SGPR) implementation from GPflow gpflow (35) and a GPflow-based
implementation of the clustered-data Gaussian process using stochastic maximum
marginal likelihood training described in Section B.2. For Hutchinson’s trace
estimator, we used 10 probe vectors. For the two-dimensional East Africa land
surface temperature dataset, we split the dataset into training and testing
sets of size $55884$ and $27525$, respectively. Both Gaussian process models
were initialized with $0.1$ likelihood noise and configured with a squared
exponential kernel with automatic relevance determination. The length scales
of the kernel were initialized to $1.0$. The inducing points were set using a
cover tree with the spatial resolution $0.030.09$. To facilitate comparisons,
in the sparse Gaussian process baseline, the inducing points were excluded
from the trainable parameter set. For both models hyperparameters were trained
using mini-batch stochastic optimization via the Adam optimizer, with constant
learning rate $0.01$ and batch size $1000$. We plot the posterior mean and
standard deviation from the clustered-data Gaussian process in Figure 8.

$\varepsilon=0.09$, $M=902$ |  |
---|---|---
$\varepsilon=0.06$, $M=1934$ |  |
$\varepsilon=0.03$, $M=6851$ |  |
|  Mean |  Standard Deviation
Figure 8: Predictions from the clustered-data Gaussian process approximation,
with varied spatial resolutions $\varepsilon$. For each spatial resolution, we
also display the resulting number of inducing points $M$.