# Pathwise Conditioning of Gaussian Processes

mes T. Wilson$*$<EMAIL_ADDRESS>
Imperial College London acheslav Borovitskiy$*$
<EMAIL_ADDRESS>
St. Petersburg State University and
St. Petersburg Department of Steklov Mathematical Institute of Russian Academy
of Sciences exander Terenin$*$<EMAIL_ADDRESS>
Imperial College London ter Mostowsky$*$<EMAIL_ADDRESS>
St. Petersburg State University rc Peter Deisenroth<EMAIL_ADDRESS>
Centre for Artificial Intelligence, University College London

###### Abstract

As Gaussian processes are used to answer increasingly complex questions,
analytic solutions become scarcer and scarcer. Monte Carlo methods act as a
convenient bridge for connecting intractable mathematical expressions with
actionable estimates via sampling. Conventional approaches for simulating
Gaussian process posteriors view samples as draws from marginal distributions
of process values at finite sets of input locations. This distribution-centric
characterization leads to generative strategies that scale cubically in the
size of the desired random vector. These methods are prohibitively expensive
in cases where we would, ideally, like to draw high-dimensional vectors or
even continuous sample paths. In this work, we investigate a different line of
reasoning: rather than focusing on distributions, we articulate Gaussian
conditionals at the level of random variables. We show how this _pathwise_
interpretation of conditioning gives rise to a general family of
approximations that lend themselves to efficiently sampling Gaussian process
posteriors. Starting from first principles, we derive these methods and
analyze the approximation errors they introduce. We, then, ground these
results by exploring the practical implications of pathwise conditioning in
various applied settings, such as global optimization and reinforcement
learning.

$*$Equal contribution.

Keywords: Gaussian processes, approximate posteriors, efficient sampling.

##  1 Introduction

In machine learning, the narrative of Gaussian processes (GPs) is dominated by
talk of distributions rasmussen06 (67). This view is often helpful and
convenient: a Gaussian process is a random function; however, seeing as we may
trivially marginalize out arbitrary subsets of this function, we can simply
focus on its behavior at a finite number of input locations. When dealing with
regression and classification problems, this reduction simplifies discourse
and expedites implementation by allowing us to work with joint distributions
at training and test locations instead of random functions.

Model-based learning and prediction generally service broader goals. For
example, when making decisions in the face of uncertainty, models enable us to
simulate the consequences of our actions. Decision-making, then, amounts to
optimizing the expectation of a simulated quantity of interest, such as a cost
or a reward. Be it for purposes of safety or for balancing trade-offs between
long-term and short-term goals, it is crucial that these simulations
faithfully portray both knowledge and uncertainty. Gaussian processes are
known to make accurate, well-calibrated predictions and, therefore, stand as
the model-of-choice in fields such as Bayesian optimization
shahriari2015taking (73), uncertainty quantification bect2012sequential (2),
and model-based reinforcement learning Deisenroth2015 (17).

Unfortunately, marginal distributions and simulations do not always go hand in
hand. When the quantity of interest is a function of a process value
$f(\boldsymbol{x}_{*})$ at an individual input location $\boldsymbol{x}_{*}$,
its expectation can sometimes be obtained analytically. Conversely, when this
quantity is a function of process values
$\boldsymbol{f}_{*}=f(\mathbf{X}_{*})$ at multiple locations $\mathbf{X}_{*}$,
its expectation is generally intractable. Rather than solving these integrals
directly in terms of marginal distributions $p(\boldsymbol{f}_{*})$, we
therefore estimate them by averaging over many simulations of
$\boldsymbol{f}_{*}$. Drawing $\boldsymbol{f}_{*}$ from
$p(\boldsymbol{f}_{*})$ takes $\mathcal{O}(*^{3})$ time, where
$*=|\mathbf{X}_{*}|$ is the number of input locations. Hence, distribution-
based approaches to sampling $\boldsymbol{f}_{*}$ quickly become untenable as
this number increases. In these cases, we may be better off thinking about GPs
from a perspective that naturally lends itself to sampling

In the early 1970s, one such view surfaced in the then nascent field of
geostatistics journel1978mining (37, 12). Instead of emphasizing the
statistical properties of Gaussian random variables, “conditioning by Kriging”
encourages us to think in terms of the variables themselves. We study the
broader implications of this paradigm shift to develop a general framework for
conditioning Gaussian processes at the level of random functions. Formulating
conditioning in terms of sample paths, rather than distributions, allows us to
separate out the effect of the prior from that of the data. By leveraging this
property, we can use _pathwise conditioning_ to efficiently approximate
function draws from GP posteriors. As we will see, working with sample paths
enables us to simulate process values $\boldsymbol{f}_{*}$ in $\mathcal{O}(*)$
time and brings with it a host of additional benefits.

The structure of the remaining text is as follows. Section 2 and Section 3
introduce pathwise conditioning of Gaussian random vectors and processes,
respectively. Section 4 surveys strategies for approximating function draws
from GP priors, while Section 5 discusses methods for mapping from prior to
posterior random variables. Section 6 studies the behavior of errors
introduced by different approximation techniques, and Section 7 complements
this theory with a bit of empiricism by exploring a number of examples.
Section 8 concludes.

#### Notation

By way of example, we denote matrices as $\mathbf{A}$ and vectors as
$\boldsymbol{a}$. We write $\boldsymbol{x}=\boldsymbol{a}\oplus\boldsymbol{b}$
for the direct sum (i.e. concatenation) of vectors $\boldsymbol{a}$ and
$\boldsymbol{b}$. Throughout, we use $|\cdot|$ to denote the cardinality of
sets and dimensionality of vectors. When dealing with covariance matrices
$\mathbf{\Sigma}=\operatorname{Cov}(\boldsymbol{x},\boldsymbol{x})$, we use
subscripts to identify corresponding blocks. For example,
$\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}=\operatorname{Cov}(\boldsymbol{a},\boldsymbol{b})$.
As shorthand, we denote the evaluation of a function
$f:\mathcal{X}\to\mathbb{R}$ at a finite set of locations
$\mathbf{X}_{*}\subset\mathcal{X}$ by the vector $\boldsymbol{f}_{*}$. Putting
these together, when dealing with random variables
$\boldsymbol{f}_{*}=f(\mathbf{X}_{*})$ and
$\boldsymbol{f}_{n}=f(\mathbf{X}_{n})$, we write
$\mathbf{K}_{*,n}=\operatorname{Cov}(\boldsymbol{f}_{*},\boldsymbol{f}_{n})$.

##  2 Conditioning Gaussian distributions and random vectors

A random vector $\boldsymbol{x}=(x_{1},\ldots,x_{n})\in\mathbb{R}^{n}$ is said
to be Gaussian if there exists a matrix $\mathbf{L}$ and vector
$\boldsymbol{\mu}$ such that

$\displaystyle\boldsymbol{x}$
$\displaystyle\overset{\mathop{}\\!\mathrm{d}}{=}\boldsymbol{\mu}+\mathbf{L}\boldsymbol{\zeta}$
$\displaystyle\boldsymbol{\zeta}$
$\displaystyle\sim\mathcal{N}(\boldsymbol{0},\mathbf{I}),$ (1)

where $\mathcal{N}(\boldsymbol{0},\mathbf{I})$ is the standard (multivariate)
normal distribution, the probability density function of which is given below.
Each such distribution is uniquely identified by its first two moments: its
mean $\boldsymbol{\mu}=\operatorname{\mathbb{E}}(\boldsymbol{x})$ and its
covariance
$\mathbf{\Sigma}=\operatorname{\mathbb{E}}[(\boldsymbol{x}-\boldsymbol{\mu})(\boldsymbol{x}-\boldsymbol{\mu})^{\top}]$.
Assuming it exists, the corresponding density function is defined as

$p(\boldsymbol{x})=\mathcal{N}(\boldsymbol{x}\mid\boldsymbol{\mu},\mathbf{\Sigma})=\frac{1}{\sqrt{\left|2\pi\mathbf{\Sigma}\right|}}\exp\left(-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^{\top}\mathbf{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\right).$
(2)

The representation of $\boldsymbol{x}$ given by (1) is commonly referred to as
its _location-scale_ form and stands as the most widely used method for
generating Gaussian random vectors. Since $\boldsymbol{\zeta}$ has identity
covariance, any matrix square root of $\mathbf{\Sigma}$, such as its Cholesky
factor $\mathbf{L}$ with $\mathbf{\Sigma}=\mathbf{L}\mathbf{L}^{\top}$, may be
used to draw $\boldsymbol{x}$ as prescribed by (1).

Here, we focus on multivariate cases $n>1$ and investigate different ways of
reasoning about random variables
$\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta}$ for non-trivial
partitions $\boldsymbol{x}=\boldsymbol{a}\oplus\boldsymbol{b}$.

###  2.1 Distributional conditioning

The quintessential approach to deriving the distribution of $\boldsymbol{a}$
subject to the condition $\boldsymbol{b}=\boldsymbol{\beta}$ begins by
employing the usual set of matrix identities to factor $p(\boldsymbol{b})$
from $p(\boldsymbol{a},\boldsymbol{b})$. Applying Bayes’ rule,
$p(\boldsymbol{b})$ then cancels out and
$p(\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta})$ is identified as the
remaining term—namely, the Gaussian distribution
$\mathcal{N}(\boldsymbol{\mu}_{\boldsymbol{a}\mid\boldsymbol{\beta}},\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{a}\mid\boldsymbol{\beta}})$
with moments

$\displaystyle\boldsymbol{\mu}_{\boldsymbol{a}\mid\boldsymbol{\beta}}$
$\displaystyle=\boldsymbol{\mu}_{\boldsymbol{a}}+\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{-1}(\boldsymbol{\beta}-\boldsymbol{\mu}_{\boldsymbol{b}})$
$\displaystyle\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{a}\mid\boldsymbol{\beta}}$
$\displaystyle=\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{a}}-\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}^{\vphantom{{-1}}}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{-1}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{a}}^{\vphantom{{-1}}}.$
(3)

Having obtained this conditional distribution, we can now generate
$\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta}$ by computing a matrix
square root of
$\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{a}\mid\boldsymbol{\beta}}$ and
constructing a location-scale transform (1).

Due to their emphasis of conditional distributions, we refer to methods that
represent or generate a random variable
$\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta}$ by way of
$p(\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta})$ as being
_distributional_ in kind. This approach to conditioning is not only standard,
but particularly natural when quantities of interest may be derived
analytically from $p(\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta})$.
Many quantities, such as expectations of nonlinear functions, cannot be
deduced analytically from
$p(\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta})$ alone, however. In
these case, we must instead work with realizations of
$\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta}$. Since the cost of
obtaining a matrix square root of
$\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{a}\mid\boldsymbol{\beta}}$ scales
cubically in $|\boldsymbol{a}|$, distributional approaches to evaluating these
quantities struggle to accommodate high-dimensional random vectors. To address
this issue, we now consider Gaussian conditioning in another light.

###  2.2 Pathwise conditioning

Figure 1: Visualization of Matheron’s update rule for a bivariate normal
distribution with correlation coefficient $\rho=0.75$. _Left:_ Draws from
$p(\boldsymbol{a},\boldsymbol{b})$ are shown alongside the marginal
distributions of $\boldsymbol{a}$ and $\boldsymbol{b}$. _Right:_ Theorem 1 is
used to update samples shown on the left subject to the condition
$\boldsymbol{b}=\boldsymbol{\beta}$. This process is illustrated in full for a
particular draw. _Top right:_ the empirical distribution of the update samples
is compared with $p(\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta})$.

Instead of taking a _distribution-first_ stance on Gaussian conditionals, we
may think of conditioning directly in terms of random variables. In this
_variable-first_ paradigm, we will explicitly map samples from the prior to
draws from a posterior and let the corresponding relationship between
distributions follow implicitly. Throughout this work, we investigate this
notion of _pathwise conditioning_ through the lens of the following result.

######  Theorem 1 (Matheron’s Update Rule)

Let $\boldsymbol{a}$ and $\boldsymbol{b}$ be jointly Gaussian, centered random
variables. Then, the random variable $\boldsymbol{a}$ conditional on
$\boldsymbol{b}=\boldsymbol{\beta}$ may be expressed as

$(\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta})\overset{\mathop{}\\!\mathrm{d}}{=}\boldsymbol{a}+\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}^{\vphantom{-1}}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{-1}(\boldsymbol{\beta}-\boldsymbol{b}).$
(4)

Proof  Comparing the mean and covariance on both sides immediately affirms the
result

$\displaystyle\operatorname{\mathbb{E}}\big{(}\boldsymbol{a}+\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}^{\vphantom{-1}}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{-1}(\boldsymbol{\beta}-\boldsymbol{b})\big{)}$
$\displaystyle\operatorname{Cov}\big{(}\boldsymbol{a}+\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}^{\vphantom{-1}}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{-1}(\boldsymbol{\beta}-\boldsymbol{b})\big{)}$
$\displaystyle\quad=\boldsymbol{\mu}_{\boldsymbol{a}}^{\vphantom{-1}}+\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}^{\vphantom{-1}}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{-1}(\boldsymbol{\beta}-\boldsymbol{\mu}_{\boldsymbol{b}}^{\vphantom{-1}})$
$\displaystyle\quad=\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{a}}^{\vphantom{-1}}\\!+\\!\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}^{\vphantom{-1}}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{-1}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{\vphantom{-1}}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{-1}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{a}}^{\vphantom{-1}}\\!-\\!2\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}^{\vphantom{-1}}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{-1}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{a}}^{\vphantom{-1}}$
(5)
$\displaystyle\quad=\operatorname{\mathbb{E}}(\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta})$
$\displaystyle\quad=\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{a}}^{\vphantom{-1}}\\!-\\!\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}^{\vphantom{-1}}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{-1}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{a}}^{\vphantom{-1}}=\operatorname{Cov}(\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta}).$

This observation leads to a straightforward, alternative recipe for generating
$\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta}$: first, draw
$\boldsymbol{a},\boldsymbol{b}\sim p(\boldsymbol{a},\boldsymbol{b})$; then,
update this sample according to (4). Compared to the location-scale approach
discussed in Section 2.1, a key difference is that we now sample _before_
conditioning, rather than after. Figure 1 visualizes the deterministic process
of updating previously generated draws from the prior subject to the condition
$\boldsymbol{b}=\boldsymbol{\beta}$.

At first glance, Matheron’s update rule may seem more like an interesting
footnote than a valuable tool. Indeed, the conventional strategy for sampling
$\boldsymbol{a},\boldsymbol{b}$ (which requires us to take a matrix square
root of $\mathbf{\Sigma}$) is more expensive than that for generating
$\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta}$. We will discuss this
matter in detail in the later sections. For now, however, let us strengthen
our intuition by delving deeper into this theorem’s function-analytic origins.

###  2.3 Deriving pathwise conditioning via conditional expectations

Here, we overview the precise formalism that gives rise to the pathwise
approach to conditioning Gaussian random variables and show how to _derive_
this result from first principles. Throughout this section, we take
$\boldsymbol{a}\in\mathbb{R}^{m}$ and $\boldsymbol{b}\in\mathbb{R}^{n}$ to be
centered random vectors defined on the same probability space.

The core idea is to decompose $\boldsymbol{a}$ as the sum of two independent
terms—one that depends on $\boldsymbol{b}$ and one that does not—and represent
$\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta}$ by conditioning both
terms on $\boldsymbol{b}=\boldsymbol{\beta}$. We first prove that conditioning
this additive decomposition of $\boldsymbol{a}$ is simple and intuitive.

###### Lemma 2

Consider three random vectors $\boldsymbol{a}\in\mathbb{R}^{m}$,
$\boldsymbol{b}\in\mathbb{R}^{n}$, $\boldsymbol{c}\in\mathbb{R}^{m}$ such that

$\boldsymbol{a}\overset{\mathop{}\\!\mathrm{d}}{=}f(\boldsymbol{b})+\boldsymbol{c},$
(6)

where $f$ is a measurable function of $\boldsymbol{b}$ and where
$\boldsymbol{b}$ is independent of $\boldsymbol{c}$. Then,

$\big{(}\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta}\big{)}\overset{\mathop{}\\!\mathrm{d}}{=}f(\boldsymbol{\beta})+\boldsymbol{c}.$
(7)

Proof  Let $\pi_{\boldsymbol{x}}$ denote the distribution of a generic random
variable $\boldsymbol{x}$. Further, let
$\pi_{\boldsymbol{a}\mid\boldsymbol{b}}(\cdot\mid\cdot)$ be the (regular)
conditional probability measure given by disintegration111See discussion and
details on disintegration by chang1997conditioning (9, 38). of
$(\boldsymbol{a},\boldsymbol{b})$, such that

$\int_{B}\pi_{\boldsymbol{a}\mid\boldsymbol{b}}(A\mid\boldsymbol{\beta})\mathop{}\\!\mathrm{d}\pi_{\boldsymbol{b}}(\boldsymbol{\beta})=\operatorname{\mathbb{P}}\left(\boldsymbol{a}\in
A,\boldsymbol{b}\in B\right)$ (8)

for measurable sets $A\subseteq\mathbb{R}^{m}$, $B\subseteq\mathbb{R}^{n}$.
When $\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta}$ is represented per
(7), we have

$\displaystyle\begin{split}\int_{B}\operatorname{\mathbb{P}}\left(f(\boldsymbol{\beta})+\boldsymbol{c}\in
A\right)\mathop{}\\!\mathrm{d}\pi_{\boldsymbol{b}}(\boldsymbol{\beta})&=\int_{B}\left(\int_{\mathbb{R}^{m}}\mathbbold{1}_{\left\\{f(\boldsymbol{\beta})+\boldsymbol{\varsigma}\in
A)\right\\}}\mathop{}\\!\mathrm{d}\pi_{\boldsymbol{c}}(\boldsymbol{\varsigma})\right)\mathop{}\\!\mathrm{d}\pi_{\boldsymbol{b}}(\boldsymbol{\beta})\\\
&=\int_{\mathbb{R}^{m}\times\mathbb{R}^{n}}\mathbbold{1}_{\left\\{f(\boldsymbol{\beta})+\boldsymbol{\varsigma}\in
A,\boldsymbol{\beta}\in
B)\right\\}}\mathop{}\\!\mathrm{d}\pi_{\boldsymbol{b},\boldsymbol{c}}(\boldsymbol{\beta},\boldsymbol{\varsigma})\\\
&=\operatorname{\mathbb{P}}\left(f(\boldsymbol{b})+\boldsymbol{c}\in
A,\boldsymbol{b}\in B\right)=\operatorname{\mathbb{P}}\left(\boldsymbol{a}\in
A,\boldsymbol{b}\in B\right),\end{split}$ (9)

where we have begun by expressing probabilities as integrals of indicator
functions, before using Tonelli’s theorem and independence to express the
iterated integral as the double integral over the joint probability measure
$\pi_{\boldsymbol{b},\boldsymbol{c}}(\boldsymbol{\beta},\boldsymbol{\varsigma})$.
Comparing the left-hand sides of (8) and (9) affirms the claim.

In words, Lemma 2 tells us that for suitably chosen functions $f$, the act of
conditioning $\boldsymbol{a}$ on $\boldsymbol{b}=\boldsymbol{\beta}$ amounts
to adding an random variable $\boldsymbol{c}$ to a deterministic
transformation $f(\boldsymbol{\beta})$ of the outcome $\boldsymbol{\beta}$.
For this statement to hold, we require the _residual_
$\boldsymbol{c}=\boldsymbol{a}-f(\boldsymbol{b})$ induced by $f$ to be
independent of $\boldsymbol{b}$. Fortunately, such a function $f$ is well-
known in the special case of jointly Gaussian random variables—namely, the
_conditional expectation_
$f:\boldsymbol{b}\mapsto\operatorname{\mathbb{E}}(\boldsymbol{a}\mid\boldsymbol{b})$.

For square-integrable random variables, the conditional expectation of
$\boldsymbol{a}$ given $\boldsymbol{b}$ is defined as the (almost surely)
unique solution to the minimization problem

$\displaystyle\operatorname{\mathbb{E}}(\boldsymbol{a}\mid\boldsymbol{b})$
$\displaystyle=\operatorname*{\arg\min}_{f\in\mathcal{F}}\operatorname{\mathbb{E}}\left\lVert\boldsymbol{a}-f(\boldsymbol{b})\right\rVert^{2},$
(10)

where $\mathcal{F}$ denotes the set of all Borel-measurable functions
$f:\mathbb{R}^{n}\to\mathbb{R}^{m}$ (kallenberg2006, 38, Chapter 6). Put
simply, $\operatorname{\mathbb{E}}(\boldsymbol{a}\mid\boldsymbol{b})$ is the
measurable function of $\boldsymbol{b}$ that best predicts $\boldsymbol{a}$ in
the sense of minimizing the mean-square error (10). This characterization of
the conditional expectation is equivalent to defining it as the orthogonal
projection of $\boldsymbol{a}$ onto the $\sigma$-algebra generated by
$\boldsymbol{b}$, denoted $\sigma(\boldsymbol{b})$. Consequently, a necessary
and sufficient condition for
$\operatorname{\mathbb{E}}(\boldsymbol{a}\mid\boldsymbol{b})\in\mathcal{F}$ to
uniquely solve (10) is that the residual
$\boldsymbol{c}=\boldsymbol{a}-\operatorname{\mathbb{E}}(\boldsymbol{a}\mid\boldsymbol{b})$
be orthogonal to all $\sigma(\boldsymbol{b})$-measurable random variables
(luenberger1997optimization, 51, 50). Here, _orthogonality_ can be understood
as the absence of correlation, which (for jointly Gaussian random variables)
implies independence. As a result, we may satisfy the assumptions of Lemma 2
by writing

$\boldsymbol{a}=\operatorname{\mathbb{E}}(\boldsymbol{a}\mid\boldsymbol{b})+\boldsymbol{c},$
(11)

such that $\boldsymbol{a}$ decomposes into a function of $\boldsymbol{b}$ and
an independent variable
$\boldsymbol{c}=\boldsymbol{a}-\operatorname{\mathbb{E}}(\boldsymbol{a}\mid\boldsymbol{b})$.

As a final remark, we may also use these principles to concisely derive the
conditional expectation for jointly Gaussian random variables. For now,
suppose that the conditional expectation is a linear function of
$\boldsymbol{b}$, i.e. that
$\operatorname{\mathbb{E}}(\boldsymbol{a}\mid\boldsymbol{b})=\mathbf{S}\boldsymbol{b}$
for some matrix $\mathbf{S}\in\mathbb{R}^{m\times n}$. To satisfy the
orthogonality condition of (10), we require
$\operatorname{Cov}(\boldsymbol{a}-\mathbf{S}\boldsymbol{b},\boldsymbol{b})=\mathbf{0}$,
implying that
$\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}-\mathbf{S}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}=\mathbf{0}$.
Rearranging terms and solving for $\mathbf{S}$ gives
$\mathbf{S}=\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}^{\vphantom{-1}}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{-1}\boldsymbol{b}$.
With this expression in hand, to show that linearity was assumed without loss
of generality, write
$\boldsymbol{a}=\mathbf{S}\boldsymbol{b}+\boldsymbol{a}-\mathbf{S}\boldsymbol{b}$,
which we may express as as
$\boldsymbol{a}=\mathbf{S}\boldsymbol{b}+\boldsymbol{c}$. Taking the
conditional expectation of both sides, we may directly calculate
$\operatorname{\mathbb{E}}(\boldsymbol{a}\mid\boldsymbol{b})$ by writing

$\operatorname{\mathbb{E}}(\boldsymbol{a}\mid\boldsymbol{b})=\operatorname{\mathbb{E}}(\mathbf{S}\boldsymbol{b}+\boldsymbol{c}\mid\boldsymbol{b})=\underbracket{\operatorname{\mathbb{E}}(\mathbf{S}\boldsymbol{b}\mid\boldsymbol{b})}_{\smash{\mathbf{S}\boldsymbol{b}}}+\underbracket{\operatorname{\mathbb{E}}(\boldsymbol{c})}_{\smash{\boldsymbol{0}}}=\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}^{\vphantom{-1}}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{-1}\boldsymbol{b},$
(12)

where we have used linearity of conditional expectation, followed by
independence of $\boldsymbol{c}$ and $\boldsymbol{b}$ to go from the second to
the third expression. We now revisit Theorem 1.

######  Theorem 1 (Matheron’s Update Rule)

Let $\boldsymbol{a}$ and $\boldsymbol{b}$ be jointly Gaussian, centered random
vectors. Then, the random vector $\boldsymbol{a}$ conditional on
$\boldsymbol{b}=\boldsymbol{\beta}$ may be expressed as

$(\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta})\overset{\mathop{}\\!\mathrm{d}}{=}\boldsymbol{a}+\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}^{\vphantom{-1}}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{-1}(\boldsymbol{\beta}-\boldsymbol{b}).$
(4)

Proof  With
$\smash{\boldsymbol{c}=\boldsymbol{a}-\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}^{\vphantom{-1}}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{-1}\boldsymbol{b}}$,
begin by writing

$\boldsymbol{a}=\operatorname{\mathbb{E}}(\boldsymbol{a}\mid\boldsymbol{b})+\left(\boldsymbol{a}-\operatorname{\mathbb{E}}(\boldsymbol{a}\mid\boldsymbol{b})\right)=\smash{\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}^{\vphantom{-1}}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{-1}\boldsymbol{\boldsymbol{b}}+\boldsymbol{c}.}$
(13)

Since $\boldsymbol{b}$ and $\boldsymbol{c}$ are jointly Gaussian but
uncorrelated, it follows that they are independent. Setting
$f(\boldsymbol{b})=\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}^{\vphantom{-1}}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{-1}\boldsymbol{\boldsymbol{b}}$
and using Lemma 2 to condition both sides on
$\boldsymbol{b}=\boldsymbol{\beta}$ gives

$(\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta})\overset{\mathop{}\\!\mathrm{d}}{=}\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}^{\vphantom{-1}}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{-1}\boldsymbol{\boldsymbol{\beta}}+\big{(}\boldsymbol{a}-\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}^{\vphantom{-1}}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{-1}\boldsymbol{b}\big{)}=\boldsymbol{a}+\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}^{\vphantom{-1}}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{-1}(\boldsymbol{\beta}-\boldsymbol{b}).$
(14)

Hence, the claim follows.

In summary, we have shown that Matheron’s update rule (Theorem 1) is a direct
consequence of the fact that a Gaussian random variable $\boldsymbol{a}$
conditioned on the outcome $\boldsymbol{\beta}$ of another (jointly) Gaussian
random variable $\boldsymbol{b}$ may be expressed as the sum of two
independent terms: the conditional expectation
$\operatorname{\mathbb{E}}(\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta})$
evaluated at $\boldsymbol{\beta}$ and the residual
$\boldsymbol{c}=\boldsymbol{a}-\operatorname{\mathbb{E}}(\boldsymbol{a}\mid\boldsymbol{b})$.
Rearranging these terms gives (4).

With these ideas in mind, we are now ready to explore this work’s primary
theme: Matheron’s update rule enables us to decompose
$\boldsymbol{a}\mid\boldsymbol{b}=\boldsymbol{\beta}$ into the prior random
variable $\boldsymbol{a}$ and a data-driven update
$\mathbf{\Sigma}_{\boldsymbol{a},\boldsymbol{b}}^{\vphantom{-1}}\mathbf{\Sigma}_{\boldsymbol{b},\boldsymbol{b}}^{-1}(\boldsymbol{\beta}-\boldsymbol{b})$
that explicitly corrects for the error in the coinciding value of
$\boldsymbol{b}$ given the condition $\boldsymbol{b}=\boldsymbol{\beta}$.
Hence, Theorem 1 provides an explicit means of separating out the influence of
the prior from that of the data. We now proceed to investigate the
implications of pathwise conditioning for Gaussian processes.

##  3 Conditioning Gaussian processes and random functions

A Gaussian process (GP) is a random function $f:\mathcal{X}\to\mathbb{R}$,
such that, for any finite collection of points $\mathbf{X}\subset\mathcal{X}$,
the random vector $\boldsymbol{f}=f(\mathbf{X})$ follows a Gaussian
distribution. Such a process is uniquely identified by a mean function
$\mu:\mathcal{X}\to\mathbb{R}$ and a positive semi-definite kernel
$k:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$. Hence, if
$f\sim\mathcal{GP}(\mu,k)$, then
$\boldsymbol{f}\sim\mathcal{N}(\boldsymbol{\mu},\mathbf{K})$ is multivariate
normal with mean $\boldsymbol{\mu}=\mu(\mathbf{X})$ and covariance
$\mathbf{K}=k(\mathbf{X},\mathbf{X})$.

Throughout this section, we investigate different ways of reasoning about the
random variable $\boldsymbol{f}_{*}\mid\boldsymbol{f}_{n}=\boldsymbol{y}$ for
some non-trivial partition
$\boldsymbol{f}=\boldsymbol{f}_{n}\oplus\boldsymbol{f}_{*}$. Here,
$\boldsymbol{f}_{n}=f(\mathbf{X}_{n})$ are process values at a set of training
locations $\mathbf{X}_{n}\subset\mathbf{X}$ where we would like to introduce a
condition $\boldsymbol{f}_{n}=\boldsymbol{y}$, while
$\boldsymbol{f}_{*}=f(\mathbf{X}_{*})$ are process values at a set of test
locations $\mathbf{X}_{*}\subset\mathbf{X}$ where we would like to obtain a
random variable $\boldsymbol{f}_{*}\mid\boldsymbol{f}_{n}=\boldsymbol{y}$.
Mirroring Section 2, we begin by reviewing distributional conditioning, before
examining its pathwise counterpart.

###  3.1 Distributional conditioning

As in finite-dimensional cases, we may obtain
$\boldsymbol{f}_{*}\mid\boldsymbol{y}$ by first finding its conditional
distribution. Since process values $(\boldsymbol{f}_{n},\boldsymbol{f}_{*})$
are defined as jointly Gaussian, this procedure closely resembles that of
Section 2.1: we factor out the marginal distribution of $\boldsymbol{f}_{n}$
from the joint distribution $p(\boldsymbol{f}_{n},\boldsymbol{f}_{*})$ and,
upon canceling, identify the remaining distribution as
$p(\boldsymbol{f}_{*}\mid\boldsymbol{y})$. Having done so, we find that the
conditional distribution is the Gaussian
$\mathcal{N}(\boldsymbol{\mu}_{*\mid\boldsymbol{y}},\mathbf{K}_{*,*\mid\boldsymbol{y}})$
with moments

$\displaystyle\boldsymbol{\mu}_{*\mid\boldsymbol{y}}$
$\displaystyle=\boldsymbol{\mu}_{*}+\mathbf{K}_{*,n}^{\vphantom{-1}}\mathbf{K}_{n,n}^{-1}(\boldsymbol{y}-\boldsymbol{\mu}_{n})$
$\displaystyle\mathbf{K}_{*,*\mid\boldsymbol{y}}^{\vphantom{-1}}$
$\displaystyle=\mathbf{K}_{*,*}-\mathbf{K}_{*,n}^{\vphantom{-1}}\mathbf{K}_{n,n}^{-1}\mathbf{K}_{n,*}^{\vphantom{-1}}.$
(15)

As before, we may now generate $\boldsymbol{f}_{*}\mid\boldsymbol{y}$ in
$\mathcal{O}(*^{3})$ time using a location-scale transform (1).

This strategy for sampling Gaussian process posteriors is subtly different
from the one given in Section 2.1. A Gaussian process is a random function,
and conditioning on $\boldsymbol{f}_{n}=\boldsymbol{y}$ does not change this
fact. Unfortunately, (conditional) distributions over infinite-dimensional
objects can be difficult to manipulate in practice. Distributional approaches,
therefore, focus on finite-dimensional subsets
$\boldsymbol{f}=\boldsymbol{f}_{n}\oplus\boldsymbol{f}_{*}$, while
marginalizing out the remaining process values. Doing so allows them to
perfectly describe the random variable $\boldsymbol{f}_{*}\mid\boldsymbol{y}$
via its mean and covariance (15).

When it comes to sampling $\boldsymbol{f}_{*}\mid\boldsymbol{y}$, however,
these approaches have clear limitations. As discussed previously, a key issue
is that their $\mathcal{O}(*^{3})$ time complexity restricts them to problems
that only require us to jointly simulate process values at a manageable number
of test locations (up to several thousand). In some senses, this condition is
fairly generous. After all, we are often only asked to generate a handful of
process values at a time. Still, other problems effectively require us to
realize $f\mid\boldsymbol{y}$ in its entirety. Similar issues arise when
$\mathbf{X}_{*}$ is not defined in advance, such as when gradient information
is used to adaptively determine the locations at which to jointly sample the
posterior. In these cases and more, we would ideally like to sample actual
functions that we can efficiently evaluate and automatically differentiate at
arbitrary test locations. To this end, we now examine the direct approach to
conditioning draws of $f\sim\mathcal{GP}(\mu,k)$.

###  3.2 Pathwise Conditioning

Figure 2: Visual guide for pathwise conditioning of Gaussian processes.
_Left:_ The residual $\boldsymbol{y}-\boldsymbol{f}_{n}$ (dashed black) of a
draw $f\sim\mathcal{GP}(0,k)$, shown in orange, given observations
$\boldsymbol{y}$ (black). _Middle:_ A pathwise update (purple) is constructed
in accordance with Corollary 4. _Right:_ Prior and update are combined to
represent conditional (blue). Empirical moments (light blue) of $10^{5}$
conditioned paths are compared with those of the model (dashed black). The
sample average, which matches the posterior mean, has been omitted for
clarity.

Examining the pathwise update given by Theorem 1, it is natural to suspect
that an analogous statement holds for Gaussian processes. A quick check
confirms this hypothesis.

###### Corollary 4

For a Gaussian process $f\sim\mathcal{GP}(\mu,k)$ with marginal
$\boldsymbol{f}_{n}=f(\mathbf{X}_{n})$, the process conditioned on
$\boldsymbol{f}_{n}=\boldsymbol{y}$ may be expressed as

$\underbracket{(f\mid\boldsymbol{y})(\cdot)\vphantom{\mathbf{K}_{n,n}^{-1}}}_{\operatorname{conditional}}\overset{\mathop{}\\!\mathrm{d}}{=}\underbracket{f(\cdot)\vphantom{\mathbf{K}_{n,n}^{-1}}}_{\operatorname{prior}}+\underbracket{k(\cdot,\mathbf{X}_{n})\mathbf{K}_{n,n}^{-1}(\boldsymbol{y}-\boldsymbol{f}_{n})}_{\operatorname{update}}.$
(16)

Proof  Follows by applying Theorem 1 to an arbitrary set of locations.
Figure 2 acts a visual guide to Corollary 4. From left to right, we begin by
generating a realization of $f\sim\mathcal{GP}(\mu,k)$ using methods that will
soon be introduced in Section 4. Having obtained a sample path, we then use
the pathwise update (16) to define a function
$k(\cdot,\mathbf{X}_{n})\mathbf{K}^{-1}_{n,n}(\boldsymbol{y}-\boldsymbol{f}_{n})$
to account for the residual $\boldsymbol{y}-\boldsymbol{f}_{n}$. Adding these
two functions together produces a draw from a GP posterior, the behavior of
which is shown on the right. Whereas distributionally conditioning on
$\boldsymbol{f}_{n}=\boldsymbol{y}$ in (15) tells us how the GP’s statistic
properties change, pathwise conditioning (16) tells us what happens to
individual sample paths. This paradigm shift echoes the running theme:
Gaussian (process) conditionals can be directly viewed in terms of random
variables. The power of Corollary 4 is that it impacts _how_ we think about
Gaussian process posteriors and, therefore, _what_ we do with them.

Having said this, there are several hurdles that we must overcome in order to
use the pathwise update (16) in the real world. First, we are typically unable
to practically sample functions $f\sim\mathcal{GP}(\mu,k)$ from (non-
degenerate) Gaussian process priors exactly. A Gaussian process can generally
be written as a linear combination of elementary _basis functions_. When the
requisite number of basis functions is infinite, however, evaluating this
linear combination is usually impossible. In Section 4, we will therefore
investigate different ways of approximating $f(\cdot)$ using a finite number
of operations.

Second, we incur $\mathcal{O}(n^{3})$ time complexity when naïvely carrying
out (16), due to the need to solve the linear system of equations
$\mathbf{K}_{n,n}\boldsymbol{v}=\boldsymbol{y}-\boldsymbol{f}_{n}$ for a
vector $\boldsymbol{v}\in\mathbb{R}^{n}$ such that

$(f\mid\boldsymbol{y})(\cdot)\overset{\mathop{}\\!\mathrm{d}}{=}f(\cdot)+\underbracket{\sum_{i=1}^{n}v_{i}k(\cdot,\boldsymbol{x}_{i})}_{\mathclap{n\mathrel{\text{missing}}{-dimensionalbasis}}}.$
(17)

Here, we have re-expressed the matrix-vector product in (16) as an expansion
with respect to the canonical basis functions $k(\cdot,\boldsymbol{x}_{i})$
centered at training locations $\boldsymbol{x}_{i}\in\mathbf{X}_{n}$. For
large training sets $\left(\boldsymbol{x}_{i},y_{i}\right)_{i=1}^{n}$, direct
application of (16) may prove prohibitively expensive. By the same token, the
stated pathwise update does not hold when outcomes $\boldsymbol{y}$ are not
defined as realizations of process values $\boldsymbol{f}_{n}$. In Section 5,
we will consider various means of resolving these challenges and ones like
them.

###  3.3 Historical remarks

Prior to continuing, we pause to reflect on the historical developments that
have paved the way for this work. In a 2005 tribute to geostatistics pioneer
Georges Matheron, chiles2005prediction (13) comment that

> [Matheron’s update rule] is nowhere to be found in Matheron’s entire
> published works, as he merely regarded it as an immediate consequence of the
> orthogonality of the [conditional expectation] and the [residual process].

As if to echo this very sentiment, doucet10 (20) begins a much appreciated
technical note on the subject of Theorem 1 with the remark

> This note contains no original material and will never be submitted anywhere
> for publication. However it might be of interest to people working with
> [Gaussian processes] so I am making it publicly available.

The presiding opinion, therefore, seems to be that Matheron’s update rule is
_too simple_ to warrant extended study. Indeed, Theorem 1 is exceedingly
straightforward to verify. As is often the case, however, this result is
harder to discover if one is not already aware of its existence. This dilemma
may help to explain why Matheron’s update rule is absent from standard machine
learning texts. By deriving this result from first principles in Section 2.3,
we hope to encourage fellow researchers to explore the strengths (and
weaknesses) of the pathwise viewpoint espoused here.

We are not the first to have realized the practical implications of pathwise
conditioning for GPs. Corollary 4 is relatively well-known in geostatistics
journel1978mining (37, 16, 23, 12). Similarly, oliver1996conditional (58)
discusses Matheron’s update rule for Gaussian likelihoods (Section 5.1). Along
the same lines, closely related ideas were rediscovered in the 1990s with
applications to astrophysics. In particular, hoffman91 (35) propose the use of
spectral approximations to stationary priors (Section 4.2) in conjunction with
canonical pathwise updates (17).

Nevertheless, these formulae are seldom seen in machine learning. We hope to
systematically organize these findings (along with our own) and communicate
them to a general audience of theorists and practitioners alike. The following
sections therefore catalog various notable approaches to representing Gaussian
process priors and pathwise updates.

##  4 Sampling functions from Gaussian process priors

The pathwise representation of GP posteriors described in the Section 3.2
allows us to represent $f\mid\boldsymbol{y}$ by transforming a draw of
$f\sim\mathcal{GP}(0,k)$. When interpreted as a generative strategy, this
approach to sampling can only be deemed _efficient_ if the tasks of realizing
the prior and performing the update both scale favorably in the total number
of locations $|\mathbf{X}|=|\mathbf{X}_{n}|+|\mathbf{X}_{*}|$. Half of the
battle is, therefore, to obtain faithful but affordable draws of $f$.
Fortunately, GP priors often exhibit convenient mathematical properties not
present in their posteriors, which can be utilized to sample them efficiently.

We focus on methods for generating random _functions_ that we may evaluate at
arbitrary locations $\boldsymbol{x}\in\mathcal{X}$ in $\mathcal{O}(1)$ time
and whose marginal distributions approximate those of
$f\sim\mathcal{GP}(0,k)$. Conceptually, techniques discussed throughout this
section will approximate GP priors as random linear combinations of suitably
chosen basis functions $\boldsymbol{\phi}=(\phi_{1},\ldots,\phi_{\ell})$.
Specifically, we will focus on _Bayesian linear models_ with Gaussian random
weights

$\displaystyle\tilde{f}(\cdot)$
$\displaystyle=\sum_{i=1}^{\ell}w_{i}\phi_{i}(\cdot)$
$\displaystyle\boldsymbol{w}$
$\displaystyle\sim\mathcal{N}(\boldsymbol{0},\mathbf{\Sigma}_{\boldsymbol{w}}),$
(18)

where the covariance of weights $\boldsymbol{w}$ will vary by case. Notice
that, for any finite collection of points $\mathbf{X}\subset\mathcal{X}$, the
random vector $\boldsymbol{\tilde{f}}=\tilde{f}(\mathbf{X})$ follows the
Gaussian distribution
$\mathcal{N}(\boldsymbol{0},\mathbf{\Phi}^{\vphantom{\top}}\mathbf{\Sigma}_{\boldsymbol{w}}\mathbf{\Phi}^{\top})$,
where $\mathbf{\Phi}=\boldsymbol{\phi}(\mathbf{X})$ is a
$|\mathbf{X}|\times\ell$ matrix of features. By design then, $\tilde{f}$ is a
Gaussian process. rasmussen06 (67) refer to (18) as the _weight-space_ view of
GPs.

From this perspective, the task of efficiently sampling the prior $\tilde{f}$
reduces to one of generating random weights $\boldsymbol{w}$. In practice,
$\mathbf{\Sigma}_{\boldsymbol{w}}$ is typically diagonal, thereby enabling us
to sample $\tilde{f}$ in $\mathcal{O}(\ell)$ time. We stress that, for any
draw of $\boldsymbol{w}$, the corresponding realization of $\tilde{f}$ is
simply a deterministic function. In particular, we incur $\mathcal{O}(1)$ cost
for evaluating $\tilde{f}(\boldsymbol{x})$ and may readily differentiate this
term with respect to $\boldsymbol{x}$ (or other parameters of interest).

Below, we review popular strategies for obtaining Bayesian linear models such
that $\tilde{f}\overset{\mathop{}\\!\mathrm{d}}{\approx}f$. Our presentation
is intended to communicate different angles for attacking this problem and is
by no means exhaustive. To set the scene for these approaches, we begin by
recounting some properties of the gold standard: location-scale methods.

###  4.1 Location-scale transformations

Location-scale methods $\eqref{eqn:location_scale}$ are the most widely used
approach for generating Gaussian random vectors. These generative strategies
are _exact_ (up to machine precision). Given locations $\mathbf{X}$, we may
simulate $\boldsymbol{f}=f(\mathbf{X})$ in location-scale fashion

$\displaystyle f(\mathbf{X})$
$\displaystyle\overset{\mathop{}\\!\mathrm{d}}{=}\mathbf{K}^{\nicefrac{{1}}{{2}}}\boldsymbol{\zeta}$
$\displaystyle\boldsymbol{\zeta}$
$\displaystyle\sim\mathcal{N}(\boldsymbol{0},\mathbf{I})$ (19)

by multiplying a square root covariance matrix
$\mathbf{K}^{\nicefrac{{1}}{{2}}}$ by a standard normal vector
$\boldsymbol{\zeta}$.

While (19) rightfully stands as the method of choice for many problems, it is
not without shortcoming. Chief among these issues is the fact that algorithms
for obtaining a matrix square root of $\mathbf{K}$ scale cubically in
$|\mathbf{X}|$. In most cases, this limits the use of location-scale
approaches to cases where the length of the desired Gaussian random vector is
manageable (up to several thousand). This overhead can be interpreted to mean
that we incur $\mathcal{O}(i^{2})$ cost for realizing the $i$-th element of
$\boldsymbol{f}$, which leads us to our second issue: reusing a draw of
$\boldsymbol{f}_{n}$ to efficiently generate the remainder of
$\boldsymbol{f}=\boldsymbol{f}_{n}\oplus\boldsymbol{f}_{*}$ requires us to
sample from the conditional distribution

$\boldsymbol{f}_{*}\mid\boldsymbol{f}_{n}\sim\mathcal{N}\left(\boldsymbol{\mu}_{*}+\mathbf{K}_{*,n}\mathbf{K}_{n,n}^{-1}(\boldsymbol{f}_{n}-\boldsymbol{\mu}_{n}),\mathbf{K}_{*,*}-\mathbf{K}_{*,n}\mathbf{K}_{n,n}^{-1}\mathbf{K}_{n,*}\right).$
(20)

Despite matching asymptotic costs, iterative approaches to sampling
$\boldsymbol{f}$ are substantially slower than simultaneous ones. In applied
settings, however, test locations $\mathbf{X}_{*}$ are often determined
adaptively, forcing location-scale-based methods for generating
$\boldsymbol{f}$ to repeatedly compute (20). Further refining this
predicament, we arrive at a final challenge: pathwise derivatives.

Differentiation is a linear operation. The gradient of a Gaussian process $f$
with respect to a location $\boldsymbol{x}$ is, therefore, another Gaussian
process $f^{\prime}$. By construction, these GPs are correlated. Using
gradient information to maneuver along a sample path—for example, to identify
its extrema—therefore requires us to re-condition both processes on the
realized values of $f(\boldsymbol{x})$ and $f^{\prime}(\boldsymbol{x})$ at
each successive step of gradient descent.

Prior to continuing, it is worth noting that the limitations of location-scale
methods can be avoided in certain cases. In particular, the otherwise cubic
costs for computing a square root in $\eqref{eqn:location_scale_gp}$ can be
dramatically reduced by exploiting structural assumptions regarding covariance
matrices $\mathbf{K}$. Well-known examples of structured matrices include
banded and sparse ones in the context of one-dimensional Gaussian processes
and Gauss–Markov random fields rue2005gaussian (68, 21, 49), block-Toeplitz
Toeplitz-block ones when evaluating stationary product kernels on regularly-
spaced grids $\mathbf{X}\subset\mathcal{X}$ zimmerman1989computationally (101,
99, 19), and kernel-interpolation-based ones wilson2015kernel (96, 61). When
the task at hand permits their usage, these methods are highly effective.

The following sections survey different approaches to overcoming the
challenges put forth above by approximating Gaussian process priors as finite-
dimensional Bayesian linear models.

###  4.2 Stationary covariances

Stationary covariance functions
$k(\boldsymbol{x},\boldsymbol{x}^{\prime})=k(\boldsymbol{x}-\boldsymbol{x}^{\prime})$,
such as the Matérn family’s limiting squared exponential kernel, give rise to
a significant portion of GP priors in use today. For centered priors
$f\sim\mathcal{GP}(0,k)$, stationarity encodes the belief that the
relationship between process values $f(\boldsymbol{x}_{i})$ and
$f(\boldsymbol{x}_{j})$ is solely determined by the difference
$\boldsymbol{x}_{i}-\boldsymbol{x}_{j}$ between locations $\boldsymbol{x}_{i}$
and $\boldsymbol{x}_{j}$. Simple but expressive, stationarity is the go-to
modeling assumption in many applied settings.

These kernels exhibit a variety of special properties that greatly facilitate
the construction of efficient, approximate priors. Here, we restrict attention
to kernels admitting a spectral density $\rho$, and focus on the class of
estimators formed by discretizing the spectral representation of $k$

$\displaystyle k(\boldsymbol{x}-\boldsymbol{x}^{\prime})$
$\displaystyle=\int_{\mathbb{R}^{d}}e^{2\pi
i\boldsymbol{\omega}^{\top}(\boldsymbol{x}-\boldsymbol{x}^{\prime})}\rho(\boldsymbol{\omega})\mathop{}\\!\mathrm{d}\boldsymbol{\omega}$
$\displaystyle\rho(\boldsymbol{\omega})$
$\displaystyle=\int_{\mathbb{R}^{d}}e^{-2\pi
i\boldsymbol{\omega}^{\top}\boldsymbol{x}}k(\boldsymbol{x})\mathop{}\\!\mathrm{d}\boldsymbol{x}.$
(21)

By the _kernel trick_ scholkopf01 (71), a kernel $k$ can be written as the
inner product in a corresponding reproducing kernel Hilbert space (RKHS)
$\mathcal{H}_{k}$ equipped with a feature map
$\varphi:\mathcal{X}\rightarrow\mathcal{H}_{k}$. In many cases, this inner
product can be approximated by

$k(\boldsymbol{x},\boldsymbol{x}^{\prime})=\left\langle\varphi(\boldsymbol{x}),\varphi(\boldsymbol{x}^{\prime})\right\rangle_{\mathcal{H}_{k}}\approx\boldsymbol{\phi}(\boldsymbol{x})^{\top}\
\overline{\boldsymbol{\phi}(\boldsymbol{x}^{\prime})},$ (22)

where $\boldsymbol{\phi}:\mathcal{X}\rightarrow\mathbb{C}^{\ell}$ is some
finite-dimensional feature map and
$\overline{\boldsymbol{\phi}(\boldsymbol{x}^{\prime})}$ denotes the complex
conjugate. Based on this idea, the method of _random Fourier features_
rahimi08 (64) constructs a Monte Carlo estimate to a stationary kernel by
representing the right-hand side of (22) with $\ell$ complex exponential basis
functions $\phi_{j}(\boldsymbol{x})=\ell^{\nicefrac{{-1}}{{2}}}\exp(2\pi
i\boldsymbol{\omega}_{j}^{\top}\boldsymbol{x})$, whose parameters
$\boldsymbol{\omega}_{j}$ are sampled proportional to the corresponding
spectral density $\rho(\boldsymbol{\omega}_{j})$.222Using elementary
trigonometric identities, we may also derive a related family of basis
functions $\phi:\mathcal{X}\to\mathbb{R}^{\ell}$ with
$\phi_{j}(\boldsymbol{x})=\sqrt{\nicefrac{{2}}{{\ell}}}\cos(2\pi\boldsymbol{\omega}_{j}^{\top}\boldsymbol{x}+\tau_{j})$,
where $\tau_{j}\sim\mathcal{U}(0,2\pi)$.

Given an $\ell$-dimensional basis
$\boldsymbol{\phi}=(\phi_{1},\ldots,\phi_{\ell})$, we may now proceed to
approximate the true prior according to the Bayesian linear model

$\displaystyle\tilde{f}(\cdot)$
$\displaystyle=\sum_{i=1}^{\ell}w_{i}\phi_{i}(\cdot)$ $\displaystyle w_{i}$
$\displaystyle\sim\mathcal{N}(0,1).$ (23)

Under this approximation, $\tilde{f}$ is a random function satisfying
$\boldsymbol{\tilde{f}}_{n}\sim\mathcal{N}(\boldsymbol{0},\mathbf{\Phi}_{n}^{\vphantom{\top}}\mathbf{\Phi}_{n}^{\top})$,
where $\mathbf{\Phi}_{n}=\boldsymbol{\phi}(\mathbf{X}_{n})$ is an
$n\times\ell$ matrix of features. Per the beginning of this section, then,
$\tilde{f}$ is a Gaussian process whose covariance approximates that of $f$.

The random Fourier feature approach is particularly appealing since its
position as a Monte Carlo estimator implies that the error introduced by the
$\ell$-dimensional basis $\boldsymbol{\phi}$ decays at the _dimension-free_
rate $\ell^{\nicefrac{{-1}}{{2}}}$ sutherland15 (83). This property enables us
to balance accuracy and cost by choosing $\ell$ to suite the task at hand.

###  4.3 Karhunen–Loève expansions

While exploitation of stationarity is arguably the most common route when
constructing approximate priors, it is neither unique nor optimal. A powerful
alternative is to utilize the _Karhunen–Loève expansion_ of a Gaussian process
prior castro1986principal (8, 27).

We begin by considering the family of $\ell$-dimensional Bayesian linear
models $\tilde{f}(\cdot)=\boldsymbol{\phi}(\cdot)^{\top}\boldsymbol{w}$
consisting of orthonormal basis functions $\phi_{i}:\mathcal{X}\to\mathbb{R}$
on a compact space $\mathcal{X}$. Following standard theory
fukunaga2013introduction (27), the _optimal_ $\tilde{f}$ for approximating a
Gaussian process $f$ (in the sense of minimizing mean square error) is found
by truncating its Karhunen–Loève expansion

$\displaystyle f(\cdot)$
$\displaystyle=\sum_{i=1}^{\infty}w_{i}\phi_{i}(\cdot)$ $\displaystyle w_{i}$
$\displaystyle\sim\mathcal{N}(0,\lambda_{i}),$ (24)

where $\phi_{i}$ and $\lambda_{i}$ are, respectively, the $i$-th eigenfunction
and eigenvalue of the covariance operator
$\psi\mapsto\int_{\mathcal{X}}\psi(\boldsymbol{x})k(\boldsymbol{x},\cdot)\mathop{}\\!\mathrm{d}\boldsymbol{x}$,
written in decreasing order of $\lambda_{i}$.333These eigenvalues are well-
ordered and countable as consequence of the compactness of $\mathcal{X}$.
Truncated versions of these expansions are used as both bases for constructing
optimal approximate GPs zhu1997gaussian (100, 80) and modeling tools in their
own right krainski19 (43). Depending on the case, eigenfunctions $\phi_{i}$
are either derived from first principles krainski19 (43) or obtained by
numerical methods lindgren11 (47, 50, 79).

In addition to being optimal, Karhunen–Loève expansions are exceedingly
general. Even when a covariance function $k$ is non-stationary or the domain
$\mathcal{X}$ is non-Euclidean—such as when Gaussian processes are used to
represent functions on manifolds borovitskiy2020matern (5) and graphs
borovitskiy2020graph (4)—the Karhunen–Loève expansion often exists.

Widespread use of truncated eigensystems is largely impeded by their frequent
lack of convenient, analytic forms. This issue is compounded by the fact that
efficient, numerical methods for obtaining (24) typically require us to
manipulate bespoke mathematical properties of specific kernels. These
properties are often closely related to the differential-equation-based
perspectives of Gaussian processes introduced in the following section.

###  4.4 Stochastic partial differential equations

Many Gaussian process priors, such as the Matérn family, can be expressed as
solutions of _stochastic partial differential equations_ (SPDEs). SPDEs are
common in fields such as physics, where they describe natural phenomena (such
as diffusion and heat transfer); many of which share a deep connection with
the squared exponential kernel grigoryan2009heat (29). Additionally, SPDEs are
often the starting point when designing non-stationary GP priors krainski19
(43). Below, we detail how the _Galerkin finite element method_ evans10 (24,
47, 50) can be used to construct Bayesian linear models that approximate GP
priors capable of being represented as SPDEs.

Suppose a Gaussian process $f\sim\mathcal{GP}(0,k)$ satisfies
$\mathcal{L}f=\mathcal{W}$, where $\mathcal{L}$ is a linear differential
operator and $\mathcal{W}$ is a Gaussian white noise process
lifshits2012lectures (46). Here, we demonstrate how to derive a Gaussian
process $\tilde{f}$ that approximately satisfies this SPDE. To begin, we
express $\mathcal{L}f=\mathcal{W}$ in its weak form444One typically integrates
$(\mathcal{L}f)(\boldsymbol{x})g(\boldsymbol{x})$ by parts, either by
necessity or due to affordances of the basis $\phi_{i}$. We suppress this to
ease notation.

$\int_{\mathcal{X}}(\mathcal{L}f)(\boldsymbol{x})g(\boldsymbol{x})\mathop{}\\!\mathrm{d}\boldsymbol{x}=\int_{\mathcal{X}}g(\boldsymbol{x})\mathop{}\\!\mathrm{d}\mathcal{W}(\boldsymbol{x}),$
(25)

where $g$ is an arbitrary element of an appropriate class of test functions.
Next, we proceed by approximating both the desired solution $f$ and the test
function $g$ with respect to a finite-dimensional basis as
$\tilde{f}(\cdot)=\sum_{i=1}^{\ell}w_{i}\phi_{i}(\cdot)$ and
$\tilde{g}(\cdot)=\sum_{j=1}^{\ell}v_{j}\phi_{j}(\cdot)$. Substituting these
terms into (25) and differentiating both sides with respect to the
coefficients of $\tilde{g}$, we obtain the following expression for each
$j=1,\ldots,\ell$

$\sum_{i=1}^{\ell}w_{i}\underbracket{\int_{\mathcal{X}}(\mathcal{L}\phi_{i})(\boldsymbol{x})\phi_{j}(\boldsymbol{x})\mathop{}\\!\mathrm{d}\boldsymbol{x}}_{A_{ij}}=\underbracket{\int_{\mathcal{X}}\phi_{j}(\boldsymbol{x})\mathop{}\\!\mathrm{d}\mathcal{W}(\boldsymbol{x})}_{b_{j}}.$
(26)

Defining $\mathbf{M}=\operatorname{Cov}(\boldsymbol{b})$, where
$\operatorname{Cov}(b_{i},b_{j})=\left\langle\phi_{i},\phi_{j}\right\rangle$
coincides with the finite-element mass matrix, allows us to rearrange this
system of random linear equations in matrix-vector form by writing
$\mathbf{A}\boldsymbol{w}=\boldsymbol{b}$. The basis coefficients of the
random function $\tilde{f}$ are, therefore, distributed as
$\boldsymbol{w}\sim\mathcal{N}\left(\boldsymbol{0},\mathbf{A}^{-1}\mathbf{M}\mathbf{A}^{-\top}\right)$.
As in the previous sections, $\tilde{f}$ can be seen as the weight-space view
of a corresponding Gaussian process.

A popular choice is to employ compactly supported basis functions $\phi_{i}$
lindgren11 (47). The matrices $\mathbf{A}$ and $\mathbf{M}$ are then sparse,
and the resulting linear systems can be solved efficiently. For example, the
family of piecewise linear basis functions is a simple but effective choice
for second order differential operators $\mathcal{L}$ evans10 (24, 50).555A
second order differential operator gives rise to a first-order bilinear form
when integrated by parts, which matches with piecewise linear basis functions
which are once differentiable almost everywhere. For higher-order operators, a
piecewise polynomial basis may be used instead.

###  4.5 Discussion

This section has focused on identifying finite-dimensional bases with which to
construct Bayesian linear models
$\tilde{f}(\cdot)=\boldsymbol{\phi}(\cdot)^{\top}\boldsymbol{w}$. These model
can be seen as _weight-space_ interpretations rasmussen06 (67) of
corresponding Gaussian process priors $\tilde{f}\sim\mathcal{GP}(0,\tilde{k})$
with covariance functions
$\tilde{k}(\boldsymbol{x},\boldsymbol{x}^{\prime})=\boldsymbol{\phi}(\boldsymbol{x})^{\top}\mathbf{\Sigma}_{\boldsymbol{w}}\boldsymbol{\phi}(\boldsymbol{x}^{\prime})$.
Since $\boldsymbol{w}$ and
$\boldsymbol{\tilde{f}}_{n}=\mathbf{\Phi}_{n}\boldsymbol{w}$ are jointly
normal, Theorem 1 implies that we may enforce the condition
$\boldsymbol{\tilde{f}}_{n}=\boldsymbol{y}$ by writing666Practical variants of
(27) avoid inverting
$\mathbf{\Phi}_{n}^{\vphantom{\top}}\mathbf{\Phi}_{n}^{\top}$ by employing,
e.g., Gaussian likelihoods (Section 5.1).

$\boldsymbol{\phi}(\cdot)^{\top}(\boldsymbol{w}\mid\boldsymbol{y})\overset{\mathop{}\\!\mathrm{d}}{=}\boldsymbol{\phi}(\cdot)^{\top}\left(\boldsymbol{w}+\mathbf{\Phi}_{n}^{\top}(\mathbf{\Phi}_{n}^{\vphantom{\top}}\mathbf{\Phi}_{n}^{\top})^{-1}(\boldsymbol{y}-\mathbf{\Phi}_{n}\boldsymbol{w})\right).$
(27)

This result encourages us to approximate posteriors in much the same way as we
have priors. After all, if we have chosen a basis $\boldsymbol{\phi}$ that
encodes our prior knowledge for $f$ (such as how smooth we believe this
function to be), then it is reasonable to think that $\boldsymbol{\phi}$ will
further enable us to efficiently approximate $f\mid\boldsymbol{y}$. To the
extent that this approach may seem like the natural evolution of ideas
discussed in this section, we argue for the benefits of _decoupling_ the
representation of the prior from that of the data.

The trouble with using a finite set of homogeneous basis functions
$\boldsymbol{\phi}=(\phi_{1},\ldots,\phi_{\ell})$ to represent both the prior
and the data is that these two tasks focus on different things. To accurately
approximate a prior is to faithfully describe a random function $f$ on a
domain $\mathcal{X}$. Consequently, parsimonious approximations $\tilde{f}$
employ global basis functions that vary non-trivially everywhere on
$\mathcal{X}$. This is largely why, e.g., Fourier features are an attractive
choice for approximating stationary priors. But what of the data?

Conditioning on observations $\boldsymbol{y}$ requires us to convey how our
understanding of $f$ has changed. In most cases, we choose priors (and
likelihoods) that reflect the belief that an observation $y_{i}$ only informs
us about the process $f$ in the immediate vicinity of a point
$\boldsymbol{x}_{i}$. Updating $f$ to account for $\boldsymbol{y}$, therefore,
typically focuses on process values corresponding to specific regions of
$\mathcal{X}$. Rather than global basis functions, the data is best
characterized by local ones that have near-zero values outside of the
aforementioned regions. Not coincidentally, the canonical basis functions
$k(\cdot,\boldsymbol{x})$ fit this description perfectly when the chosen prior
implies that $y_{i}$ is only locally informative.

A key property of pathwise conditioning is that it not only provides us with a
natural decomposition of GP posteriors—as sums of prior random variables and
data-driven updates—but enables us to represent these terms in separate bases.
Similar ideas can be found in recent works that explore alternative
decompositions of Gaussian processes, such as separation of mean and
covariance functions cheng17 (11, 69) or decoupling of RKHS subspaces and
their orthogonal complements shi19 (74). Unlike these works, however, we
stress decoupling in the sense of using different classes of basis functions
to represent different aspects of GP posteriors. While this type of decoupling
is not unique to pathwise approaches lazaro2009inter (44, 31), they
drastically simplify the process by eliminating the need to analytically solve
for sufficient statistics.

Figure 3: Overview of variance starvation when conditioning on
$n\in\\{10,100,1000\\}$ observations of the form
$y_{i}\sim\mathcal{N}(f_{i},10^{-5})$ located within the gray shaded region.
_Top:_ Comparison of pathwise updates to a single draw from an approximate
prior $\tilde{f}(\cdot)=\boldsymbol{\phi}(\cdot)^{\top}\boldsymbol{w}$,
constructed using $\ell=1000$ Fourier features $\phi$. Updates defined using
the same Fourier basis $\boldsymbol{\phi}(\cdot)$ and the canonical basis
functions $k(\cdot,\mathbf{X})$ are shown in blue and dashed-black,
respectively. _Bottom:_ Mean and two standard deviations of the empirical
posteriors formed by applying the aforementioned updates to $10^{5}$ draws
from the approximate prior.

This line of reasoning also helps to explain why finite-dimensional GPs
constructed from homogeneous basis functions often produce poorly-calibrated
posteriors. For now, we restrict our attention to the issue of _variance
starvation_ wang18 (94, 55, 7) and return this topic in Section 5.5. Figure 3
demonstrates what happens as the number of observations $n=|\boldsymbol{y}|$
approaches the number of random Fourier features $\ell=1000$ used to
approximate a squared exponential kernel. In general, the approximate
posteriors produce extrapolations which become increasingly erratic. Note that
the rate at which these defects materialize depends upon the choice of kernel
and likelihood. In the figure, posteriors yielded by pathwise updates in
canonical and Fourier bases (all other things being held equal) diverge as the
number of observations $n$ approaches the number of random Fourier features
$\ell$. This pattern emerges because the Fourier basis is better at describing
stationary priors than non-stationary posteriors. Fourier features excel at
capturing the global properties of the prior, but struggle to portray the
localized effects of the data.

Of course, different types of data impose different kinds of conditions on the
process $f$. We now examine various pathwise updates that enforce prominent
types of conditions.

##  5 Conditioning via pathwise updates

Building off of the foundation prepared in Section 3, we now adapt Corollary 4
to accommodate different types of conditions and computational budgets.
Throughout this section, we use $\boldsymbol{\gamma}$ to denote the random
variable realized by observations $\boldsymbol{y}$ under the chosen
likelihood.

###  5.1 Gaussian updates

Corollary 4 treats observations $\boldsymbol{y}$ as a realization of process
values $\boldsymbol{f}_{n}=f(\mathbf{X}_{n})$. Hence, the conditions it
imposes manifest as the equality constraint
$\boldsymbol{f}_{n}=\boldsymbol{y}$. In the real world, however, we seldom
observe $\boldsymbol{f}_{n}$ directly. To account for this nuance, an
observation $y$ is modeled by a _likelihood_ $p\left(y\mid
f(\boldsymbol{x})\right)$. Viewed from this perspective, the equality
constraint $\boldsymbol{f}_{n}=\boldsymbol{y}$ correspond to the limit where
$p$ contracts to a point mass. Seeing as $y$ usually fails to fully
disambiguate the true value of $f(\boldsymbol{x})$, we typically employ
likelihoods that induce weaker conditions than strict equalities.

For regression problems, the most common choice is to employ a Gaussian
likelihood $p(y\mid f(\boldsymbol{x}))=\mathcal{N}(y\mid
f(\boldsymbol{x}),\sigma^{2})$, the log of which penalizes the squared
Euclidean distance of $f(\boldsymbol{x})$ from $y$. Under the corresponding
observation model
$\boldsymbol{\gamma}=\boldsymbol{f}_{n}+\boldsymbol{\varepsilon}$ with
$\boldsymbol{\varepsilon}\sim\mathcal{N}(\boldsymbol{0},\sigma^{2}\mathbf{I})$,
$f$ and $\boldsymbol{y}$ are jointly Gaussian. By Corollary 4 then, we may
condition $f$ on $\boldsymbol{\gamma}=\boldsymbol{y}$ by writing

$\displaystyle\begin{split}(f\mid\boldsymbol{\gamma}=\boldsymbol{y})(\cdot)&\overset{\mathop{}\\!\mathrm{d}}{=}f(\cdot)+k(\cdot,\mathbf{X})(\mathbf{K}_{n,n}+\sigma^{2}\mathbf{I})^{-1}(\boldsymbol{y}-\boldsymbol{f}_{n}-\boldsymbol{\varepsilon}).\end{split}$
(28)

Rather than exactly passing through observations $\boldsymbol{y}$, the
conditioned path $f\mid\boldsymbol{y}$ now smoothly interpolates between them.
In cases where $\boldsymbol{\gamma}$ is not a Gaussian random variable,
additional tools are needed.

Figure 4: Visual comparison of different pathwise updates. _Left and middle:_
Variational inference is used to learn sparse updates at $m=10$ inducing
locations $\mathbf{Z}$ (circles). _Right:_ preconditioned conjugate gradients
is used to iteratively solve for Gaussian updates. In all cases, 1000
observations $\boldsymbol{y}$ are evenly spaced in the shaded region. Dashed
lines denote mean and two standard deviations of ground truth posteriors,
colored regions and thicker lines denote those of empirical ones. Middle and
right plots illustrate regression with a Gaussian likelihood
$\mathcal{N}(y_{i}\mid f_{i},10^{-3})$. The left plot shows binary
classification with a Bernoulli likelihood and probit link function $g$; every
tenth label is shown as a small, vertical bar.

###  5.2 Non-Gaussian updates

In the general setting, where the random variable $\boldsymbol{\gamma}$ is
arbitrarily distributed under the chosen likelihood, $\boldsymbol{\gamma}$
relates to process values $\boldsymbol{f}$ by way of the non-conjugate prior

$p(\boldsymbol{\gamma},\boldsymbol{f})=p\left(\boldsymbol{\gamma}\mid
g^{-1}(\boldsymbol{f})\right)\mathcal{N}(\boldsymbol{f}\mid\boldsymbol{\mu},\mathbf{K}),$
(29)

where the _link function_ $g:\mathcal{Y}\to\mathbb{R}$ maps from the space of
predictions $\mathcal{Y}\subset\mathbb{R}$ to the range of $f$. For binary
classification problems, popular choices for $g:[0,1]\to\mathbb{R}$ include
logit and probit functions rasmussen06 (67). The left column of Figure 4
illustrates this scenario using methods described below.

Even under a non-conjugate prior (29), the conditional expectation
$\operatorname{\mathbb{E}}(f\mid\boldsymbol{\gamma})$ and the residual
$f-\operatorname{\mathbb{E}}(f\mid\boldsymbol{\gamma})$ it induces are
uncorrelated (see Section 2.3). Since $p(\boldsymbol{f},\boldsymbol{\gamma})$
may not be Gaussian, however, it no longer follows that this lack of
correlation implies independence—hence, the pathwise update (16) may not hold.

Exact Bayesian inference and prediction are typically intractable when dealing
with non-conjugate priors. Strategies for circumventing this issue generally
approximate the true posterior by introducing an auxiliary random variable
$\boldsymbol{u}\sim q(\boldsymbol{u})$ such that $f\mid\boldsymbol{u}$
resembles $f\mid\boldsymbol{y}$ according to a chosen measure of similarity
nickisch2008approximations (57, 33). For practical reasons, $\boldsymbol{u}$
is typically assumed to be jointly Gaussian with $\boldsymbol{f}$.777Note
that, in the special case where $p(\boldsymbol{f},\boldsymbol{\gamma})$ is
Gaussian, the optimal $q$ is also Gaussian titsias09b (88). Consequently, non-
conjugate priors $p(\boldsymbol{f},\boldsymbol{\gamma})$ are replaced by
conjugate ones $p(\boldsymbol{f},\boldsymbol{u})$ to aid in the construction
of approximate posteriors, whereupon Matheron’s update rule holds once more.
The following section explores these _sparse_ approximations in greater
detail.

###  5.3 Sparse updates

Approximations to GP posteriors frequently revolve around conditioning a
process $f$ on a random variable
$\boldsymbol{u}=(u_{1},\ldots,u_{m})\in\mathbb{R}^{m}$. Per the previous
section, this may be because the outcome variable $\boldsymbol{\gamma}$ is
non-Gaussian nickisch2008approximations (57, 86, 33). Alternatively, the
$\mathcal{O}(n^{3})$ cost for directly conditioning on all
$n=|\boldsymbol{y}|$ observations may be prohibitive titsias09a (87, 32). In
these cases and more, we would like to infer a distribution
$q(\boldsymbol{u})$ such that $f\mid\boldsymbol{u}$ explains the data.
Defining (approximate) posteriors in this way not only avoids potential issues
arising from non-Gaussianity of $\boldsymbol{\gamma}$, but associates the
computational cost of conditioning with $\boldsymbol{u}$. As discussed below,
this leads to pathwise updates that run in $\mathcal{O}(m^{3})$ time.

Comprehensive treatment of different approaches to learning _inducing
distributions_ $q(\boldsymbol{u})$ is beyond the scope of this work. In
general, however, these procedures operate by finding an approximate posterior
$q(\boldsymbol{f},\boldsymbol{u})$ within a tractable family of approximating
distributions $\mathcal{Q}$. For reasons that will soon become clear, this
family of distributions typically includes an additional set of parameters
$\mathbf{Z}$, which help to define the joint distribution
$p(\boldsymbol{f},\boldsymbol{u})$. To help streamline presentation, we focus
on the simplest and most widely used abstraction for inducing variables
$\boldsymbol{u}$: namely, _pseudo-data_.

The noise-free pseudo-data framework snelson06 (77, 63, 87) treats each draw
of a random vector $\boldsymbol{u}\sim q(\boldsymbol{u})$ as a realization of
process values $\boldsymbol{f}_{m}=f(\mathbf{Z})$ at a corresponding set of
tunable locations $\mathbf{Z}\in\mathcal{X}^{m}$. This paradigm gets its name
from the intuition that the (random) collection of pseudo-data
$\left(\boldsymbol{z}_{j},u_{j}\right)_{j=1}^{m}$ mimics the effect of a
noise-free data set $\left(\boldsymbol{x}_{i},f_{i}\right)_{i=1}^{n}$ on $f$.
By construction, $\boldsymbol{u}$ is jointly Gaussian with $f$.888This
condition holds when $\boldsymbol{u}$ relates to $f$ via a linear map
lazaro2009inter (44). Appealing to Corollary 4, we define the _sparse pathwise
update_ as

$(f\mid\boldsymbol{u})(\cdot)\overset{\mathop{}\\!\mathrm{d}}{=}f(\cdot)+\underbracket{\sum_{i=1}^{m}v_{i}k(\cdot,\boldsymbol{z}_{i})}_{\mathclap{m\mathrel{\text{missing}}{-dimensionalbasis}}},$
(30)

where
$\boldsymbol{v}=\mathbf{K}_{m,m}^{-1}\left(\boldsymbol{u}-\boldsymbol{f}_{m}\right)$.
This formula is identical to the one given by Corollary 4, save for the fact
that we now sample $\boldsymbol{u}\sim q(\boldsymbol{u})$ and solve for a
linear system involving the $m\times m$ covariance matrix
$\mathbf{K}_{m,m}=k(\mathbf{Z},\mathbf{Z})$ at $\mathcal{O}(m^{3})$ cost. The
middle column of Figure 4 illustrates the sparse update induced by Gaussian
$\boldsymbol{u}\sim\mathcal{N}(\boldsymbol{\mu}_{\boldsymbol{u}},\mathbf{\Sigma}_{\boldsymbol{u}})$
with learned moments $\boldsymbol{\mu}_{\boldsymbol{u}}$ and
$\mathbf{\Sigma}_{\boldsymbol{u}}$.

Just as we can imitate process values $\boldsymbol{f}_{n}$, we can also
emulate (Gaussian) observations $\boldsymbol{y}$. This intuition leads to the
Gaussian pseudo-data family of inducing distributions, whose moments

$\displaystyle\boldsymbol{\mu}_{u}$
$\displaystyle=\mathbf{K}_{m,m}(\mathbf{K}_{m,m}+\mathbf{\Lambda})^{-1}\boldsymbol{\tilde{y}}$
$\displaystyle\mathbf{\Sigma}_{\boldsymbol{u}}$
$\displaystyle=(\mathbf{K}_{m,m}^{-1}+\mathbf{\Lambda})^{-1}$ (31)

are parameterized by _pseudo-observations_
$\boldsymbol{\tilde{y}}\in\mathbb{R}^{m}$ and _pseudo-noise_
$\boldsymbol{\tilde{\sigma}}\in\mathbb{R}_{+}^{m}$, where
$\mathbf{\Lambda}=\operatorname{diag}(\boldsymbol{\tilde{\sigma}}^{2})$. This
choice of parameterization is motivated by the observation that, given $n\leq
m$ Gaussian random variables
$\boldsymbol{\gamma}\sim\mathcal{N}(\boldsymbol{f}_{n},\sigma^{2}\mathbf{I})$,
the family of distributions it generates contains the optimal $q$ despite
housing only $\mathcal{O}(m)$ free terms seeger1999bayesian (72, 59).999We
recover the true posterior by, e.g., taking
$(\tilde{y}_{i},\tilde{\sigma}_{i})=(y_{i},\sigma)$ for all $i\leq n$ and
sending $\tilde{\sigma}_{i}\to\infty$ otherwise. Using the Gaussian pathwise
update (28), we may express $\boldsymbol{u}$ itself as

$\displaystyle\boldsymbol{u}$
$\displaystyle\overset{\mathop{}\\!\mathrm{d}}{=}\boldsymbol{f}_{m}+\mathbf{K}_{m,m}(\mathbf{K}_{m,m}+\mathbf{\Lambda})^{-1}(\boldsymbol{\tilde{y}}-\boldsymbol{f}_{m}-\boldsymbol{\tilde{\varepsilon}})$
$\displaystyle\boldsymbol{\tilde{\varepsilon}}$
$\displaystyle\sim\mathcal{N}(\boldsymbol{0},\mathbf{\Lambda}).$ (32)

Here, despite the fact that $\boldsymbol{f}_{m}$ and
$\boldsymbol{\tilde{\varepsilon}}$ generate $\boldsymbol{u}$, it remains the
case that
$\operatorname{Cov}(\boldsymbol{f}_{m}+\boldsymbol{\tilde{\varepsilon}},\boldsymbol{u})=\mathbf{0}$.
Substituting this expression into (30) and simplifying gives the pathwise
update101010This same line of reasoning leads to a _rank-1 pathwise update_
for cases where conditions arrive online.

$(f\mid\boldsymbol{u})(\cdot)\overset{\mathop{}\\!\mathrm{d}}{=}f(\cdot)+k(\cdot,\mathbf{Z})\left(\mathbf{K}_{m,m}+\mathbf{\Lambda}\right)^{-1}\left(\boldsymbol{\tilde{y}}-\boldsymbol{f}_{m}-\boldsymbol{\tilde{\varepsilon}}\right).$
(33)

Hence, while sampling $\boldsymbol{u}$ is more complicated in the Gaussian
pseudo-data case, the resulting pathwise update is straightforward. This
family of inducing distributions is particularly advantageous in the large $m$
setting, both because it contains only $\mathcal{O}(m)$ free parameters and
for reasons discussed in the following section.

In rough analogy to methods discussed in Section 4, we may think of the sparse
updates introduced here as using an $m$-dimensional basis
$k(\cdot,\mathbf{Z})$ to approximate functions defined in terms of the
$n$-dimensional basis $k(\cdot,\mathbf{X}_{n})$. In practice, this basis is
often efficient because neighboring training locations give rise to similar
basis functions. Kernel basis functions at appropriately chosen sets of $m\ll
n$ locations $\mathbf{Z}$ exploit this redundancy to produce a sparser, more
cost-efficient representation. burt20 (6) study this problem in detail and
derive bounds on the quality of variational approximations to GP posteriors as
$m\to n$.

###  5.4 Iterative solvers

Throughout this section, we have focused on the high-level properties of
pathwise updates in relation to various problem settings. We have said little,
however, regarding the explicit means of executing such an update. In all
cases discussed here, pathwise updates have amounted to solutions to system of
linear equations. For example, the update originally featured in Corollary 4
solves the system
$\mathbf{K}_{n,n}\boldsymbol{v}=\boldsymbol{y}-\boldsymbol{f}_{n}$ for a
vector of coefficients $\boldsymbol{v}$, which define how the same realization
of $f$ changes when subjected to the condition
$\boldsymbol{f}_{n}=\boldsymbol{y}$. Given a reasonable number of conditions
$n$ (up to several thousand), we may obtain $\boldsymbol{v}$ by first
computing the Cholesky factor
$\mathbf{L}_{n,n}=\mathbf{K}^{\nicefrac{{1}}{{2}}}_{n,n}$ and then solving for
a pair of triangular systems
$\mathbf{L}_{n,n}\boldsymbol{\bar{v}}=\boldsymbol{u}-\boldsymbol{f}_{n}$ and
$\mathbf{L}_{n,n}^{\top}\boldsymbol{v}=\boldsymbol{\bar{v}}$. For large $n$,
however, the $\mathcal{O}(n^{3})$ time complexity for carrying out this recipe
is typically prohibitive.

Rather than solving for coefficients $\boldsymbol{v}$ directly, we may instead
employ an _iterative solver_ that constructs a sequence of estimates
$\boldsymbol{v}^{(1)},\boldsymbol{v}^{(2)},\dotsc$ to $\boldsymbol{v}$, such
that $\boldsymbol{v}^{(j)}$ converges to the true $\boldsymbol{v}$ as $j$
increases. Depending on the numerical properties of the linear system in
question, it is possible (or even likely) that a high-quality estimate
$\boldsymbol{v}^{(j)}$ will be obtained after only $j\ll n$ iterations. This
line of reasoning features prominently in a number of recent works, where
iterative solvers have been shown to be highly competitive for purposes of
approximating GP posteriors pleiss2018constant (61, 28, 93). The right column
of Figure 4 visualizes an iterative solution to the Gaussian pathwise update
(28) obtained using preconditioned conjugate gradients gardner2018gpytorch
(28).

In these cases, posterior sampling via pathwise conditioning enjoys an
important advantage over distributional approaches: it allows us to solve for
linear system of the form $\mathbf{K}^{-1}_{n,n}\boldsymbol{v}$ rather than
working with $\mathbf{K}^{\nicefrac{{1}}{{2}}}_{*,*\mid n}\boldsymbol{\zeta}$.
Whereas the former amounts to a standard solve, the latter often requires
special considerations pleiss2020fast (62) and can be difficult to work with
when typical square root decompositions prove impractical parker12 (60).

Lastly, we note that these techniques can be combined with sparse
approximations for improved scaling in $m$ and faster convergence of iterative
solves. As a concrete example, we return to the Gaussian pseudo-data
variational family (31). By construction, the corresponding pathwise update
(33) closely resembles the original Gaussian update (28). In general, however,
pseudo-noise variances $\tilde{\sigma}_{i}^{2}$ are often significantly larger
than the true noise variance $\sigma^{2}$. The resulting linear system
$(\mathbf{K}_{m,m}+\mathbf{\Lambda})^{-1}\boldsymbol{v}$ is, therefore,
substantially better-conditioned than that of the exact alternative—implying
that it can be solved in far fewer iterations.

###  5.5 Discussion

In Section 4.5, we discussed finite-dimensional approximations of Gaussian
process posteriors. There, we explored how the globality of the prior
reinforces the use of basis functions $\phi_{i}:\mathcal{X}\to\mathbb{R}$ that
inform us about $f$ over the entire domain $\mathcal{X}$, while the localized
effects of the data encourages the use of $\phi_{i}$ that only tell us about
$f$ on subsets of $\mathcal{X}$. This conflict hinders our ability to
efficiently represent both the prior and the data (i.e., the posterior) using
a single class of basis functions. That discussion ended with a demonstration
of what happens when $\boldsymbol{\phi}=(\phi_{1},\ldots,\phi_{\ell})$ solely
consists of global basis functions, specifically random Fourier features. Most
works, however, have focused on the use of canonical basis functions
$k(\cdot,\boldsymbol{x})$, which are typically local. This section, therefore,
aims to fill in the gaps.

At the end of Section 4.5, we saw how trouble conveying the data in global
bases led to approximate posteriors that were starved for variance (Figure 3).
Writing the update rules—for a draw from an approximate prior
$\tilde{f}(\cdot)=\boldsymbol{\phi}(\cdot)^{\top}\boldsymbol{w}$ subject to
the condition $\boldsymbol{\tilde{f}}_{n}=\boldsymbol{y}$—in both unified and
decoupled bases side-by-side helps to highlight their key differences

$\displaystyle\underbracket{\tilde{f}(\cdot)+\boldsymbol{\phi}(\cdot)^{\top}\mathbf{\Phi}_{n}^{\top}\big{(}\mathbf{\Phi}_{n}^{\vphantom{\top}}\mathbf{\Phi}_{n}^{\top}\big{)}^{-1}\big{(}\boldsymbol{y}-\boldsymbol{\tilde{f}}_{n}\big{)}\vphantom{\tilde{f}(\cdot)+k(\cdot,\mathbf{X}_{n})\mathbf{K}_{n,n}^{-1}(\boldsymbol{y}-\boldsymbol{\tilde{f}}_{n})}}_{\mathrel{\text{missing}}{unifiedapproximateposterior}}$
$\displaystyle\underbracket{\tilde{f}(\cdot)+k(\cdot,\mathbf{X}_{n})\mathbf{K}_{n,n}^{-1}\big{(}\boldsymbol{y}-\boldsymbol{\tilde{f}}_{n}\big{)}}_{\mathrel{\text{missing}}{decoupledapproximateposterior}}.$
(34)

On the right, the cross-covariance term
$\boldsymbol{\phi}(\cdot)^{\top}\mathbf{\Phi}_{n}^{\top}=\boldsymbol{\phi}(\cdot)^{\top}\boldsymbol{\phi}(\mathbf{X}_{n})$
is replaced by $k(\cdot,\mathbf{X}_{n})$. Seeing as the former is often chosen
to approximate the latter in a way that converges when an appropriate limit is
taken, for instance in (22), it comes as no surprise that
$k(\cdot,\mathbf{X}_{n})$ more accurately represents data. Moreover, the
matrix inverse
$\big{(}\mathbf{\Phi}_{n}^{\vphantom{\top}}\mathbf{\Phi}_{n}^{\top}\big{)}^{-1}$
appearing on the left is often ill-conditioned and, therefore, amplifies
numerical errors. Finite-dimensional GPs constructed from local basis
functions exhibit similar issues, albeit for essentially the opposite reason.
Rather than failing to adequately represent the data, local basis functions
struggle to reproduce the prior.

Many approaches to approximating Gaussian processes $f\sim\mathcal{GP}(0,k)$
revolve around representing the data in terms of $m$-dimensional canonical
bases $k(\cdot,\mathbf{Z})$; for a review, see quinonero2007approximation
(63). Early iterations of this strategy silverman1985some (76, 92, 85),
typically used $k(\cdot,\mathbf{Z})$ to define degenerate Gaussian processes
rasmussen06 (67). Here, the term _degenerate_ emphasizes the fact that the
covariance function

$\tilde{k}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})=k(\boldsymbol{x}_{i},\mathbf{Z})k(\mathbf{Z},\mathbf{Z})^{-1}k(\mathbf{Z},\boldsymbol{x}_{j})$
(35)

of such a process has a finite number of non-zero eigenvalues. From the
weight-space perspective, degenerate GPs are Bayesian linear models
$\tilde{f}(\cdot)=k(\cdot,\mathbf{Z})\boldsymbol{w}$, which makes it clear
that $\tilde{f}(\cdot)$ goes to zero as
$k(\cdot,\mathbf{Z})\to\boldsymbol{0}$. This behavior is particularly
troublesome if all $\boldsymbol{z}\in\mathbf{Z}$ are positioned near training
locations $\mathbf{X}_{n}$: since $k(\boldsymbol{x}_{*},\mathbf{Z})$ typically
vanishes as $\boldsymbol{x}_{*}$ retreats from $\mathbf{Z}$, both the prior
and the posterior collapse to point masses away from the data.

Instead of focusing on the data, one idea is to start by finding a basis
$k(\cdot,\mathbf{Z})$ capable of accurately reproducing the prior.
Accomplishing this feat will require us to use a relatively large number of
basis functions, since $\mathbf{Z}$ will need to effectively cover the
(compact) domain $\mathcal{X}$. As mentioned in Section 4.1, certain kernels
produce special kinds of matrices when evaluated on particular sets.
Exploiting these special properties—e.g., by taking the Toeplitz matrices
formed when evaluating a stationary product kernel $k$ on a regularly spaced
grid $\mathbf{Z}$ and embedding them inside of circulant ones
wood1994simulation (99, 19)—enables us to drastically reduce the cost of
expensive matrix operations, such as multiplies, decompositions, and inverses.
Especially when $\mathcal{X}$ is low dimensional, then, we can use the
canonical basis to efficiently approximate the prior.

Kernel interpolation methods wilson2015kernel (96, 61) take this idea a step
further. Given a set of $m$ inducing locations $\mathbf{Z}$, let
$\boldsymbol{\xi}:\mathcal{X}\to\mathbb{R}^{m}$ be a _weight function_
silverman1984spline (75) mapping locations $\boldsymbol{x}_{i}$ onto (sparse)
weight vectors $\boldsymbol{\xi}_{i}$ such that
$k(\boldsymbol{x}_{i},\mathbf{Z})\approx\boldsymbol{\xi}_{i}^{\top}k(\mathbf{Z},\mathbf{Z})$.
By applying this technique to (35), we can define another Gaussian process
$g\sim\mathcal{GP}(0,c)$ with degenerate covariance
$c(\boldsymbol{x}_{i},\boldsymbol{x}_{j})=\boldsymbol{\xi}_{i}^{\top}k(\mathbf{Z},\mathbf{Z})\boldsymbol{\xi}_{j}^{\vphantom{\top}}$.
As a Bayesian linear model, we have
$g(\cdot)=\boldsymbol{\xi}(\cdot)^{\top}\boldsymbol{g}_{m}$. Notice that
process values $\boldsymbol{g}_{m}=g(\mathbf{Z})$ now play the role of random
weights $\boldsymbol{w}$ and fully determine the behavior of the random
function $g$. Assuming $\mathbf{Z}$ was chosen so that
$k(\mathbf{Z},\mathbf{Z})$ admits convenient structure, random vectors
$\boldsymbol{g}_{m}\mid\boldsymbol{y}$ and, hence, random functions
$(g\mid\boldsymbol{y})(\cdot)$ can be obtained cheaply pleiss2018constant
(61). When $\mathbf{Z}$ is sufficiently dense in $\mathcal{X}$ (so as to be
reasonably close to $\boldsymbol{x}_{*}$), this strategy provides an
alternative means of efficiently sampling from GP posteriors.

###  5.6 An empirical study

Figure 5:  Accuracy and cost of different methods for sampling from GP
posteriors given $n$ observations. Draws from the prior are generated using
either location-scale (_top_) or $\ell=4096$ random Fourier features
(_bottom_). We denote Gaussian updates by black dots, sparse updates by blue
stars, CG updates by orange and red triangles, and RFF updates by green
diamonds. Sparse and RFF updates both utilized $m=\frac{n}{4}$ basis
functions. All results are reported as medians and interquartile ranges
measured over 32 independent trials. _Left:_ 2-Wasserstein distances of
empirical distributions of $10^{5}$ samples from the ground truth GP
posterior. _Middle and right:_ Time taken to generate a draw of
$(\boldsymbol{f}_{*}\mid\cdot\,)\in\mathbb{R}^{1024}$ with and without caching
of terms that are independent of $\mathbf{X}_{*}$.

By now, we have explored a variety of techniques for sampling from GP
posteriors. Each of these methods is well suited for a particular type of
problem. To help shed light on their respective niches, we conducted a simple
controlled experiment.

Here, our goal is to better understand how different methods balance the
tradeoff of cost and accuracy. We measured cost in terms of runtimes and
accuracy in terms of 2-Wasserstein distances between empirical distributions
and true posterior (see Section 6). To eliminate confounding variables, we
assumed a known Matérn-$\nicefrac{{5}}{{2}}$ prior on random functions
$f:\mathbb{R}^{4}\to\mathbb{R}$. All trials began by sampling this prior at
$n$ training locations $\mathbf{X}_{n}$ and 1024 test locations
$\mathbf{X}_{*}$, using either location-scale transforms or random Fourier
features. We then used the various update rules explored in this section to
condition on $n$ observations
$\boldsymbol{y}\sim\mathcal{N}\left(\boldsymbol{f}_{n},10^{-3}\mathbf{I}\right)$.

Sparse updates were constructed using $m=\frac{n}{4}$ inducing variables
$\boldsymbol{u}$, whose distributions $q(\boldsymbol{u})$ and inducing
locations $\mathbf{Z}$ were obtained by minimizing Kullback–Leibler
divergences. Conjugate-gradient-based updates were carried out by, first,
computing partial pivoted Cholesky decompositions in order to precondition
linear systems
$(\mathbf{K}_{n,n}+\sigma^{2}\mathbf{I})\boldsymbol{v}=(\boldsymbol{y}-\boldsymbol{f}_{n}-\boldsymbol{\varepsilon})$.
We then iteratively solved for Gaussian pathwise updates using the method of
conjugate gradients. Stopping conditions for both the partial pivoted Cholesky
decomposition and conjugate gradient solver were chosen to match those of
gardner2018gpytorch (28). Prior to discussing trends in Figure 5, we would
like to point out that curves associated with Gaussian updates (black) are
heavily obscured: in the left column, by CG-based ones (orange and red) and in
top middle and top right plots by RFF-based ones (green).

Comparing the rows of Figure 5, we see that random Fourier feature (RFF)
approximations to priors introduce modest amounts of error in exchange for
large cost reductions. These savings are particularly dramatic in cases where
test inputs $\mathbf{X}_{*}$ significantly outnumber training locations
$\mathbf{X}_{n}$. Echoing discussion in Section 4.5, however, $m$-dimensional
random Fourier bases struggle to represent the data. All other things being
held equal, sparse updates performed in the canonical basis consistently
outperform RFF-based ones. These sparse methods are also considerable faster
than competing approaches when $m\ll n$.

Direct comparison of sparse and CG updates is difficult, since both methods
are sensitive to various design choices. In our experiments, CG-based updates
behaved tantamount to exact ones—with two important caveats. First, CG-based
updates were initially slower than exact ones but outpace them as $n$
increased. Second, naïvely computing pathwise updates using CG is highly
inefficient when it comes to caching. When repeatedly conditioning on
(potentially different realizations of) $\boldsymbol{\gamma}=\boldsymbol{y}$,
one option is to use CG to precompute the matrix inverse
$(\mathbf{K}_{n,n}+\sigma^{2}\mathbf{I})^{-1}$. This CG+ variant is
significantly more cache-friendly, but also much more susceptible to round-off
error—see dashed red curves in Figure 5.

These empirical results help to characterize the behaviors of errors
introduced by different approximation schemes, but leave many questions
unanswered. In order to fill in some of the remaining gaps, we now analyze
various types of approximation error in details.

##  6 Error analysis

Over the course of this section, we will analyze the different types of error
introduced by pathwise approximations. Speaking about these errors requires us
to agree upon a suitable notion of similarity between Gaussian processes.
Ultimately, we are interested in understanding how these approximations
influence Monte Carlo estimators. We therefore focus on $2$-Wasserstein
distances between true and approximate posteriors, since they control
downstream Monte Carlo errors.111111$2$-Wasserstein distances majorize
$1$-Wasserstein distances, which regulate expectations of Lipschitz
functionals by Kantorovich–Rubinstein duality villani08 (91). These distances
measure the similarity of Gaussian processes $\tilde{f}$ and $f$ as the
expectation of a metric $d\big{(}\tilde{f},f\big{)}$ under the best possible
_coupling_ of the two processes. Formally, we have

$W_{2,d}\big{(}\tilde{f},f\big{)}=\left[\inf_{\pi\in\Pi(\tilde{\mu},\mu)}\operatorname{\mathbb{E}}_{\pi}d\big{(}\tilde{f},f\big{)}^{2}\right]^{\nicefrac{{1}}{{2}}},$
(36)

where $\Pi(\tilde{\mu},\mu)$ denotes the set of valid couplings
mallasto2017learning (52), i.e. joint measures whose marginals correspond with
the Gaussian measures $\tilde{\mu}$ and $\mu$ induced by processes $\tilde{f}$
and $f$, respectively. Below, we employ $L^{2}$ and supremum norms as the
underlying metrics used to define $2$-Wasserstein distances.

For the remainder of this section, we assume that the domain $\mathcal{X}$ is
a compact subset of some metric measure space $\mathcal{M}$ and that
$\mathcal{X}$ has finite measure. As a straightforward example, the domain may
be a $d$-dimensional hypercube $\mathcal{X}=[a,b]^{d}$ within
$\mathcal{M}=\mathbb{R}^{d}$.

Lastly, let us introduce some additional notation to simplify material
presented below. First, we will use $\tilde{f}\mid\boldsymbol{y}$ and
$\tilde{f}\mid\boldsymbol{u}$ to denote pathwise conditioning of an
approximate prior $\tilde{f}$ via canonical (16) and sparse (30) update rules,
respectively. These constructions should not be confused with the approximate
posteriors discussed in Sections 4.5 and 5.5. Second, we will superscript
covariance functions $k$ to convey their corresponding processes. For example,
$k^{(\tilde{f})}$ will denote the kernel of the approximation prior
$\tilde{f}$. Third and finally, given a set of $n$ training locations
$\mathbf{X}_{n}\subset\mathcal{X}$, define the _weight function_
$\boldsymbol{\xi}:\mathcal{X}\to\mathbb{R}^{n}$ as

$\boldsymbol{\xi}(\cdot)=k(\mathbf{X}_{n},\mathbf{X}_{n})^{-1}k(\mathbf{X}_{n},\cdot).$
(37)

Variants of this function have been extensively studied in the context of
regression; see silverman1984spline (75, 81) and references contained therein.

###  6.1 Posterior approximation errors

This section adapts the results of wilson20 (97) to study the error in the
decoupled approximate posterior

$(\tilde{f}\mid\boldsymbol{y})(\cdot)\overset{\mathop{}\\!\mathrm{d}}{=}\tilde{f}(\cdot)+k(\cdot,\mathbf{X}_{n})\mathbf{K}_{n,n}^{-1}(\boldsymbol{y}-\boldsymbol{\tilde{f}})=\tilde{f}(\cdot)+\boldsymbol{\xi}(\cdot)^{\top}(\boldsymbol{y}-\boldsymbol{\tilde{f}})$
(38)

formed by updating an $\ell$-dimensional approximate priors
$\tilde{f}(\cdot)=\boldsymbol{\phi}(\cdot)^{\top}\boldsymbol{w}$ via an
$n$-dimensional canonical basis $k(\cdot,\mathbf{X}_{n})$ so as to satisfy the
condition imposed by $n$ noise-free observations $\boldsymbol{y}$.

###### Proposition 5

Assume that $\mathcal{X}\subset\mathbb{R}^{d}$ is compact and that the
stationary kernel $k$ is sufficiently regular for $f\sim\mathcal{GP}(\mu,k)$
to be almost surely continuous. Accordingly, if we define
$C_{1}=\sqrt{2}\operatorname{diam}(\mathcal{X})^{\nicefrac{{d}}{{2}}}\big{(}1+\lVert
k\rVert_{\mathcal{C}(\mathcal{X}^{2})}^{2}\lVert\mathbf{K}_{n,n}^{-1}\rVert_{L(\ell^{\infty};\ell^{1})}^{2}\big{)}^{\nicefrac{{1}}{{2}}},$
then we have

$W_{2,L^{2}(\mathcal{X})}\big{(}\tilde{f}\mid\boldsymbol{y},f\mid\boldsymbol{y}\big{)}=\left(\inf_{\pi\in\Pi(\tilde{\mu},\mu)}\operatorname{\mathbb{E}}_{\pi}\left\lVert(\tilde{f}\mid\boldsymbol{y})-(f\mid\boldsymbol{y})\right\rVert_{L^{2}(\mathcal{X})}^{2}\right)^{\nicefrac{{1}}{{2}}}\leq
C_{1}W_{2,\mathcal{C}(\mathcal{X})}\big{(}\tilde{f},f\big{)},$ (39)

where $W_{2,L^{2}(\mathcal{X})}$ and $W_{2,\mathcal{C}(\mathcal{X})}$
respectively denote 2-Wasserstein distances over the Lebesgue space
$L^{2}(\mathcal{X})$ and the space of continuous functions
$\mathcal{C}(\mathcal{X})$ equipped with the supremum norm,
$\left\lVert\cdot\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}$ is the supremum
norm over continuous functions, and
$\lVert\cdot\rVert_{L(\ell^{\infty};\ell^{1})}$ is the operator norm between
$\ell^{\infty}$ and $\ell^{1}$ spaces.

Proof  We begin by considering the term inside the expectation in (39).
Applying Matheron’s rule followed by Hölder’s inequality ($p=1$, $q=\infty$),
we have

$\displaystyle\begin{split}\left|(\tilde{f}\mid\boldsymbol{y})(\boldsymbol{x})-(f\mid\boldsymbol{y})(\boldsymbol{x})\right|^{2}&\leq
2\left|\tilde{f}(\boldsymbol{x})-f(\boldsymbol{x})\right|^{2}+2\left|\boldsymbol{\xi}(\boldsymbol{x})^{\top}(\boldsymbol{\tilde{f}}_{n}-\boldsymbol{f}_{n})\right|^{2}\\\
&\leq
2\left\lVert\tilde{f}-f\right\rVert_{L^{\infty}(\mathcal{X})}^{2}+2\left\lVert\boldsymbol{\xi}(\boldsymbol{x})\right\rVert_{\ell^{1}}^{2}\left\lVert\boldsymbol{\tilde{f}}_{n}-\boldsymbol{f}_{n}\right\rVert_{\ell^{\infty}}^{2}.\end{split}$
(40)

Continuing from the second line, the definition of the operator norm implies
that

$\displaystyle\begin{split}\left|(\tilde{f}\mid\boldsymbol{y})(\boldsymbol{x})-(f\mid\boldsymbol{y})(\boldsymbol{x})\right|^{2}&\leq
2\left(1+\left\lVert
k(\boldsymbol{x},\mathbf{X}_{n})\right\rVert_{\ell^{\infty}}^{2}\left\lVert\mathbf{K}_{n,n}^{-1}\right\rVert_{L(\ell^{\infty};\ell^{1})}^{2}\right)\left\lVert\tilde{f}-f\right\rVert_{L^{\infty}(\mathcal{X})}^{2}\\\
&\leq
2\left(1+\big{\lVert}k\big{\rVert}_{\mathcal{C}(\mathcal{X}^{2})}^{2}\left\lVert\mathbf{K}_{n,n}^{-1}\right\rVert_{L(\ell^{\infty};\ell^{1})}^{2}\right)\left\lVert\tilde{f}-f\right\rVert_{L^{\infty}(\mathcal{X})}^{2}\\\
&=\underbracket{2\left(1+\big{\lVert}k\big{\rVert}_{\mathcal{C}(\mathcal{X}^{2})}^{2}\left\lVert\mathbf{K}_{n,n}^{-1}\right\rVert_{L(\ell^{\infty};\ell^{1})}^{2}\right)}_{C_{0}}\left\lVert\tilde{f}-f\right\rVert_{\mathcal{C}(\mathcal{X})}^{2},\end{split}$
(41)

where, in the final line, we have used continuity of sample paths to replace
$\left\lVert\cdot\right\rVert_{L^{\infty}(\mathcal{X})}$ with
$\left\lVert\cdot\right\rVert_{\mathcal{C}(\mathcal{X})}$. We now lift this
bound between sample paths to one on 2-Wasserstein distances by integrating
both sides with respect to the optimal coupling $\pi\in\Pi(\tilde{\mu},\mu)$

$\displaystyle\begin{split}W_{2,L^{2}(\mathcal{X})}\big{(}\tilde{f}\mid\boldsymbol{y},f\mid\boldsymbol{y}\big{)}&=\left(\inf_{\pi\in\Pi(\tilde{\mu},\mu)}\operatorname{\mathbb{E}}_{\pi}\left\lVert(\tilde{f}\mid\boldsymbol{y})-(f\mid\boldsymbol{y})\right\rVert_{L^{2}(\mathcal{X})}^{2}\right)^{\nicefrac{{1}}{{2}}}\\\
&\leq\left(C_{0}\operatorname{vol}(\mathcal{X})\inf_{\pi\in\Pi(\tilde{\mu},\mu)}\operatorname{\mathbb{E}}_{\pi}\left\lVert\tilde{f}-f\right\rVert_{\mathcal{C}(\mathcal{X})}^{2}\right)^{\nicefrac{{1}}{{2}}}\\\
&\leq
C_{1}W_{2,\mathcal{C}(\mathcal{X})}\big{(}\tilde{f},f\big{)},\end{split}$ (42)

where $\operatorname{vol}(\mathcal{X})$ denotes the Lebesgue measure of
$\mathcal{X}$. Hence, the claim follows.121212Note that, since $f$ is sample-
continuous and $\mathcal{C}(\mathcal{X})$ is a separable metric space,
$W_{2,\mathcal{C}(\mathcal{X})}$ is a proper metric.

###### Proposition 6

With the same assumptions, let
$C_{2}=n\left(1+\lVert\mathbf{K}^{-1}_{n,n}\rVert_{\mathcal{C}(\mathcal{X}^{2})}\lVert
k\rVert_{\mathcal{C}(\mathcal{X}^{2})}\right)^{2}$. Then,

$\mathbb{E}_{\boldsymbol{\phi}}\big{\lVert}k^{(\tilde{f}\mid\boldsymbol{y})}-k^{(f\mid\boldsymbol{y})}\big{\rVert}_{\mathcal{C}(\mathcal{X}^{2})}\leq
C_{2}\mathbb{E}_{\boldsymbol{\phi}}\big{\lVert}k^{(\tilde{f})}-k\big{\rVert}_{\mathcal{C}(\mathcal{X}^{2})}.$
(43)

Moreover, when $\tilde{f}$ is a random Fourier feature approximation of the
prior, it follows that

$\mathbb{E}_{\boldsymbol{\phi}}\big{\lVert}k^{(\tilde{f}\mid\boldsymbol{y})}-k^{(f\mid\boldsymbol{y})}\big{\rVert}_{\mathcal{C}(\mathcal{X}^{2})}\leq\ell^{\nicefrac{{-1}}{{2}}}C_{2}C_{3},$
(44)

where $C_{3}$ is one of several possible constants given by sutherland15 (83).

Proof  Let
$M_{k}:\mathcal{C}(\mathcal{X}\times\mathcal{X})\to\mathcal{C}(\mathcal{X}\times\mathcal{X})$
be the bounded linear operator given by

$\displaystyle(M_{k}c)(\boldsymbol{x},\boldsymbol{x}^{\prime})$
$\displaystyle=c(\boldsymbol{x},\boldsymbol{x}^{\prime})-c(\boldsymbol{x},\mathbf{X}_{n})\boldsymbol{\xi}(\boldsymbol{x}^{\prime})-\boldsymbol{\xi}(\boldsymbol{x})^{\top}c(\mathbf{X}_{n},\boldsymbol{x}^{\prime})+\boldsymbol{\xi}(\boldsymbol{x})^{\top}c(\mathbf{X}_{n},\mathbf{X}_{n})\boldsymbol{\xi}(\boldsymbol{x}^{\prime}).$
(45)

Henceforth, we omit the subscript from $M_{k}$ to ease notation. Note that, by
construction,

$\displaystyle
k^{(f\mid\boldsymbol{y})}(\boldsymbol{x},\boldsymbol{x}^{\prime})$
$\displaystyle=(Mk)(\boldsymbol{x},\boldsymbol{x}^{\prime})$ $\displaystyle
k^{(\tilde{f}\mid\boldsymbol{y})}(\boldsymbol{x},\boldsymbol{x}^{\prime})$
$\displaystyle=(Mk^{(\tilde{f})})(\boldsymbol{x},\boldsymbol{x}^{\prime}).$
(46)

Focusing on the integrand on the left-hand side of (43), we begin by
separating out the operator norm $\left\lVert
M\right\rVert_{L(\mathcal{C}(\mathcal{X}^{2});\mathcal{C}(\mathcal{X}^{2}))}$
as

$\left\lVert
k^{(\tilde{f}\mid\boldsymbol{y})}-k^{(f\mid\boldsymbol{y})}\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}=\left\lVert
Mk^{(\tilde{f})}-Mk\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}\leq\left\lVert
M\right\rVert_{L(\mathcal{C}(\mathcal{X}^{2});\mathcal{C}(\mathcal{X}^{2}))}\left\lVert
k^{(\tilde{f})}-k\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}.$ (47)

Refining this inequality requires us to upper bound $\left\lVert
M\right\rVert_{L(\mathcal{C}(\mathcal{X}^{2});\mathcal{C}(\mathcal{X}^{2}))}$.
To do so, we write

$\displaystyle\left\lVert Mc\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}$
$\displaystyle\leq\left\lVert
c\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}+2\left\lVert
c(\cdot,\mathbf{X}_{n})\boldsymbol{\xi}(\cdot)\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}+\left\lVert\boldsymbol{\xi}(\cdot)^{\top}c(\mathbf{X}_{n},\mathbf{X}_{n})\boldsymbol{\xi}(\cdot)\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}.$
(48)

We now use Hölder’s inequality ($p=1$, $q=\infty$) followed by the definition
of the operator norm
$\left\lVert\cdot\right\rVert_{L(\ell^{\infty};\ell^{1})}$ to bound the second
and third terms on the right as

$\displaystyle\begin{split}\left\lVert
c(\cdot,\mathbf{X}_{n})\boldsymbol{\xi}(\cdot)\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}&=\sup_{\boldsymbol{x},\boldsymbol{x}^{\prime}\in\mathcal{X}}\left[c(\boldsymbol{x},\mathbf{X}_{n})\boldsymbol{\xi}(\boldsymbol{x}^{\prime})\right]\\\
&\leq\sup_{\boldsymbol{x},\boldsymbol{x}^{\prime}\in\mathcal{X}}\left[\left\lVert
c(\boldsymbol{x},\mathbf{X}_{n})\right\rVert_{\ell^{\infty}}\left\lVert\mathbf{K}_{n,n}^{-1}\right\rVert_{L(\ell^{\infty};\ell^{1})}\left\lVert
k(\mathbf{X}_{n},\boldsymbol{x}^{\prime})\right\rVert_{\ell^{\infty}}\right]\\\
&\leq\left\lVert
c\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}\left\lVert\mathbf{K}^{-1}_{n,n}\right\rVert_{L(\ell^{\infty};\ell^{1})}\left\lVert
k\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}\end{split}$ (49)

and

$\left\lVert\boldsymbol{\xi}(\cdot)^{\top}c(\mathbf{X}_{n},\mathbf{X}_{n})\boldsymbol{\xi}(\cdot)\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}\leq
n\left\lVert
c\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}\left\lVert\mathbf{K}^{-1}_{n,n}\right\rVert_{L(\ell^{\infty};\ell^{1})}^{2}\left\lVert
k\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}^{2}.$ (50)

Returning to (48), we may now bound $\left\lVert
Mc\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}$ by writing

$\displaystyle\begin{split}\left\lVert
Mc\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}&\leq\left\lVert
c\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}\Big{(}1+2\left\lVert\mathbf{K}^{-1}_{n,n}\right\rVert_{L(\ell^{\infty};\ell^{1})}\left\lVert
k\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}+n\left\lVert\mathbf{K}^{-1}_{n,n}\right\rVert_{L(\ell^{\infty};\ell^{1})}^{2}\left\lVert
k\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}^{2}\Big{)}\\\ &\leq\left\lVert
c\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}\Big{(}n\left[1+\left\lVert\mathbf{K}^{-1}_{n,n}\right\rVert_{L(\ell^{\infty};\ell^{1})}\left\lVert
k\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}\right]^{2}\Big{)},\end{split}$
(51)

which immediately implies that

$\left\lVert
M\right\rVert_{L(\mathcal{C}(\mathcal{X}^{2});\mathcal{C}(\mathcal{X}^{2}))}=\sup_{c\neq
0}\frac{\left\lVert Mc\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}}{\left\lVert
c\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}}\leq
n\left[1+\left\lVert\mathbf{K}^{-1}_{n,n}\right\rVert_{L(\ell^{\infty};\ell^{1})}\left\lVert
k\right\rVert_{\mathcal{C}(\mathcal{X}^{2})}\right]^{2}.$ (52)

Note that, since this bound is independent of the particular realization of
the $\ell$-dimensional random Fourier basis $\boldsymbol{\phi}$ used to
construct the approximate prior $\tilde{f}$, it is constant with respect to
the expectation (43). Finally, sutherland15 (83) have shown that there exists
a constant $C_{3}$ such that

$\operatorname{\mathbb{E}}_{\boldsymbol{\phi}}\big{\lVert}k^{(\tilde{f})}-k\big{\rVert}_{\mathcal{C}(\mathcal{X}^{2})}\leq\ell^{\nicefrac{{-1}}{{2}}}C_{3}.$
(53)

Combining this inequality with the preceding ones gives the result.

Together, Propositions 5 and 6 show that error in the approximate prior
$\tilde{f}$ controls the error in the resulting approximate posterior
$\tilde{f}\mid\boldsymbol{y}$. These bounds are not tight, seeing as constants
$C_{1}$ and $C_{2}$ both depend on and may grow with the number of
observations $n$. Based on this observation, it is tempting to think that the
error in $\tilde{f}\mid\boldsymbol{y}$ therefore increases in $n$.
Empirically, however, the opposite trend is observed: the error in
$\tilde{f}\mid\boldsymbol{y}$ actually diminishes as $n$ grows wilson20 (97).
To better understand this behavior, we now study the conditions under which a
pathwise update may counteract the error introduced by an approximate prior.

###  6.2 Contraction of approximate posteriors with noise-free observations

This section formalizes the following syllogism: (i) the true posterior
$f\mid\boldsymbol{y}$ and the approximate posterior
$\tilde{f}\mid\boldsymbol{y}$ have the same mean; (ii) as $n$ increases, both
posteriors contract to their respective means; (iii) therefore, as $n$
increases, the error introduced by the approximate prior $\tilde{f}$ washes
out.

To begin, let $\boldsymbol{\phi}:\mathcal{M}\to\mathbb{R}^{\ell}$ be an
$\ell$-dimensional feature map on an ambient space $\mathcal{M}$ consisting of
linearly independent basis functions $\phi_{i}$. We will say that $\tilde{f}$
is a _standard normal Bayesian linear model_ if it admits the representation

$\displaystyle\tilde{f}(\cdot)$
$\displaystyle=\sum_{i=1}^{\ell}w_{i}\phi_{i}(\cdot)$ $\displaystyle w_{i}$
$\displaystyle\sim\mathcal{N}(0,1).$ (54)

This description includes the Karhunen–Loève and Fourier feature
approximations described in Section 4. As before, let
$\mathbf{\Phi}_{n}=\boldsymbol{\phi}(\mathbf{X}_{n})$ be an $n\times\ell$
feature matrix and $\mathcal{H}_{k}$ be the reproducing kernel Hilbert space
associated with a kernel $k$. We say that a function $\phi_{i}$ _lies locally_
in $\mathcal{H}_{k}$ for a compact $\mathcal{X}\subseteq\mathcal{M}$ if there
exists a function $\psi_{j}\in\mathcal{H}_{k}$ that agrees with $\phi_{i}$ on
$\mathcal{X}$, i.e.
${\left.\kern-1.2pt\phi_{i}\vphantom{\big{|}}\right|_{\mathcal{X}}}={\left.\kern-1.2pt\psi_{j}\vphantom{\big{|}}\right|_{\mathcal{X}}}$.

When $\mathcal{M}$ is a compact metric space, the eigenfunctions $\phi_{i}$
used to construct (truncated) Karhunen–Loève expansions belong to
$\mathcal{H}_{k}$ by construction. More generally, assessing whether or not
$\phi_{i}$ lies locally in $\mathcal{H}_{k}$ is often straightforward for
kernels with known reproducing kernel Hilbert spaces. As a concrete example,
the RKHS of a Matérn-$\nu$ kernel is the Sobolev space of order
$\kappa=\nu+\nicefrac{{d}}{{2}}$. For integer values of $\kappa$, this is the
space of square-integrable functions with $\kappa$ square-integrable weak
derivatives. Trigonometric basis functions
$\phi_{i}(\boldsymbol{x})=\cos(2\pi\boldsymbol{\omega}_{i}^{\top}\boldsymbol{x}+\tau_{i})$
can readily be adapted to satisfy this requirement. Specifically, we may
multiply them by a suitably chosen, infinitely-differentiable function that
ensures they decay to zero outside of $\mathcal{X}$, such that the resulting
basis functions (and their derivatives) are square-integrable.

We are now ready to state and prove the primary claim. In the following,
Proposition 7 and Corollary 8 will demonstrate that
$\tilde{f}\mid\boldsymbol{y}$ contracts at the same rate as
$f\mid\boldsymbol{y}$. Subsequently, Corollary 9 will show that the error in
$\tilde{f}\mid\boldsymbol{y}$ vanishes as $n\to\infty$ in any reasonable limit
where the variance of the true posterior contracts to zero everywhere on
$\mathcal{X}$.

###### Proposition 7

Suppose $\mathcal{X}\subseteq\mathcal{M}$ is compact and that each of the
$\ell$ basis functions $\phi_{i}$ used to construct the standard normal
Bayesian linear model $\tilde{f}$ lies locally in $\mathcal{H}_{k}$. If the
points $\mathbf{X}_{n}\subset\mathcal{X}$ used to condition the approximate
posterior $\tilde{f}\mid\boldsymbol{y}$ are chosen such that
$f\mid\boldsymbol{y}$ satisfies
$\sup_{\boldsymbol{x}\in\mathcal{X}}k^{(f\mid\boldsymbol{y})}(\boldsymbol{x},\boldsymbol{x})\leq\varepsilon$,
then it follows that131313This result holds even when the weights are not
assumed i.i.d., albeit with a slightly different constant.

$\sup_{\boldsymbol{x}\in\mathcal{X}}\big{|}k^{(\tilde{f}\mid\boldsymbol{y})}(\boldsymbol{x},\boldsymbol{x})\big{|}\leq
C_{4}\varepsilon,$ (55)

where we have defined
$C_{4}=\ell\max_{i}\inf\left\\{\lVert\psi_{i}\rVert_{\mathcal{H}_{k}}^{2}:\psi_{i}|_{\mathcal{X}}=\phi_{i}|_{\mathcal{X}},\forall\psi_{i}\in\mathcal{H}_{k}\right\\}$.

Proof  Recall from (38) we can use the weight function
$\boldsymbol{\xi}(\cdot)=k(\mathbf{X}_{n},\mathbf{X}_{n})^{-1}k(\mathbf{X}_{n},\cdot)$
to express the approximate posterior as
$(\tilde{f}\mid\boldsymbol{y})(\cdot)\,\smash{\overset{\mathop{}\\!\mathrm{d}}{=}}\,\boldsymbol{\phi}(\cdot)^{\top}\boldsymbol{w}-\boldsymbol{\xi}(\cdot)^{\top}(\boldsymbol{y}-\mathbf{\Phi}_{n}\boldsymbol{w})$.
Under this notation, it is clear that we may immediately upper bound the
variance of the $\tilde{f}\mid\boldsymbol{y}$ as

$\begin{split}\operatorname{Var}\left((\tilde{f}\mid\boldsymbol{y})(\cdot)\right)=\operatorname{\mathbb{E}}\left[\left(\boldsymbol{\phi}(\cdot)^{\top}-\boldsymbol{\xi}(\cdot)^{\top}\mathbf{\Phi}_{n}\right)\boldsymbol{w}\right]^{2}\leq\ell\max_{i}\left(\phi_{i}(\cdot)-\boldsymbol{\xi}(\cdot)^{\top}\phi_{i}(\mathbf{X}_{n})\right)^{2},\end{split}$
(56)

where, on the right, we have used the fact that
$\operatorname{\mathbb{E}}\left\lVert\boldsymbol{w}\right\rVert^{2}=\ell$. By
further denoting $\mathcal{G}=\\{g\in\mathcal{H}_{k}:\left\lVert
g\right\rVert_{\mathcal{H}_{k}}=1\\}$, we may now exploit the dual
representation of the RKHS norm to write

$\displaystyle\begin{split}\left|\phi_{i}(\boldsymbol{x}_{*})-\boldsymbol{\xi}(\boldsymbol{x}_{*})^{\top}\phi_{i}(\mathbf{X}_{n})\right|&\leq\left\lVert\phi_{i}\right\rVert_{\mathcal{H}_{k}}\sup_{g\in\mathcal{G}}\left|g(\boldsymbol{x}_{*})-\boldsymbol{\xi}(\boldsymbol{x}_{*})^{\top}g(\mathbf{X}_{n})\right|\\\
&=\left\lVert\phi_{i}\right\rVert_{\mathcal{H}_{k}}\left\lVert
k(\cdot,\boldsymbol{x}_{*})-\boldsymbol{\xi}(\boldsymbol{x}_{*})^{\top}\mathbf{K}_{n,*}\right\rVert_{\mathcal{H}_{k}}\\\
&=\left\lVert\phi_{i}\right\rVert_{\mathcal{H}_{k}}\underbracket{\sqrt{k(\boldsymbol{x}_{*},\boldsymbol{x}_{*})-\mathbf{K}_{*,n}\mathbf{K}_{n,n}^{-1}\mathbf{K}_{n,*}}}_{\mathcal{P}_{\mathbf{X}}(\boldsymbol{x}_{*})},\end{split}$
(57)

where, because $\phi_{i}$ lies locally in $\mathcal{H}_{k}$, we may replace it
with any $\psi_{i}\in\mathcal{H}_{k}:\
\psi_{i}|_{\mathcal{X}}=\phi_{i}|_{\mathcal{X}}$. Noting that
$\mathcal{P}_{\mathbf{X}}(\cdot)=\sqrt{k^{(f\mid\boldsymbol{y})}(\cdot,\cdot)}$
and collecting terms gives the result.

###### Corollary 8

With the same assumptions, as
$\sup_{\boldsymbol{x}\in\mathcal{X}}k^{(f\mid\boldsymbol{y})}(\boldsymbol{x},\boldsymbol{x})\rightarrow
0$, it follows that

$\sup_{\boldsymbol{x},\boldsymbol{x}^{\prime}\in\mathcal{X}}\big{|}k^{(\tilde{f}\mid\boldsymbol{y})}(\boldsymbol{x},\boldsymbol{x}^{\prime})-k^{(f\mid\boldsymbol{y})}(\boldsymbol{x},\boldsymbol{x}^{\prime})\big{|}\rightarrow
0.$ (58)

Proof  Begin by applying the triangle inequality to the above and,
subsequently, use the Cauchy-Schwartz inequality to bound
$k(\boldsymbol{x},\boldsymbol{x}^{\prime})\leq\sqrt{k(\boldsymbol{x},\boldsymbol{x})}\sqrt{k(\boldsymbol{x}^{\prime},\boldsymbol{x}^{\prime})}$,
which gives

$\displaystyle\begin{split}\sup_{\boldsymbol{x},\boldsymbol{x}^{\prime}\in\mathcal{X}}\big{|}k^{(\tilde{f}\mid\boldsymbol{y})}(\boldsymbol{x},\boldsymbol{x}^{\prime})-k^{(f\mid\boldsymbol{y})}(\boldsymbol{x},\boldsymbol{x}^{\prime})\big{|}&\leq\sup_{\boldsymbol{x},\boldsymbol{x}^{\prime}\in\mathcal{X}}\big{|}k^{(\tilde{f}\mid\boldsymbol{y})}(\boldsymbol{x},\boldsymbol{x}^{\prime})\big{|}+\sup_{\boldsymbol{x},\boldsymbol{x}^{\prime}\in\mathcal{X}}\big{|}k^{(f\mid\boldsymbol{y})}(\boldsymbol{x},\boldsymbol{x}^{\prime})\big{|}\\\
&\leq\sup_{\boldsymbol{x}\in\mathcal{X}}\big{|}k^{(\tilde{f}\mid\boldsymbol{y})}(\boldsymbol{x},\boldsymbol{x})\big{|}+\sup_{\boldsymbol{x}\in\mathcal{X}}\big{|}k^{(f\mid\boldsymbol{y})}(\boldsymbol{x},\boldsymbol{x})\big{|}.\end{split}$
(59)

In the final expression, convergence of the former term is given by
Proposition 7, while the latter goes to zero by assumption.

###### Corollary 9

With the same assumptions, as
$\sup_{\boldsymbol{x}\in\mathcal{X}}k^{(f\mid\boldsymbol{y})}(\boldsymbol{x},\boldsymbol{x})\rightarrow
0$, it follows that

$W_{2,L^{2}(\mathcal{X})}(f\mid\boldsymbol{y},\tilde{f}\mid\boldsymbol{y})\rightarrow
0.$ (60)

Proof  Since $L^{2}(\mathcal{X})$ is a normed space and
$\operatorname{\mathbb{E}}(f\mid\boldsymbol{y})=\operatorname{\mathbb{E}}(\tilde{f}\mid\boldsymbol{y})$,
we have that

$W_{2,L^{2}(\mathcal{X})}(f\mid\boldsymbol{y},\tilde{f}\mid\boldsymbol{y})=W_{2,L^{2}(\mathcal{X})}\bigg{(}\underbracket{f(\cdot)-\boldsymbol{\xi}(\cdot)^{\top}f(\mathbf{X}_{n})}_{(f\mid\boldsymbol{y})_{0}},\,\underbracket{\boldsymbol{\phi}(\cdot)^{\top}\boldsymbol{w}-\boldsymbol{\xi}(\cdot)^{\top}\mathbf{\Phi}_{n}\boldsymbol{w}}_{(\tilde{f}\mid\boldsymbol{y})_{0}}\bigg{)},$
(61)

where $(f\mid\boldsymbol{y})_{0}$ and $(\tilde{f}\mid\boldsymbol{y})_{0}$
denote centered processes. Now, let $\mathbbold{0}$ be an almost surely zero
stochastic process over $\mathcal{X}$. Then, by the triangle inequality,

$W_{2,L^{2}(\mathcal{X})}\Big{(}(f\mid\boldsymbol{y})_{0},(\tilde{f}\mid\boldsymbol{y})_{0}\Big{)}\leq
W_{2,L^{2}(\mathcal{X})}\Big{(}(f\mid\boldsymbol{y})_{0},\mathbbold{0}\Big{)}+W_{2,L^{2}(\mathcal{X})}\Big{(}(\tilde{f}\mid\boldsymbol{y})_{0},\mathbbold{0}\Big{)}.$
(62)

Expanding the definition of Wasserstein distances $W_{2,L^{2}(\mathcal{X})}$
before using Tonelli’s theorem to change the order of integration gives

$\displaystyle
W_{2,L^{2}(\mathcal{X})}\Big{(}(f\mid\boldsymbol{y})_{0},(\tilde{f}\mid\boldsymbol{y})_{0}\Big{)}$
$\displaystyle\leq\Big{(}\operatorname{\mathbb{E}}\big{\lVert}(f\mid\boldsymbol{y})_{0}-\mathbbold{0}\big{\rVert}_{L^{2}(\mathcal{X})}^{2}\Big{)}^{\nicefrac{{1}}{{2}}}+\Big{(}\operatorname{\mathbb{E}}\big{\lVert}(\tilde{f}\mid\boldsymbol{y})_{0}-\mathbbold{0}\big{\rVert}_{L^{2}(\mathcal{X})}^{2}\Big{)}^{\nicefrac{{1}}{{2}}}$
$\displaystyle=\Big{(}\int_{\mathcal{X}}k^{(f\mid\boldsymbol{y})}(\boldsymbol{x},\boldsymbol{x})\mathop{}\\!\mathrm{d}\boldsymbol{x}\Big{)}^{\nicefrac{{1}}{{2}}}+\Big{(}\int_{\mathcal{X}}k^{(\tilde{f}\mid\boldsymbol{y})}(\boldsymbol{x},\boldsymbol{x})\mathop{}\\!\mathrm{d}\boldsymbol{x}\Big{)}^{\nicefrac{{1}}{{2}}},$
(63)

where both terms in the final expression converge to zero by compactness of
$\mathcal{X}$ together with Proposition 7.
Together, these claims demonstrate that the decoupled approximate posterior
$\tilde{f}\mid\boldsymbol{y}$, formed by using the canonical basis
$k(\cdot,\mathbf{X}_{n})$ to update a well-specified approximate prior
$\tilde{f}$, _inherits_ the contractive properties of the true posterior
$f\mid\boldsymbol{y}$.

Per the beginning of this section, approximate priors $\tilde{f}$ defined as
standard normal Bayesian linear models with basis functions that lie locally
in $\mathcal{H}_{k}$ are well-specified. The following counterexample helps
clarify what can happen when $\tilde{f}$ is misspecified. Consider an
approximate prior $\tilde{f}\sim\mathcal{GP}(0,\delta)$ equipped with the
Kronecker delta kernel $\delta$ such that
$\operatorname{Cov}\big{(}\tilde{f}(\boldsymbol{x}_{i}),\tilde{f}(\boldsymbol{x}_{j})\big{)}=1$
if $\boldsymbol{x}_{i}=\boldsymbol{x}_{j}$ and $0$ otherwise. Given a finite
set of test locations
$\mathbf{X}_{*}\subset\mathcal{X}\setminus\mathbf{X}_{n}$, let
$\mathbf{\Xi}=\boldsymbol{\xi}(\mathbf{X}_{*})^{\top}$. Applying the pathwise
update (17) to $\tilde{f}$, the posterior covariance is then

$\operatorname{Cov}\big{(}\boldsymbol{\tilde{f}}_{*}\mid\boldsymbol{y}\big{)}=\operatorname{Cov}\big{(}\boldsymbol{\tilde{f}}_{*}\big{)}+\mathbf{\Xi}\operatorname{Cov}\big{(}\boldsymbol{\tilde{f}}_{n}\big{)}\mathbf{\Xi}^{\top}-2\operatorname{Cov}\big{(}\boldsymbol{\tilde{f}}_{*},\boldsymbol{\tilde{f}}_{n}\big{)}\mathbf{\Xi}^{\top}=\mathbf{I}+\mathbf{K}_{*,n}^{\vphantom{-2}}\mathbf{K}_{n,n}^{-2}\mathbf{K}_{n,*}^{\vphantom{-2}}.$
(64)

Since the second of the two terms on the right is guaranteed non-negative, the
variance of the resulting posterior is bounded from below by $1$. For this
choice of $\tilde{f}$, then, the approximation error inherent to
$\tilde{f}\mid\boldsymbol{y}$ does not diminish as $n$
increases.141414Contraction of the true posterior is well-studied and has
strong ties to the literature on kernel methods. kanagawa2018gaussian (40)
reviews these connections in greater detail: there, Theorem 5.4 shows how the
_power function_ $\mathcal{P}_{\mathbf{X}}$ can be bounded in terms of the
fill distance
$h(\mathbf{X}_{n})=\sup_{\boldsymbol{x}_{*}\in\mathcal{X}}\mathop{\mathrm{inf}\vphantom{\mathrm{sup}}}_{\boldsymbol{x}\in\mathbf{X}_{n}}\left\lVert\boldsymbol{x}_{*}-\boldsymbol{x}\right\rVert$.

###  6.3 Sparse approximation errors

We now examine the error introduced by using a sparse pathwise update (30) to
construct an approximate posterior. As notation, we write
$f\mid\boldsymbol{u}$ and $\tilde{f}\mid\boldsymbol{u}$ for the approximate
posteriors formed by applying the sparse update to the true prior $f$ and to
the approximate prior $\tilde{f}$, respectively. Results discussed here mirror
those presented by wilson20 (97). Appealing to the triangle inequality, we
have

$\displaystyle\begin{split}W_{2,L^{2}(\mathcal{X})}\big{(}\tilde{f}\mid\boldsymbol{u},f\mid\boldsymbol{y}\big{)}&\leq\underbracket{W_{2,L^{2}(\mathcal{X})}\big{(}\tilde{f}\mid\boldsymbol{u},f\mid\boldsymbol{u}\big{)}}_{\hskip
39.18051pt\mathclap{\mathrel{\text{missing}}{errorinapproximateprior}}\hskip
39.18051pt}+\underbracket{W_{2,L^{2}(\mathcal{X})}\big{(}f\mid\boldsymbol{u},f\mid\boldsymbol{y}\big{)}}_{\hskip
37.67357pt\mathclap{\mathrel{\text{missing}}{errorinsparseupdate}}\hskip
37.67357pt}\\\\[-2.69096pt]
\mathbb{E}_{\boldsymbol{\phi}}\big{\lVert}k^{(\tilde{f}\mid\boldsymbol{u})}-k^{(f\mid\boldsymbol{y})}\big{\rVert}_{\mathcal{C}(\mathcal{X}^{2})}&\leq\overbracket{\mathbb{E}_{\boldsymbol{\phi}}\big{\lVert}k^{(\tilde{f}\mid\boldsymbol{u})}-k^{(f\mid\boldsymbol{u})}\big{\rVert}_{\mathcal{C}(\mathcal{X}^{2})}}_{\smash{\hskip
78.36102pt}}+\overbracket{\big{\lVert}k^{(f\mid\boldsymbol{u})}-k^{(f\mid\boldsymbol{y})}\big{\rVert}_{\mathcal{C}(\mathcal{X}^{2})}\vphantom{\mathbb{E}_{\boldsymbol{\phi}}\big{\lVert}k^{(\tilde{f}\mid\boldsymbol{u})}-k^{(f\mid\boldsymbol{u})}\big{\rVert}_{\mathcal{C}(\mathcal{X}^{2})}}}_{\smash{\hskip
75.34714pt}}.\end{split}$ (65)

From here, any of the previously presented propositions enable us to control
the total error. For the first terms on the right, the same arguments as
before lead to the same results; however, the constants involved will change,
since the sparse update now assumes the role of the canonical one. The latter
terms do not involve the approximate prior and are therefore beyond the scope
of our present analysis. Note that similar statements hold for the Gaussian
pathwise update (28).

As a final remark, note that we may reduce the total error (65) by
incorporating additional basis functions $k(\cdot,\mathbf{X})$ into the sparse
update. Conceptually, the act of _augmenting_ a sparse update amounts to
replacing $\boldsymbol{u}\sim q(\boldsymbol{u})$ with
$\boldsymbol{u}^{\prime}\sim
q(\boldsymbol{u}^{\prime})=p(\boldsymbol{f}\mid\boldsymbol{u})q(\boldsymbol{u})$,
where $\boldsymbol{f}$ are process values at centers $\mathbf{X}$ rasmussen05
(66, 63). By construction, $q(\boldsymbol{u})$ and
$q(\boldsymbol{u}^{\prime})$ induce the same posterior on $f$. However,
because the augmented update utilizes additional basis functions, the error in
the induced distribution of $\boldsymbol{\tilde{f}}_{*}$ diminishes. This
result follows from the same line of reasoning as before: since
$\operatorname{\mathbb{E}}\big{(}\boldsymbol{f}_{*}\mid\boldsymbol{u}^{\prime}\big{)}=\operatorname{\mathbb{E}}\big{(}\boldsymbol{\tilde{f}}_{*}\mid\boldsymbol{u}^{\prime}\big{)}$,
$f\mid\boldsymbol{u}^{\prime}$ and $\tilde{f}\mid\boldsymbol{u}^{\prime}$
contract to the same function as $|\boldsymbol{u}^{\prime}|\to\infty$. Hence,
the approximate prior washes out and the total error decreases.

##  7 Applications

This section examines the practical consequences of pathwise conditioning in
terms of a curated set of representative tasks. Throughout, we focus on how
pathwise methods for efficiently generating function draws from GP posteriors
enable us to overcome common obstacles and open doors for new research. We
provide a general framework for pathwise conditioning of Gaussian processes
based on GPflow matthews2017gpflow (53).151515Code is available online at
https://github.com/j-wilson/GPflowSampling.

###  7.1 Optimizing black-box functions

Figure 6: Median performances and interquartile ranges of Thompson sampling
methods and popular baselines when optimizing function draws from known GP
priors on $d=\dim(\mathcal{X})$ dimensional domains. Location-scale Thompson
sampling performs well in low-dimensional settings (left), but struggles as
$d$ increase due to its inability to efficiently utilize gradient information.
RFF posteriors enable us to generate function draws, but demand many more
basis functions $b=\ell+n$ than data points $n$ (middle vs. right). Decoupled
approaches using canonical basis functions $k(\cdot,\boldsymbol{x})$ to update
RFF priors $\tilde{f}$ avoids these pitfalls and consistently match or
outperform competing strategies.

Global optimization revolves around the challenge of efficiently identifying a
global minimizer

$\displaystyle\boldsymbol{x}_{\mathrel{\text{missing}}{min}}$
$\displaystyle\in\mathbf{X}_{\mathrel{\text{missing}}{min}}$
$\displaystyle\mathbf{X}_{\mathrel{\text{missing}}{min}}=\operatorname*{\arg\min}_{\boldsymbol{x}\in\mathcal{X}}f(\boldsymbol{x})$
(66)

of a black-box function $f:\mathcal{X}\to\mathbb{R}$. Since $f$ is a black
box, our understanding of its behavior is limited to a set of observations
$\boldsymbol{y}$ at locations $\mathbf{X}_{n}$. Gaussian processes are a
natural and widely used way of representing possible functions
$f\mid\boldsymbol{y}$ movckus1975bayesian (54, 82, 26). In these cases, we
reason about global minimizers (66) in terms of a belief over the random set

$\mathbf{X}_{\mathrel{\text{missing}}{min}}^{(f\mid\boldsymbol{y})}=\operatorname*{\arg\min}_{\boldsymbol{x}\in\mathcal{X}}\,(f\mid\boldsymbol{y})(\boldsymbol{x}).$
(67)

Approaches to these problems are often characterized as striking a balance
between two competing agendas: the need to learn about the function’s global
behavior by _exploring_ the domain $\mathcal{X}$ and the need to obtain
(potentially local) minimizers by _exploiting_ what is already known.

Thompson sampling is a classic decision-making strategy that balances the
tradeoff between exploration and exploitation by sampling actions
$\boldsymbol{x}\in\mathcal{X}$ in proportion to the probability that
$\boldsymbol{x}\in\mathbf{X}_{\mathrel{\text{missing}}{min}}^{(f\mid\boldsymbol{y})}$
thompson33 (84). At first glance, this task may seem daunting, since
$\mathbf{X}_{\mathrel{\text{missing}}{min}}^{(f\mid\boldsymbol{y})}$ is
random. For a given draw of $f\mid\boldsymbol{y}$, however,
$\mathbf{X}_{\mathrel{\text{missing}}{min}}^{(f\mid\boldsymbol{y})}$ is
deterministic. Accordingly, we may Thompson sample an action by generating a
function $f\mid\boldsymbol{y}$ and, subsequently, finding a pathwise global
minimizer.

Thompson sampling’s relative simplicity makes it a natural test bed for
evaluating different sampling strategies, while its real-world performance
chapelle11 (10) assures its ongoing relevance in applied settings. A key
strength of these methods is that they support embarassingly-parallel batch
selection hernandez2017parallel (34, 41). While many GP-based search
strategies allow us to choose $\kappa>1$ queries at a time snoek2012practical
(78, 98), their compute costs tend to scale aggressively in $\kappa$.
Especially when evaluations can be carried out in parallel, then, Thompson
sampling provides an affordable alternative to comparable approaches.

We considered three different variants of Thompson sampling, corresponding
with different approaches to sampling from GP posteriors. The first approach
samples random vectors $\boldsymbol{f}_{*}\mid\boldsymbol{y}$ using location-
scale transforms (19); the second approximates posteriors with Bayesian linear
models; and, the third updates function draws from $\ell$-dimensional
approximate priors $\tilde{f}=\boldsymbol{\phi}(\cdot)^{\top}\boldsymbol{w}$
using canonical basis functions centered at the $n$ training
locations.161616Equation (34) highlights the difference between the second and
third approaches. For fair comparison, we allocate $b=\ell+n$ random Fourier
basis functions to Bayesian linear models employed by the second approach.

At each round of Thompson sampling, we began by sampling process values
$f_{i}\mid\boldsymbol{y}$ independently on a randomly generated discretization
of $\mathcal{X}=[0,1]^{d}$. Next, we constructed a candidate set
$\mathbf{X}_{*}$ using the locations that produces the smallest realizations
of $f_{i}\mid\boldsymbol{y}$. Under a location-scale approach, we then jointly
sampled process values at $|\mathbf{X}_{*}|=2048$ candidates. For both of the
alternatives, we instead used $|\mathbf{X}_{*}|=32$ candidates to initialize
multi-start gradient descent. In all three cases, queries were chosen as
minimizers of the resulting vector $\boldsymbol{f}_{*}\mid\boldsymbol{y}$.
Batches of queries were obtained using $\kappa$ independent runs of this
algorithm.

To eliminate confounding variables, we experimented with black-box functions
drawn from a known Matérn-$\nicefrac{{5}}{{2}}$ prior with an isotropic length
scale $l=\sqrt{\nicefrac{{d}}{{100}}}$ and Gaussian observations
$y\sim\mathcal{N}\left(f(\boldsymbol{x}),10^{-3}\right)$. We set $\kappa=d$,
but this choice was not found to significantly influence our results. Below,
we focus on comparing each Thompson sampling variant’s behavior for different
amounts of design variables $d$ and basis functions $\ell$.

Figure 6 reports key findings based on 32 independent trials; for extended
results, see wilson20 (97). First, location-scale methods’ inability to use
gradient information to efficiently find pathwise minimizers causes its
performance to wane as $d$ increases. In contrast, both of the alternative
variants of Thompson sampling rely on pathwise-differentiable function draws
and, therefore, scale more gracefully in $d$. Second, RFF-based Bayesian
linear models struggle to represent posteriors due to variance starvation
(Section 4.5). As the number of observations $n$ increases relative to the
number of basis functions $b=\ell+n$, the function draws they produce come to
inadequately characterize the true posterior, causing Thompson sampling to
falter. Decoupled approaches to updating $\tilde{f}$ avoid this issue by,
e.g., associating the data with the $n$-dimensional canonical basis
$k(\cdot,\mathbf{X}_{n})$.

###  7.2 Generating boundary-constrained sample paths

Figure 7: Pathwise conditioning of samples from Matérn priors subject to
observations $\boldsymbol{y}$ (black dots) and Dirichlet boundary conditions
$f\big{|}_{\partial\mathcal{X}}=0$. From left to right, the first three
columns show a draw from the prior, a pathwise update, and the corresponding
realization of the posterior. The final two columns communicate the empirical
mean and standard deviation of the posterior, respectively. _Top:_
Illustration of a rectangular domain for which Laplacian eigenpairs are
calculated analytically. _Bottom:_ A non-trivial domain for which the
eigenpairs are approximated numerically.

This section illustrates how techniques introduced in the preceding sections
can be used to efficiently sample Gaussian process posteriors subject to
boundary conditions solin2019know (79). whittle63 (95) showed that a Matérn GP
$f$ defined over $\mathbb{R}^{d}$ satisfies the stochastic partial
differential equation

$\left(\frac{2\nu}{\kappa^{2}}-\Delta\right)^{\frac{\nu}{2}+\frac{d}{4}}f=\mathcal{W},$
(68)

where $\mathcal{W}$ is a (rescaled) white noise process, and $\Delta$ is the
Laplacian. Following solin2019know (79) and rue2005gaussian (68), we restrict
(68) onto a (well-behaved) compact domain $\mathcal{X}\subset\mathbb{R}^{d}$
and impose Dirichlet boundary conditions
${\left.\kern-1.2ptf\vphantom{\big{|}}\right|_{\partial\mathcal{X}}}=0$ to
define a boundary-constrained Matérn Gaussian process over $\mathcal{X}$.
solin2019know (79) demonstrate that such a prior admits the Karhunen–Loève
expansion

$\displaystyle f(\cdot)$
$\displaystyle=\sum_{i=1}^{\infty}w_{i}\phi_{i}(\cdot)$ $\displaystyle
w_{i}\sim\mathcal{N}\bigg{(}0,\frac{\sigma^{2}}{C_{\nu}}\left(\frac{2\nu}{\kappa^{2}}+\lambda_{i}\right)^{-\nu-\frac{d}{2}}\bigg{)},$
(69)

where $\phi_{i}$ are eigenfunctions of the _boundary-constrained_ Laplacian.
We truncate this expansion to obtain the $\ell$-dimensional Bayesian linear
model $\tilde{f}$, which we use together with a pathwise update to construct
the posterior.

Figure 7 visualizes function draws from boundary-constrained priors and
posterior for two choices of boundaries on $\mathbb{R}^{2}$, a rectangle and
the symbol for infinity. Note that eigenfunctions for rectangular regions of
Euclidean domains are available analytically, while those of the infinity
symbol are obtained numerically by solving a Helmholtz equation. Examining
this figure, we see that the sample paths respect the Dirichlet boundary
condition
${\left.\kern-1.2ptf\vphantom{\big{|}}\right|_{\partial\mathcal{X}}}=0$.
Karhunen–Loève expansions enable boundary-constrained GPs, an important class
of non-stationary priors, to be used within the pathwise conditioning
framework.

###  7.3 Simulating dynamical systems

Figure 8: Model-based simulations of a stochastic FitzHugh–Nagumo neuron.
_Left:_ Phase portrait of the true drift function subject to a fixed current
$a=0.5$. _Middle:_ Empirical medians and interquartile ranges of simulated
voltage traces driven by a sinusoidal current (dotted black); ground truth
quartiles are shown in dashed gray. Trajectories generated via location-scale
transforms are summarized on the top in orange, while those produced by
decoupled drift functions are portrayed on the bottom in blue. _Top right:_
Comparison of simulation runtimes. _Bottom right:_ Sinkhorn estimates
cuturi2013sinkhorn (14) to 2-Wasserstein distances between model-based and