Department of Mathematics
Imperial College London
a dissertation submitted for the degree of
Doctor of Philosophy

CHAPTER: DECLARATION

This dissertation and all research contained in it are a product of my original work, except where indicated otherwise by explicit statement or reference.
All ideas, quotations, and data originating from the works of others, published or otherwise, are fully acknowledged according to standard practices of the academic discipline.

CHAPTER: COPYRIGHT

The copyright of this thesis rests with the author. Unless otherwise indicated, its contents are licensed under a Creative Commons Attribution 4.0 International License (CC BY).

Under this license, you may copy and redistribute the material in any medium or format for both commercial and non-commercial purposes.
You may also create and distribute modified versions of the work.
This is on the condition that you credit the author.

When reusing or sharing this work, ensure you make the license terms clear to others by naming the license and linking to the license text.
Where a work has been adapted, you should indicate that the work has been changed and describe those changes.

Please seek permission from the copyright holder for uses of this work that are not included in this license or permitted under UK Copyright Law.

CHAPTER: ACKNOWLEDGMENTS

This thesis is dedicated to the several hundred Twitter users who thought what I was working on was interesting enough to give it a read, or a comment.
Thanks to you, my experiment in writing a thesis in an open-source manner visible to the public was a resounding success.
Getting the opportunity to write this thesis has been a dream come true, and I am delighted to share the experience with you.

The original acknowledgments I had planned to write during most of my doctoral work would have started on a far less positive note.
The sheer volume of kind feedback I've received in the weeks during which I was writing this thesis, almost all of it from people I've never met, has convinced me to leave the past behind, and not write about certain people, actions, and events that do not deserve to be remembered.
Instead, I choose to thank those due to whom I had the chance.

Firstly, I am profoundly grateful to Marc Deisenroth for taking me as a student, and giving me a third chance to complete a Ph.D., and to Seth Flaxman for taking me as a co-supervisee in order to make this possible.
Most people who do not succeed the first time, whatever the reason may be, don't even get a second chance, much less a third one.
You believed in me in the darkest moments of my academic career, and anything that I ever accomplish in my career will only have been possible because of that belief.

Secondly, I am also profoundly grateful to David Draper for believing in me, and taking me as a student when nobody was willing to talk to me or take my ideas seriously.
It is thanks to you that I have a master's degree, and though I was not able to finish a doctorate at the time under your guidance, I hope this thesis shows that you did the right thing in supporting me, even when it became clear I had no academic future at the place I was at, and later in helping me return to scientific work.

Thirdly, I am grateful to Eric Xing  for hiring me when my career was otherwise in the midst of full-blown chaos.
The time I spent at Petuum convinced me that the machine learning community was the right intellectual home for me, and prompted me to look beyond the narrow scope of ideas I had previously been around.
Though I was not truly ready to be an independent researcher at that time, and therefore accomplished little, I am grateful for the opportunity and all that I learned from it.

I thank my thesis committee members, Michael Bronstein and Stefano Ermon, for examining me and for their comments during the viva.
Your careful reading helped significantly improve the quality of this thesis.
I am particularly grateful for Michael's comment that "[t]he introductory chapters, with enough effort, could probably be transformed into a book."
These are among the kindest words my work has ever received.

I also thank my large and ever-growing set of scientific collaborators.
You all give me confidence that, together, we can solve all of the hardest and most important problems!
I am constantly learning from each and every one of you, and none of the work in this thesis would have been possible without your support, ideas, and opportunity to discuss.

I am particularly grateful to Viacheslav Borovitskiy for both the opportunity to collaborate on all of the work presented in this thesis, and to reconnect with my heritage by constantly practicing my Russian.
Together, we were able to write substantially stronger work than any of us, or at least certainly I, would have alone—the awards we've won together speak for themselves.
A large part of why I have gotten better at practicing mathematics over the years—a skill of fundamental importance—is thanks to working with you.

I thank James T. Wilson, Steindór Sæmundsson, Samuel Cohen, and So Takao, for the opportunity to do research and write together.
I am grateful to Yicheng Luo, Sanket Kamthe, Hugh Salimbeni, Yasemin Bekiroğlu, Sicelukwanda Zwane, K. S. Sesh Kumar, and others in the group who have spoken to me at length about ideas during my time.
I feel fortunate to to have been your labmate, and have learned an immense amount from you over the years.

I thank Andreas Hochlehnert, Kai Biegun, and Lucas Cosier for trusting in me as a researcher enough to allow me to supervise your work.
Watching your ideas develop and grow has been a beautiful experience, and I look forward to seeing you shine in the years to come.

I am grateful to Michael Hutchinson, James-Michael Leahy, Peter Mostowsky, Iskander Azangulov, and others who I have collaborated with on papers, or spoken to at length about ideas.
I am particularly grateful to Brandon Amos both for the opportunity to collaborate, and for suggesting I speak to Marc Deisenroth in the weeks leading up to me becoming his student—without this, my path would have been very different.
I thank Nick Sharp for teaching me how to make the three-dimensional figures in this thesis.

I am grateful to Imperial College, and the Department of Mathematics in particular, for both the opportunity in funding my studies, and for support when things did not go as planned.
I am particularly thankful to Henrik Jensen: thank you for treating me with kindness in difficult times.
I am grateful to Pierre Degond, Darryl Holm, and Chris Hallsworth for your ideas in teaching me mathematics.
I joined Imperial because I wanted to get better at mathematics, and it is thanks largely to you that I feel I succeeded at this.

I am thankful to Måns Magnusson, Leif Jonsson, and Shawfeng Dong for collaborating with me in my early days as a researcher, and helping me learn and get my bearings together.
In particular, I thank Peter Drake and Raya Feldman for first setting me on this path.
I am grateful to Matthew Johnson and Chris De Sa for the counterexamples sent to me in the early times, which convinced me that proper mathematics was of utmost importance and that I had to improve at it in order trust myself to say things that are true.

I am grateful to my friends at Carnegie Mellon University, including Dominic Chen, Aurick Qiao, Willie Neiswanger, Kumail Jaffer, Ziv Scully, Sarah Allen Scully, Sol Boucher, Guillaume Didier, Stefan Muller, Priya Donti, Gabriele Farina, Noam Brown, Ben Blum, Evan Cavallo, and others who convinced me to apply to a PhD program by showing me that many that of the people spending their life studying ideas were just like me.
I am grateful to Yuanran Zhu for showing me the same, but in a different place and time.

I am grateful to friends who have supported me over the years, including to Jon Frost for almost a decade of friendship, as well as Paula Siauchó Unriza, Victor Espinosa, Sam Aragon, and others who know me well and have been there for me over the years.
Even as time had gone on, whether we last spoke days or years ago, I have always found nothing about our friendship to have changed whenever we have had the chance to speak again.

Finally, I want to thank my family, including my mom, Irina Terenina.
Without your sacrifices, I would not have spent my youth in the United States, or had almost any of the opportunities I have had.
I am grateful to my grandmother, Olga Novikova, for your wisdom and for keeping our family together.
I wish you could have seen me graduate and write this thesis.

I am grateful to my dad, Vadim Terenin: your authoritative presentation style has at times blinded me to your ideas.
I am grateful to my stepmom, Rheesa Eddings, for both your deeply insightful way of understanding the world, and for your humanity, diplomacy, kindness, helpful ideas, and support over the years.

My gratitude to family includes those who are family by virtue of their friendship and support for my entire adult life, including Jerry Hirsch, John Halper, Lois Morera, Christian Morera, and Peter Morera.
Thank you for continuing to stay in touch, for being the extended family I never otherwise had, and for being with me as my journey has unfolded.

Though this list of acknowledgments has ballooned beyond what might otherwise be expected in a dissertation, I nonetheless feel compelled to write it in full by virtue of my path being what it was.
It is also without a doubt incomplete, so I am grateful to anyone whose name should have been here but has been accidentally omitted, including whoever suggested I limit all paragraphs written to seven lines or less.
To conclude, I offer you, the reader, my exceeding gratitude for taking the time to read my work and ideas.

CHAPTER: ABSTRACT

Bayesian learning using Gaussian processes provides a foundational framework for making decisions in a manner that balances what is known with what could be learned by gathering data.
In this dissertation, we develop techniques for broadening the applicability of Gaussian processes.
This is done in two ways.

Firstly, we develop pathwise conditioning techniques for Gaussian processes, which allow one to express posterior random functions as prior random functions plus a dependent update term.
We introduce a wide class of efficient approximations built from this viewpoint, which can be randomly sampled once in advance, and evaluated at arbitrary locations without any subsequent stochasticity.
This key property improves efficiency and makes it simpler to deploy Gaussian process models in decision-making settings.

Secondly, we develop a collection of Gaussian process models over non-Euclidean spaces, including Riemannian manifolds and graphs.
We derive fully constructive expressions for the covariance kernels of scalar-valued Gaussian processes on Riemannian manifolds and graphs.
Building on these ideas, we describe a formalism for defining vector-valued Gaussian processes on Riemannian manifolds.
The introduced techniques allow all of these models to be trained using standard computational methods.

In total, these contributions make Gaussian processes easier to work with and allow them to be used within a wider class of domains in an effective and principled manner.
This, in turn, makes it possible to potentially apply Gaussian processes to novel decision-making settings.

CHAPTER: PRINCIPAL NOTATION

Principal Notation

§ GENERAL

@ p1.625cm X<-1pc[1pc][r][0pc].0pc @ >p0.625cm @
$(\Omega,\c{F},\P)$     Probability space     ntn:prob

$\E(\.)$     Expectation of a real-valued random variable     ntn:expectation-cov

$\Cov(\.,\.)$     Covariance between two real-valued random variables     ntn:expectation-cov

$y \~ \pi$     Random variable $y$ with distribution $\pi$     ntn:rv

$\c{M}_1(X)$     Space of probability measures over $X$     ntn:space-of-measures

$D_{\f{KL}}(\.||\.)$     Kullback–Leibler divergence between probability measures     ntn:kl-div

$W_{p,d}(\.,\.)$     Wasserstein distance of $p$th order over distance $d$     ntn:wass-dist

$\N$     Natural numbers, not including zero

$\Z$     Integers

$\R$     Real numbers

$\C$     Complex numbers

$\R^d$     Euclidean space of dimension $d$

$\R^X$     Space of functions from $X$ to $\R$

$\v{a}$     Euclidean vector (bold italic letters)

$\m{A}$     Matrix (bold upface letters)

$\oplus$     Direct sum of vector spaces     ntn:direct-sum

$V^*$     Dual space of a topological vector space     ntn:dual-space

$\c{A}^*$     Adjoint of a bounded linear operator     ntn:bounded-operators

$\norm{\.}$     Norm of a vector in a Banach space     ntn:banach-hilbert

$\innerprod{\.}{\.}$     Inner product between vectors in a Hilbert space     ntn:banach-hilbert

$\dualprod{\.}{\.}$     Duality pairing between topological vector spaces     ntn:dual-pair

$L(V;W)$     Banach space of bounded linear operators     ntn:bounded-operators

$C^0(X;\R)$     Banach space of real-valued continuous functions     ntn:conts-fns

$L^p(X;\R)$     Lebesgue space     ntn:leb-space

§ BAYESIAN LEARNING

@ p1.625cm X<-1pc[1pc][r][0pc].0pc @ >p0.625cm @
$(\Theta,\mathit\Theta)$     Measurable space representing quantity of interest     ntn:meas

$(Y,\c{Y})$     Measurable space representing the data     ntn:meas

$A_{(\.)}$     Measurable set     ntn:meas

$\theta, y$     Random variables     ntn:rv

$\pi_{(\.)}$     Probability measure     ntn:prob

$f_{(\.)}$     Probability density     ntn:density

$(\.\given\.)$     Jointly measurable stochastic process     ntn:jmsp

$\pi_{(\. \given \.)}$     Probability kernel     ntn:prob-ker

$f_{(\.\given\.)}$     Conditional probability density     ntn:density

§ MARKOV DECISION PROCESSES

@ p1.625cm X<-1pc[1pc][r][0pc].0pc @ >p0.625cm @
$S$     State space     ntn:mdp

$A$     Action space     ntn:mdp

$r(\.\given\.,\.)$     Reward kernel     ntn:mdp

$p(\.\given\.,\.)$     Transition kernel     ntn:mdp

$\pi$     Policy     ntn:mdp-policy

$V^{(\pi)}$     Value function with respect to policy $\pi$     ntn:mdp-value

$R^{(\pi)}$     Regret of policy $\pi$     ntn:mdp-regret

§ MULTI-ARMED BANDITS

@ p1.625cm X<-1pc[1pc][r][0pc].0pc @ >p0.625cm @
$X$     Space of arms     ntn:mab

$K$     Number of arms     ntn:mab-number-arms

$T$     Maximum time     ntn:mab-regret

$y$     Observed stochastic rewards     ntn:mab

$f(\.)$     Expected rewards     ntn:mab

$\eps$     Observation noise     ntn:mab

$R$     Cumulative stochastic regret of a bandit algorithm     ntn:mab-regret

$\Delta(\.)$     Regret of a specific arm     ntn:mab-regret-arms

$n_t$     Cumulative number of arm pulls up to time $t$     ntn:mab-regret-arms

$f^+_t(\.)$     Upper confidence bound at time $t$     ntn:mab-ucb

$\mu_t(\.)$     Empirical mean rewards at time $t$     ntn:mab-means

$\sigma_t(\.)$     Width of confidence bounds at time $t$     ntn:mab-ucb

$c_t$     Scaling constant for confidence bound width at time $t$     ntn:mab-ucb

§ GAUSSIAN PROCESSES

@ p1.625cm X<-1pc[1pc][r][0pc].0pc @ >p0.625cm @
$\f{N}(\mu,\sigma^2)$     Gaussian distribution with mean $\mu$ and variance $\sigma^2$     ntn:norm-dist

$\f{N}(\v\mu,\m\Sigma)$     Multivariate Gaussian with mean $\v\mu$ and covariance $\m\Sigma$     ntn:mvn-dist

$k(x,x')$     Positive semi-definite kernel     ntn:kernel

$\m{K}_{\v{x}\v{x}'}$     Kernel matrix with entries $k(x_i, x'_j)$     ntn:kernel

$\f{GP}(\mu,k)$     Gaussian process with mean $\mu$ and covariance kernel $k$     ntn:gp-dist

$f$     Prior Gaussian process     ntn:gp-model

$\v\eps$     Likelihood noise     ntn:gp-model

$\m\Sigma$     Covariance of likelihood noise     ntn:gp-model

$(\v{x},\v\gamma)$     Data used to condition the model $y_i = f(x_i) + \eps_i$      ntn:gp-model

$f\given\v{y}$     Posterior Gaussian process conditioned on $\v{y} = \v\gamma$     ntn:gp-model

$n$     Number of training data points     ntn:gp-complexity

$n_*$     Number of evaluation locations     ntn:gp-complexity

$\tilde{f}$     Approximate prior defined via a finite basis expansion     ntn:gp-approx-prior

$\phi_i(\.)$     Basis function used by approximate prior     ntn:gp-approx-prior

$w_i$     Random weight of basis function $\phi_i$     ntn:gp-approx-prior

$\ell$     Number of basis functions used in approximate prior     ntn:gp-approx-prior

$k(x_j,\.)$     Canonical basis function     ntn:canonical-basis-fns

$v_j$     Random weight of canonical basis function $k(x_j,\.)$     ntn:canonical-basis-fns

$H_k$     Reproducing kernel Hilbert space with kernel $k$     ntn:rkhs

$\v\phi$     Finite-dimensional approximate feature map     ntn:approx-feature-map

$m$     Number of inducing points in sparse approximation     ntn:gp-num-inducing

§ MANIFOLDS AND GRAPHS

@ p1.625cm X<-1pc[1pc][r][0pc].0pc @ >p0.625cm @
$X$     Manifold     ntn:manifold

$TX$     Tangent bundle     ntn:tangent-bdl

$T^*X$     Cotangent bundle     ntn:cotangent-bdl

$C^\infty(X)$     Space of smooth real-valued functions     ntn:smooth-fns

$\Gamma(TX)$     Space of smooth sections     ntn:smooth-sections

$d$     Dimension of $X$     ntn:manifold

$f_*$     Pushforward of a function $f$     ntn:vector-pushforward

$g(\.,\.)$     Metric tensor of a Riemannian manifold     ntn:metric-tensor

$\lap_g$     Laplace–Beltrami operator induced by $g$     ntn:laplace-beltrami

$(\lambda_n, f_n)$     Laplace–Beltrami eigenpair     ntn:laplace-beltrami-eigenpairs

$\Phi(\lap_g)$     Operator defined via functional calculus     ntn:functional-calculus

$\c{W}_g$     Gaussian white noise induced by $g$     ntn:riemannian-white-noise

$\m\lap$     Graph Laplacian     ntn:graph-laplacian

$\bc{W}$     Standard Gaussian     ntn:graph-white-noise

$\Phi(\m\lap)$     Matrix defined via functional calculus     ntn:matrix-functional-calculus

$F$     Frame over $X$     ntn:frame

$\m{P}_F(\.)$     Projection map induced by $F$     ntn:frame-projection

CHAPTER: INTRODUCTION

Learning from experience in order to change behavior is one of the defining abilities of biological systems, which differentiates them from other kinds of systems found in the world.
Replicating the processes biological systems use to learn and adapt is a fundamental goal of science and technology.
To this end, the development of mathematical formalisms rich enough to capture the notion of learning is one of the crowning achievements of statistics, machine learning, and artificial intelligence.

One such formalism is the Bayesian view of learning.
The idea behind Bayesian learning is to represent the information known about the quantity of interest using probability.
To do so, the relationship between the quantity of interest and the data is formalized as a joint probability distribution.
This gives rise to a conditional probability distribution describing what was learned about the quantity of interest by observing the data.

Bayesian learning fits naturally within a theory of decision, which describes how an abstract decision system should select actions in pursuit of a goal.
This is done by learning how different actions affect pursuit of the goal, and selecting optimal actions consistent with what was learned.
By virtue of being probabilistic, such decision systems assess and propagate uncertainty, enabling them to balance what is already known with what could be learned by taking actions—a concept known as the explore-exploit tradeoff.

The performance of a decision system can be evaluated by examining how quickly its decisions improve and become optimal.
A decision system's regret is the reduction in its quality of decisions by virtue of not knowing the quantity of interest in advance.
In most non-trivial settings, one can show that some regret is inevitable: a decision-making system must make some degree of mistakes in order to learn.
A decision system is considered optimal if its regret is within a constant factor of the best possible regret.

Decision systems with optimal or close-to-optimal regret require less data in order to solve their respective tasks, and are called data-efficient.
Data-efficiency is a key concern in practical settings, where data-collection takes time and can be expensive.
By virtue of resolving explore-exploit tradeoffs in a manner amenable to regret analysis, the Bayesian formalism gives broad tools for constructing data-efficient decision systems.

The key limitation of the Bayesian approach is that it often leads to computational problems which are intractable.
Conditional distributions generally contain more information than actually needed to make optimal decisions, yet calculating them is largely unavoidable.
Probabilistic decision systems are thus most attractive in settings where their strengths—including data-efficiency, solid technical foundations, and amenability to analysis—can shine, while computational costs are kept under control.

In my view, Gaussian processes are one such setting: they are powerful enough to model wide classes of unknown quantities of interest, yet their computational costs are generally polynomial.
Better yet, Gaussian-process-based decision systems have demonstrated excellent performance in real-world scientific applications.
Studying Gaussian processes is therefore a promising avenue towards improved understanding of Bayesian learning and Bayesian decision-making in pursuit of artificial intelligence.

The goals of this dissertation are twofold: (i) to make Gaussian processes easier to work with when used within larger decision systems, and (ii) to expand the set of settings where Gaussian processes can used, enabling construction of decision systems for applications not previously considered.
Contributions toward (i) include path-wise conditioning techniques studied in <Ref>, and contributions toward (ii) include non-Euclidean Gaussian processes studied in <Ref>.
Following these, <Ref> concludes.

To pursue these goals, it is critically important that all of the concepts described in the preceding paragraphs be made into rigorous mathematics, so that the ideas described in the sequel ultimately reduce to definitions and implications, and not metaphor or opinion.
Together, we therefore begin by defining the key mathematical notions needed.

§.§ Contributions

The work presented in this thesis is published as a series of papers.
My contributions to the individual works are primarily on the theoretical and methodological side, and are described below.

In <Ref>, we present pathwise conditioning techniques: this work is published as wilson20,wilson21.
My contributions include (i) development of the random-function-based formalism for describing pathwise conditioning of Gaussian processes, (ii) error analysis of basis-function-approximation-based pathwise sampling, (iii) re-interpretation of inducing point methods, and (iv) review of prior sampling methods.
All of these were developed jointly with the other authors.

In <Ref>, we present Matérn Gaussian processes in non-Euclidean settings: this work is published as borovitskiy20,borovitskiy21,jacquier21,hutchinson21.
My contributions include (i) developing the differential-geometric and stochastic-partial-differential-equation-based formalisms for defining these processes, and (ii) describing computational techniques for working with these models in practice.
All of these ideas were developed jointly with the other authors.

§ BAYESIAN LEARNING

The first concept we develop in depth along our path towards a mathematically precise understanding of statistical decision-making is Bayesian learning—a mathematical formalism for reasoning about unknown quantities of interest on the basis of data.
Bayesian learning is a probabilistic theory: a model specifies how the quantity of interest and the data depend on one another within a probability distribution.
Learning entails calculating how the distribution of the quantity of interest changes upon observing the data.

One of the key strengths of Bayesian theory is that it applies in wide generality, owing to the substantial scope of probability theory.
In particular, one can study learning of quantities of interest that are function-valued using observed data consisting of pointwise function evaluations.
This, however, requires a non-elementary treatment, as one cannot rely solely on probability densities to define what conditional distributions are.
We therefore begin by recalling the necessary mathematical formalism.

For an overview of Bayesian methods from a model-building perspective, see gelman14, and for an overview from a mathematical perspective see ghosal17,gine15.

§.§ Review of probability theory

We adopt the language of measure-theoretic probability, which we now describe.
To ease presentation, we state the definitions together with useful ways of thinking about them.
We use kallenberg06 as our standard reference, along with certain results from bogachev07a,bogachev07b.

We say that measurable space is a pair $(Y,\c{Y})$ consisting of a set $Y$ and a $\sigma$-algebra $\c{Y}$ over $Y$.
A $\sigma$-algebra is a set of subsets of $Y$ containing the space itself which is closed under countable unions, intersections, complements, and therefore set-theoretic monotone limits.
These can be reinterpreted as Boolean logical operations, so $\c{Y}$ can be thought of as the set of all true-false questions one can ask about elements of the set $Y$.
These questions, then, are closed under and/or/not operations and monotone limits thereof.

Product of measurable spaces
Given two measurable spaces $(Y,\c{Y})$ and $(\Theta,\mathit\Theta)$, if we form the Cartesian product $Y \x \Theta$, then we can define the product $\sigma$-algebra $\c{Y}\ox\mathit\Theta$ as the smallest $\sigma$-algebra containing all sets of the form $A_y \x A_\theta$ with $A_y \in\c{Y}$ and $A_\theta\in\mathit\Theta$.
Measurable subspace
For a measurable subset $Y' \subseteq Y$, we can define $\c{Y}' = \{A_y \^ Y' : A_y \in\c{Y}\}$, which is also a $\sigma$-algebra, thereby making $(Y',\c{Y}')$ into a measurable space.
We call $\c{Y}'$ the subset $\sigma$-algebra.

If $Y$ is a topological space, the Borel $\sigma$-algebra $\c{B}(Y)$ is defined as the smallest $\sigma$-algebra containing the open sets in the topology.
Topological spaces, in turn, are sets equipped with additional structure enabling them to admit notions such as locality and convergence—we review these spaces in additional detail in <Ref>.

Measurable function
A map $f : Y \-> Y'$ between measurable spaces $(Y,\c{Y})$ and $(Y',\c{Y}')$ is said to be measurable if its preimage defines a map $f^{-1} : \c{Y}' \-> \c{Y}$ between the respective $\sigma$-algebras.
This intuitively means that true-false questions for the space $Y'$ can be asked and answered relative to those in $Y$.
On product spaces, a map $f : Y \-> \Theta \x \Theta'$ is measurable if its components are measurable, but a map $f : Y \x Y' \-> \Theta$ need not automatically be measurable if $f(\.,y') : Y \-> \Theta$ and $f(y,\.) : Y' \-> \Theta$ are measurable for all $y$ and $y'$.

A probability measure is a non-negative countably additive map $\pi_y : \c{Y} \-> \R$ satisfying $\pi_Y(Y) = 1$.
This can be thought of as a map that takes a true/false question, and assigns a number indicating how close to true or false its answer is according to the measure—in this view, probability measures describe uncertainty.
A probability space $(\Omega,\c{F},\P)$ consists of a measurable space and probability measure.
The set $\Omega$ can be thought of as a space of abstract random numbers, with a random number generator described by $\P$.

Given a probability space, we say that a random variable is a measurable map $y : \Omega \-> Y$.
A random variable, then, maps random numbers $\omega\in\Omega$ into the space $Y$.
The distributionDistribution of a random variable is defined as the pushforward measure $\pi_y = y_* \P$ where $y_*$ is defined as $(y_* \P)(A_y) = \P(y^{-1}(A_y))$ for all $A_y\in\c{Y}$.
The probability of an event $A_y$, then, is determined by measuring the probability of random numbers under which $A_y$ occurs—measurability guarantees this is possible.

Notation for random variables
In this work, we will generally not adopt the standard convention of suppressing $\omega$-arguments of random variables from our notation, and will write such arguments explicitly in cases that other function arguments are also used.
Though this makes expressions slightly denser, it also avoids ambiguity and presents the mathematics more precisely.

The expectation $\E(y)$ of a real-valued random variable $y$, if it exists, is defined as its integral with respect to $\P$.
This can be thought of as the average value of $y$ under the random numbers generated from $\P$, with averaging performed via the algebraic structure of the reals.
This notion extends to vector-valued random variables.
The covariance of two real-valued random variables, if it exists, is defined as $\Cov(y,y') = \E(yy') - \E(y)\E(y')$.
From this, we define the covariance of a finite-dimensional random vector componentwise.

Equality of random variables
There are multiple senses in which we can say random variables are equal.
We say that $y = y'$ surely if they are equal as mathematical functions.
This notion is essentially never used, because it is exceedingly strong.
We say that $y = y'$ almost surely if $\P(y \neq y') = 0$, in which case $y$ and $y'$ are equal as functions, except possibly on a set of probability zero.
We say that $y = y'$ in distribution if $y_* \P = y'_* \P$, meaning their distributions are equal.

Existence of a random variable with given distribution
For a given probability space $(\Omega,\c{F},\P)$, a random variable $y : \Omega \-> Y$ with distribution $\pi_y$ need not exist.
However, if it does not, there always exists another probability space $(\Omega',\c{F}',\P')$ and a measurable map $i : \Omega' \-> \Omega$ such that $\P = i_* \P'$, and for which a $\pi_y$-distributed random variable $y : \Omega' \-> Y$ does exist.
We thus implicitly assume all probability spaces are large enough to ensure all random variables with prescribed distributions exist.

Probability kernel
We will be interested in probability measures extended to allow them to be parameterized by other quantities.
A probability kernel is defined to be a map $\pi_{y\given\theta} : \c{Y} \x \Theta \-> \R$ satisfying two conditions: (i) the map $\pi_{y\given\theta}(\.,\vartheta) : \c{Y} \-> \R$ is a probability measure for all $\vartheta\in\Theta$, and (ii) the map $\pi_{y\given\theta}(A_y,\.) : \Theta \-> \R$ is measurable for all $A_y\in\c{Y}$.
In particular, this means that a probability kernel can be integrated against other measures.

Conditional distribution
Random variables can also be extended to parameterize them by other quantities.
A $\c{F}\ox \mathit\Theta$-measurable map $y\given\theta : \Omega \x \Theta \-> Y$ is called a jointly measurable stochastic process—we largely eschew this terminology to emphasize the Bayesian formalism.
In particular, unlike most presentations of stochastic processes, we never think of $\Theta$ as representing time.
Define the conditional distribution of $y\given\theta$ to be $(y\given\theta)_*\P$, with the pushforward taken in the first argument—by <Ref>, this is a probability kernel.

§.§ Bayes' Rule for probability measures

The key idea behind Bayesian learning is to formalize learning using the concept of conditional probability.
This entails (i) building a model quantifying the relationship between the quantity of interest $\theta$, and the data $y$, and (ii) quantifying what is learned using Bayes' Rule.
We now begin to make these concepts precise in the language of measure-theoretic probability.

Let $Y$ and $\Theta$ be sets.
A Bayesian model is a probability measure defined on $\Theta \x Y$.

A model can be constructed in a number of different ways.
The most common technique is to specify two components: (i) a probability distribution describing what is known about the quantity of interest external to the data, and (ii) how the data relates to the quantity of interest.
These are called the prior distribution and likelihood, respectively.
For example, a simple Gaussian model is given by

\[
y_i \given \theta &\~[N](\theta, 1)
\theta &\~[N](0,1)
\]

In this model, $\theta$ is an unknown scalar assigned a standard normal prior, while $y_i\given\theta$ is assigned a Gaussian likelihood centered at $\theta$.
Our first goal is to make this notation into proper mathematics.

A Bayesian model can be constructed in the following ways.
Using measures: integrate the following two components.
Prior: a probability measure $\pi_\theta : \mathit\Theta \-> \R$.
Likelihood: a probability kernel $\pi_{y\given\theta} : \c{Y}\x\Theta \-> \R$.

Using random variables: compose the following two components.
Prior: a random variable $\theta:\Omega\->\Theta$.
Likelihood: a jointly measurable stochastic process $y\given\theta: \Omega\x\Theta\->Y$.

Moreover, if one takes the pushforward of the composition of the obtained random variables with respect to $\P\x\P$, then the two constructions coincide.

For the former, define

\[
\pi_{\theta,y}(A_\theta\x A_y) = \int_{A_\theta} \pi_{y\given\theta}(A_y\given\theta) \d\pi_\theta(\theta)
\]

which extends to the full product $\sigma$-algebra by <Ref>, giving the desired probability measure.
For the latter, define the random variable

\[
\gamma : \Omega\x\Omega &\-> \Theta \x Y
\gamma : (\omega,\omega') &\|> (\theta(\omega), (y\given\theta)(\omega',\theta(\omega)))
\]

and take $\pi_{\theta,y} = \gamma_* (\P\x\P)$.
We now prove these expressions coincide.

\[
\pi_{\theta,y}(A_\theta\x A_y) &= \int_{A_\theta} \pi_{y\given\theta}(A_y\given\theta) \d\pi_\theta(\theta)
\\
&= \int_\Omega\1_{\theta(\omega)\in A_\theta} \pi_{y\given\theta}(A_y\given\theta(\omega)) \d\bb{P}(\omega)
\\
&= \int_\Omega\1_{\theta(\omega)\in A_\theta} \int_\Omega \1_{(y\given\theta)(\omega',\theta(\omega)) \in A_y} \d\bb{P}(\omega') \d\bb{P}(\omega)
\\
&= \int_\Omega\int_\Omega \1_{(\theta(\omega),(y\given\theta)(\omega',\theta(\omega))) \in A_\theta \x A_y} \d\bb{P}(\omega') \d\bb{P}(\omega)
\\
&= (\P\x\P)(\gamma \in A_\theta \x A_y)
\]

which follows using Tonelli's Theorem and <Ref>.

These two components should be interpreted as follows: $\theta$ describes what is known about the quantity of interest external to the data, and $y\given\theta$ describes how the data $y$ relates to the quantity of interest.
More specifically, the likelihood describes how $y$ would be distributed if $\theta$ was known and equal to the conditioned value.

The condition that the pushforward of the composition of the prior and likelihood needs to be taken with respect to the product measure $\P\x\P$ should be interpreted as a kind of non-duplicate dependence condition which ensures that $y\given\theta$ depends on $\theta$ only through its second argument, and not through the abstract random numbers $\omega$.

Given a Bayesian model, we can formalize the notion of what is learned about $\theta$ from observing $y$ as its conditional probability distribution—note that this is meant in a distributional sense, not a random-variable-based sense.
The formulation we use requires topological assumptions: we say that a topological space is Polish if it is homeomorphic to a complete separable metric space.
Note that any two uncountable Polish spaces are Borel-isomorphic: see [Chapter 1]villani08, for discussion.
The key result is given as follows.

external/export next=false

[axis lines=none, height=7.5cm, axis equal, view=31022.5, xmin=-0.6, xmax=1.75, ymin=-1.75, ymax=1.75, zmin=0, zmax=1]
at (0,0.03,0.31) ;
3[no markers, domain=-1.75:1.75, smooth, samples=32, samples y=0, thick, color=black, opacity=0] (-0.5,x,0.9 * exp(-x^2) * exp(-(x + 0.5)^2));
3[no markers, very thick, color=black, line cap=round, dashed] coordinates (-0.5,-1.75,0) (-0.5,1.75,0);

[Bayesian models]Here, we illustrate a Bayesian model. For all possible values of the quantity of interest, the likelihood describes the respective distribution of the data, shown on the left. This is combined with the prior, also shown on the left, to form a joint distribution, shown on the right. From this, the posterior distribution, also shown on the right, is obtained by conditioning. Note that the two plots are shown not to scale, and that the joint and posterior live on different spaces, hence are normalized differently.

[Bayes' Rule]
Suppose that $\Theta$ is a Polish topological space, and let $\pi_y = \pi_{\theta,y}(\Theta\x\.)$.
Then for every Bayesian model $\pi_{\theta,y}$ there is a $\pi_y\ae[-]$ unique probability kernel $\pi_{\theta\given y}$ satisfying

\[
\pi_{\theta,y}(A_\theta \x A_y) = \int_{A_y} \pi_{\theta\given y}(A_\theta \given y) \d\pi_y(y)
\]

which we call the posterior distribution.

The claim follows directly from [Corollary 10.4.15]bogachev07b.
See also [Theorem 5.3 and Theorem 5.4]kallenberg06, [Theorem 5.3.1]ambrosio08, and chang97 for variations of this result.

This result is illustrated in <Ref>, and shows that given a Bayesian model, the posterior distribution describing what was learned from the data exists.
In its full abstract formulation, however, Bayes' Rule is non-constructive, and it is not at all clear how to calculate any kind of useful formula from it.
Moreover, the null sets can be problematic: to ensure that $\pi_{\theta\given y}$ is defined pointwise, one needs further properties such as continuity.
Fortunately, in many settings these issues are resolved by virtue of additional structure.

The simplest such structure occurs when $\pi_{\theta,y}$ admits a density with respect to a product measure $\lambda$.
In this case, it follows from <Ref> that $\pi_\theta$ admits the density $f_\theta$ with respect to $\lambda(\.\x Y)$, and similarly for $\pi_y$ and $f_y$.
Define a conditional density as the ratio of joint and marginal densities, for instance $f_{\theta\given y} = \frac{f_{\theta,y}}{f_y}$, and let $\propto$ denote equality up to a multiplicative proportionality constant.
Then, we have the following.

Suppose that $\pi_{\theta,y}$ admits the density $f_{\theta,y}$ with respect to $\lambda_{\theta,y}$.
Then we have

\[
f_{\theta\given y} \propto f_{y\given\theta}f_\theta
\]

We have

\[
f_{\theta,y} &= \frac{f_{\theta,y}}{f_\theta} f_\theta = f_{y\given\theta}f_\theta
f_{\theta,y} &= \frac{f_{\theta,y}}{f_y} f_y = f_{\theta\given y}f_y
\]

which, when combined, give the result.

Sometimes, the knowledge of $f_{\theta\given y}$ up to proportionality is enough to fully deduce its form.
For example, if the likelihood is Gaussian with unknown mean and known variance, and the prior is Gaussian, then posterior is also Gaussian.
In cases such as this, where the prior and posterior land in the same parameterized class of distributions, the pair $(\pi_{y\given\theta}, \pi_\theta)$ are called conjugate.

Bayes' Rule for densities is remarkable in its generality yet restrictiveness.
On the one hand, we require no direct assumptions about the spaces $\Theta$ and $Y$, and in particular allow for real spaces, discrete spaces, and Riemannian manifolds.
On the other hand, other than in the aforementioned settings, where we can employ the Lebesgue, counting, and Riemannian volume measures, respectively, finding a suitable reference measure can be difficult.

We will work with posterior distributions in infinite-dimensional vector spaces in the sequel—there, densities are either not available or not convenient.
In those cases, it is enough to calculate the posterior on arbitrary finite-dimensional marginal projections to uniquely determine its value on the full infinite-dimensional space.

Conditioning and marginalization commute: given a probability measure $\pi_{\theta,\theta',y}$, if Bayes' Rule holds, then we have $\pi_y\ae[-]$ that

\[
\pi_{\theta \given y} = \pi_{\theta,\theta' \given y}(\. \x \Theta')
\]

Let $\pi_{\theta,\theta'\given y}$ be the $\pi_y\ae[-]$ unique probability kernel given by the Disintegration Theorem satisfying

\[
\pi_{\theta,\theta',y}(A_{\theta,\theta'} \x A_y) = \int_{A_y} \pi_{\theta,\theta'\given y}(A_{\theta,\theta'} \given y) \d\pi_y(y)
\]

Plugging in $A_{\theta} \x \Theta'$ for $A_{\theta,\theta'}$ gives

\[
\pi_{\theta,y}(A_\theta \x A_y) = \int_{A_y} \pi_{\theta,\theta'\given y}(A_\theta \x \Theta' \given y) \d\pi_y(y)
\]

On the other hand, applying the Disintegration Theorem to $\pi_{\theta,y}$ gives a $\pi_y\ae[-]$ unique probability kernel satisfying

\[
\pi_{\theta,y}(A_\theta \x A_y) = \int_{A_y} \pi_{\theta\given y}(A_\theta \given y) \d\pi_y(y)
\]

Since both $\pi_{\theta,\theta'\given y}$ and $\pi_{\theta,y}$ are defined with respect to the same null sets, we conclude by uniqueness that they coincide, and the claim follows.

This result makes densities into a substantially more powerful tool than they would be otherwise, since it enables us to use them for calculating posterior distributions even where they are not directly available.
In particular, we can map an infinite-dimensional function space into finite-dimensional vector spaces induced by pointwise function evaluations at arbitrary points, enabling us to calculate posterior distributions in such settings, in spite of no suitable densities existing directly in the space of interest.

We now introduce the variational formulation of Bayes' Rule, which expresses the posterior as the solution to an infinite-dimensional optimization problem.
Let $\c{M}_1(\Theta)$ be the space of all probability measures over $\Theta$.

Let $\pi_{y,\theta}$ be a Bayesian model, and assume it admits the density $f_{\theta,y}$ with respect to a product measure $\lambda$.
Then for every $\gamma\in Y$, the posterior distribution satisfies

\[
\pi_{\theta\given y}(\.\given \gamma) = \argmin_{\bb{q}_\theta \in \c{M}_1(\Theta)} D_{\f{KL}}(\bb{q}_\theta \from \pi_\theta) - \operatorname*{\E}_{\vartheta\~\bb{q}_\theta} \ln f_{y\given\theta}(\gamma\given\vartheta)
\]

where $D_{\f{KL}}$ is the Kullback–Leibler divergence between probability measures, and the minima does not depend on the choice of $\lambda$.

Write the quantity being minimized as

\[
D_{\f{KL}}(\bb{q}_\theta \from \pi_\theta) - \operatorname*{\E}_{\vartheta\~\bb{q}_\theta} \ln f_{y\given\theta}(\gamma\given\vartheta) &= \operatorname*{\E}_{\vartheta\~\bb{q}_\theta} \ln \frac{f_{\bb{q}}(\vartheta)}{f_\theta(\vartheta) f_{y\given\theta}(\gamma\given\vartheta)}
\\
&= \operatorname*{\E}_{\vartheta\~\bb{q}_\theta} \ln \frac{f_{\bb{q}}(\vartheta)}{f_{\theta\given y}(\vartheta\given\gamma) f_y(\gamma)}
\\
&= D_{\f{KL}}(\bb{q}_\theta \from \pi_{\theta\given y}(\.\given \gamma)) - \ln f_y(\gamma)
\]

Since $\ln f_y(\gamma)$ does not depend on $\bb{q}_\theta$, and since $D_{\f{KL}}(\mu \from \nu) = 0$ implies $\mu = \nu$, we conclude that the objective is minimized by taking $\bb{q}_\theta = \pi_{\theta\given y}(\.\given\gamma)$.
Since the minima does not depend on the choice of the measure $\lambda$ with respect to which the densities are defined, the claim follows.

For any given dataset, this result shows that  Bayes' Rule can be viewed in information-theoretic manner: among all probability measures, the posterior maximizes predictive power, while retaining as many bits as possible from the prior, in the sense given by the Kullback–Leibler divergence.

It's also possible to prove the result in a different way using the calculus of variations, with variations taken in the Banach space of signed measures under the total variation norm.
That argument also employs Bayes' Rule for densities for its key step, making it largely similar in spirit.

The result suggests a way to approximately compute posterior distributions: restrict the optimization problem to a suitably chosen subspace of the space of all probability measures $\c{M}_1(\Theta)$.
This strategy will be particularly fruitful in the later-described setting of Gaussian processes, where techniques for constructing such approximations with well-understood and favorable accuracy will be considered.

We always view variational approximations as minimization of Kullback–Leibler divergences, as this is mathematically sound.
We will eschew standard presentations involving evidence lower bounds: the only widely-known mathematical explanations for why maximizing these bounds should improve model performance appeal to Kullback–Leibler divergences—so, then, why talk about evidence lower bounds in the first place?

Once a posterior is calculated, the next step is to extract the relevant quantities from it.
Traditionally, this is often done by displaying summary statistics to be evaluated and interpreted by a person with statistical training, often with a focus on assessing uncertainty.
We will instead focus on settings where the posterior is given as input to an upstream decision-making algorithm: we explore these next.

§.§ Technical lemmas

We now prove a number of technical lemmas used in the preceding text, which for completeness are presented here in order to avoid disrupting the reader's flow.

Let $b:\Omega\x A\->B$ be a jointly measurable stochastic process.
Then the map $b_*\P : \c{B} \x A \-> \R$, where the pushforward is taken in the first argument, is a probability kernel.

It is clear that $b_* \P$ is a probability measure for all $a'\in A$, so we only need to prove that the map $a \|> (b(\.,a)_*\P)(A_b)$ is measurable for all $A_b\in\c{B}$.
First, write

\[
(b_*\P)(A_b) = \P(b(\.,a)^{-1}(A_b)) = \int_\Omega \1_{b(\omega,a)\in A_b} \d\bb{P}(\omega)
\]

Now, note that the map $\Omega\x A\->\R$ given by $(\omega,a) \|> \1_{b(\omega,a)\in A_b}$ is bounded measurable, since $b$ is measurable in both arguments, and indicators of measurable functions on measurable sets are bounded measurable.
Finally, since for any bounded measurable $f : \Omega \x A \-> \R$, the map $a \|> \int_\Omega f(\omega,a) \d\bb{P}(\omega)$ is measurable by Fubini's Theorem, the claim follows.

Let $\c{S}$ be a semi-algebra of sets, which is defined as a family of sets satisfying the following conditions.
$\c{S}$ contains the empty set.
$\c{S}$ is closed under finite intersections.
The complement of any set in $\c{S}$ can be written as a finite union of disjoint sets in $\c{S}$.

Then any bounded countably additive non-negative function $\mu : \c{S} \-> \R$ satisfying $\mu(\emptyset) = 0$ extends uniquely to a measure defined on the smallest $\sigma$-algebra containing $\c{S}$.

Define the algebra of sets generated by $\c{S}$ to be

\[
a(\c{S}) = \cbr{\U_{i=1}^n S_i : S_i \in \c{S} \t{disjoint}}
\]

and note that the smallest $\sigma$-algebra generated by $\c{S}$ obviously coincides with the smallest $\sigma$-algebra generated by $a(\c{S})$.
We claim every function $\mu : \c{S} \-> \R$ satisfying the given assumptions extends uniquely to a function on $a(\c{S})$.
Every element of $a(\c{S})$ can be written as a finite union of disjoint sets—using this, define

\[
\mu^{(a)} : a(\c{S}) &\-> \R
\mu^{(a)}\del{\U_{i=1}^n S_i} &= \sum_{i=1}^n \mu(S_i)
\]

To ensure this is well-defined, we check that the definition is independent of the choice of which disjoint subsets to take the union of—suppose that $\U_{i=1}^n S_i = \U_{j=1}^m T_j$.
Then we have $S_i = \U_{j=1}^m S_i \^ T_j$ and similarly $T_j = \U_{i=1}^n S_i \^ T_j$: plugging these in to $\mu$ yields a double sum of disjoint sets and affirms well-definedness.
Next, note that $\mu^{(a)}(\emptyset) = \mu(\emptyset) = 0$.
To see that $\mu$ inherits countable additivity, take a sequence $S_1,S_2,.. \in a(\c{S})$ and note that

\[
\U_{n=1}^\infty S_n = \U_{n=1}^\infty \U_{m=1}^{p_n} T_{nm}
\]

where $T_{nm} \in \c{S}$ and $p_n \in \N$.
Plugging this into $\mu^{(a)}$ and applying countable additivity of $\mu$ gives the desired property.
We have thus obtained a uniquely defined countably additive function $\mu^{(a)} : a(\c{S}) \-> \R$ defined on an algebra of sets $a(\c{S})$ which extends $\mu$.
This function satisfies the assumptions of Carathéodory's Extension Theorem—see [Theorem A1.1]kallenberg06 or [Theorem 3.1]billingsley08—applying this result gives the claim, where we note that uniqueness follows since the range of $\mu^{(a)}$ is bounded.

Let $\pi : \c{A} \x \c{B} \-> \R$ be a countably additive function with $\pi(\emptyset) = 0$.
Then $\pi$ extends uniquely to a measure on the product $\sigma$-algebra.

By <Ref>, it suffices to show that $\c{A} \x \c{B}$ is a semi-algebra of sets.
The first required property is immediate, the second and third properties follow by $(A \x B) \^ (A' \x B') = (A \^ A') \x (B \^ B')$ and $(A \x B)^c = (A^c \x B) \u (A \x B^c) \u (A^c \x B^c)$.
The claim follows.

Let $\pi_{a,b}$ be a measure which admits the density $f_{a,b}$ with respect to a product measure $\lambda_{a,b}$.
Then $\pi_a = \pi_{a,b}(\.\x B)$ admits the density $f_a$ with respect to $\lambda_a = \lambda_{a,b}(\.\x B)$, and similarly in the other argument.

By assumption and Tonelli's Theorem, we have

\[
\pi_a(A_a) &= \pi_{a,b}(A_a \x B)
\\
&= \int_{A_a \x B} f_{a,b}(\alpha,\beta) \d\lambda_{a,b}(\alpha,\beta)
\\
&= \int_{A_a} \int_B f_{a,b}(\alpha,\beta) \d\lambda_b(\beta)\d\lambda_a(\alpha)
\\
&= \int_{A_a} f_a(\alpha) \d\lambda_a(\alpha)
\]

where the desired probability density is $f_a(\alpha) = \int_B f_{a,b}(\alpha,\beta) \d\lambda_b(\beta)$.

§ STATISTICAL DECISION-MAKING

We now use the Bayesian formalism to construct a probabilistic theory of decision-making.
To begin, we formalize the very general concept of an abstract agent making decisions in an environment in pursuit of some goal.
We use sutton18,bertsekas19 as our references.

§.§ Markov decision processes

A Markov decision process is a stochastic system consisting of a set of states, actions, and transitions between states, together with a rewards that vary depending on states and actions.
This is defined as follows.

A discrete-time Markov decision process is a $4$-tuple consisting of the following.
State space: a measurable space $S$.
Action space: a measurable space $A$.
Reward: a probability kernel $r : \c{B}(\R) \x S \x A \-> \R$.
Transition kernel: a probability kernel $p : \c{S} \x S \x A \-> \R$.

This is a very broad notion—however, a number of variations are also possible.
For instance, one can consider continuous-time, purely deterministic, and partially observed analogs—each of these involve their own subtleties and deserve study in their own right, but we do not pursue them here.

external/export next=false

[shift=(-1.5,2.125), xscale=-0.0085, yscale=0.0085, thick]
[preaction = draw,ultra thick] (300, -395) – (285, -355) .. controls (285, -355) and (269.142, -342.071) .. (255, -335) .. controls (245, -330) and (227.237, -322.04) .. (215, -320) .. controls (185, -315) and (175, -295) .. (175, -295) .. controls (175, -295) and (150, -295) .. (135, -285) .. controls (129.117, -281.078) and (100, -250) .. (100, -235) .. controls (100, -215) and (105, -205) .. (115, -195) .. controls (135.274, -174.726) and (144.617, -159.069) .. (186.751, -142.257) .. controls (230, -125) and (245, -125.009) .. (255, -125) .. controls (270.117, -124.987) and (289.478, -127.154) .. (305.786, -129.965) .. controls (335, -135) and (347.49, -144.511) .. (364.446, -158.01) .. controls (378.439, -169.15) and (389.863, -183.998) .. (397.853, -200) .. controls (405.048, -214.411) and (409.778, -230.664) .. (410, -246.769) .. controls (410.144, -257.217) and (410.264, -270.263) .. (402.959, -277.316) .. controls (395, -285) and (385, -290) .. (385, -290) .. controls (385, -290) and (377.69, -310.274) .. (375.386, -315.536) .. controls (371.018, -325.511) and (360.595, -330.936) .. (350.937, -334.628) .. controls (340.799, -338.503) and (318.451, -336.817) .. (318.451, -336.817) – (326.08, -380.928) – cycle;
(300, -400) – (300, -275);
(275, -350) – (275, -250) – (225, -250);
(250, -350) – (250, -275);
(225, -325) – (225, -300) – (200, -275) – (100, -275);
(325, -350) – (325, -250) – (300, -250);
(400, -300) – (350, -300) – (350, -275) – (375, -250) – (375, -200) – (400, -200);
(350, -250) – (350, -225) – (250, -225) – (250, -175) – (225, -150) – (150, -150);
(350, -200) – (300, -200) – (300, -175);
(400, -175) – (325, -175) – (300, -150) – (300, -125);
(200, -250) – (150, -250) – (125, -225) – (75, -225);
(175, -200) – (225, -200) – (225, -225) – (175, -225);
(225, -175) – (150, -175) – (150, -200) – (100, -200);
(275, -200) – (275, -125);
[fill=white] (200, -250) circle (7.5);
[fill=white] (225, -250) circle (7.5);
[fill=white] (250, -275) circle (7.5);
[fill=white] (300, -275) circle (7.5);
[fill=white] (300, -250) circle (7.5);
[fill=white] (350, -250) circle (7.5);
[fill=white] (350, -200) circle (7.5);
[fill=white] (175, -225) circle (7.5);
[fill=white] (175, -200) circle (7.5);
[fill=white] (225, -175) circle (7.5);
[fill=white] (300, -175) circle (7.5);
[fill=white] (275, -200) circle (7.5);

[shift=(2.5,-1), scale=0.0125, very thick, line cap=round]
(2, 2) – (66, 2);
(82, 2) – (162, 2);
(50, 18) – (114, 18);
(146, 18) – (162, 18);
(18, 34) – (50, 34);
(98, 34) – (114, 34);
(2, 50) – (34, 50);
(50, 50) – (66, 50);
(114, 50) – (146, 50);
(18, 66) – (34, 66);
(66, 66) – (114, 66);
(130, 66) – (146, 66);
(34, 82) – (66, 82);
(98, 82) – (130, 82);
(2, 98) – (18, 98);
(34, 98) – (50, 98);
(82, 98) – (98, 98);
(146, 98) – (162, 98);
(50, 114) – (66, 114);
(82, 114) – (98, 114);
(2, 130) – (18, 130);
(34, 130) – (50, 130);
(98, 130) – (114, 130);
(130, 130) – (146, 130);
(18, 146) – (34, 146);
(50, 146) – (66, 146);
(114, 146) – (130, 146);
(146, 146) – (162, 146);
(2, 162) – (82, 162);
(98, 162) – (162, 162);
(2, 2) – (2, 162);
(18, 18) – (18, 34);
(18, 66) – (18, 82);
(18, 98) – (18, 114);
(18, 130) – (18, 146);
(34, 2) – (34, 18);
(34, 50) – (34, 66);
(34, 82) – (34, 114);
(50, 18) – (50, 34);
(50, 50) – (50, 82);
(50, 114) – (50, 130);
(50, 146) – (50, 162);
(66, 34) – (66, 66);
(66, 82) – (66, 114);
(66, 130) – (66, 146);
(82, 18) – (82, 50);
(82, 66) – (82, 98);
(82, 114) – (82, 162);
(98, 34) – (98, 66);
(98, 98) – (98, 114);
(98, 130) – (98, 162);
(114, 18) – (114, 34);
(114, 82) – (114, 130);
(130, 2) – (130, 34);
(130, 50) – (130, 146);
(146, 34) – (146, 50);
(146, 82) – (146, 98);
(146, 114) – (146, 130);
(162, 2) – (162, 162);

[minimum size=85] at (-3.5,0) (a) ;
[minimum size=85] at (3.5,0) (e) ;
[anchor=north] at (a.south) Agent;
[anchor=north] at (e.south) Environment;
[-latex,very thick] (a) to[bend left=22.5] node[midway, above] action (e);
[-latex,very thick] (e) to[bend left=22.5] node[midway, below] reward, next state (a);

[Markov decision processes]Illustration of the feedback loop induced by a Markov decision process. Here, the agent chooses an action, which results in the environment changing to a new state. The agent observes the new state, as well as the reward given by the previous state and action. The agent's goal is to select actions to maximize long-term rewards.

The idea behind this definition is that, at every point in time, the agent observes the current state, chooses an action $a \in A$, and obtains another state $s' \~ p(s\given a)$.
<Ref> illustrates this.
Note that time can, and often will, be part of the state $s$.
The agent's goal is to choose actions so as to control the full trajectory of states in order to obtain maximum rewards.
The choice of action in every state is called a policy, and is formalized as follows.

Define the following.
A measurable function $\pi : S \-> A$ is called a deterministic policy.
A probability kernel $\pi : \c{A} \x S \-> \R$ is called a Markov policy.

Markov policies include deterministic policies as a special case, by taking the conditional distributions to be Dirac and re-interpreting the given expressions appropriately.
As with Markov decision processes, one can also consider more general classes of policies, but we do not do so here.

Different policies yield different state trajectories, and therefore different rewards.
Of these, some obtain more rewards than others.
We therefore introduce notions for distinguishing between policies according to the rewards they obtain.

Let $T \in \N$ be the time horizon.
A policy is called optimal if it maximizes the value function

\[
V^{(\pi)}(s_0) = \E \sum_{t=0}^T r_t
\]

where $r_t \given s_t,a_t \~ r(s_t,a_t)$, $a_t \given s_t \~ \pi(s_t)$, and $s_{t+1} \given s_t,a_t \~ p(s_t, a_t)$.

Finding an optimal policy therefore amounts to selecting the best possible actions to maximize expected total rewards.
Note that closely-related alternative notions of optimality, such as minimizing infinite discounted sums, or limits of averages, are also possible.
The most important distinction between different decision problems for finding optimal policies is given by what is assumed known.

If $p$ and $r$ are known, we say we have an optimal control problem.
Otherwise, we say we have a reinforcement learning problem.

These classes differ fundamentally from one another.
Optimal control problems can be viewed as a class of structured optimization problems, where the goal is to compute $\pi$ by evaluating $r$ and $p$ as necessary.
Here, one generally proceeds by proving that $V^{(\pi)}$ and $\pi$ satisfy certain recursive equations, and developing schemes for solving them directly.

Reinforcement learning problems are more complex.
Due to the rewards or dynamics being unknown, they cannot simply be maximized and their expectation must be learned.
This forces one to consider whether to take advantage of actions known to be good, or to try others in case they might be better—this is known as the explore-exploit tradeoff.
For such settings, we need an appropriate solution concept—to obtain one, define the following.

The regret of a policy $\pi$ is defined as

\[
R^{(\pi)}(s_0) = V^*(s_0) - V^{(\pi)}(s_0)
\]

provided that the value function $V^*$ with respect to an optimal policy exists.

Minimizing regret is equivalent to maximizing value, but when $p$ and $r$ are unknown doing so directly is impossible.
Instead, the goal is to find an algorithm—that is, a way of updating the policy based on observed data—so as to limit growth of regret.

In most settings, one can prove that every algorithm which does not know $p$ and $r$ necessarily incurs some level of regret.
This is done by exhibiting a randomized set of problems and rewards over which regret is lower-bounded in expectation for any algorithm.
In such a class, actions that perform well on one problem will necessarily perform badly on another problem.
Such arguments show that some degree of regret is inevitable.

On the other hand, some algorithms incur more regret than others.
The obviously-bad algorithm which learns nothing and chooses the exact same action over and over again incurs at most linear regret.
An algorithm is said to solve a decision problem if its asymptotic regret rate with respect to $T$ matches the respective regret lower bound in the given problem class.
Finding such algorithms is of key interest.

One way to construct algorithms for solving decision problems is to employ a model-based approach, which loosely speaking works as follows.

Learn the unknown transitions and/or rewards from observed data using a supervised learning approach.
Use the learned model(s) to find a policy satisfying some criteria.

The details of such approaches depend on the setting.
We distinguish between two key kinds of reinforcement learning problems.

If $|S| = 1$, we say we have a multi-armed bandit problem.
Otherwise, we say we have a general reinforcement learning problem.

Multi-armed bandits possess no variable state $S$, and thus only require one to learn the rewards to determine optimal actions.
In general reinforcement learning, actions can influence transitions between states—these problems require long-term planning, making them much more general, difficult, and important.
One can argue that as a mathematical theory, reinforcement learning is powerful enough to describe many aspects of human and animal intelligence, making it fundamentally interesting to study and develop.

We now restrict ourselves to the bandit setting, which is substantially easier to study and so can be understood more deeply.
Here, even when $A$ is a finite set and the rewards are Gaussian, the model-based approach consisting of (i) estimating rewards using empirical risk minimization and (ii) choosing the policy which maximizes rewards is known to be non-optimal.
This approach fails to explore, and can get stuck chasing inferior rewards.

One way to fix this problem is to replace empirical risk minimization with Bayesian learning, and adopt an appropriate decision rule for selecting actions.
Doing this yields approaches which can be shown optimal in many settings.
We therefore proceed to study multi-armed bandits in more detail.

§.§ Multi-armed bandits

The multi-armed bandit problem takes its name from a casino analogy.
In the 1950s, slot machines often had levers one could pull instead of buttons one could press, and were called one-armed bandits for their ability to empty gamblers' wallets.
A multi-armed bandit is a slot machine, which, for a fixed cost, allows one to pull an arm $x \in X$ and receive a random reward whose distribution depends on $x$.
The goal is to minimize  expected loss, or, equivalently, maximize total expected rewards.

Multi-armed bandits can be viewed as discrete-time Markov decision processes with a one-element state space, but this is not necessarily the most fruitful way to think about them.
We thus begin by introducing formalism and notation better suited to the given setting, which can be viewed as a special case of the notions considered previously.
We use slivkins19,lattimore20 as references.

Let $f : X \-> \R$ be a bounded above measurable function, let $\eps : \Omega \x X \-> \R$ be a jointly measurable stochastic process such that $\E(\eps(x)) = 0$ for all $x$, and let $y(\omega,x) = f(x) + \eps(\omega,x)$.
Define the following.
We say that the Markov decision process $(\{1\},X,y_*\P,\delta_{1})$, fully defined below, is a multi-armed bandit.
We say that a probability kernel  $\pi : \c{X} \x \bigoplus_{n=1}^\infty (X \x \R)^n \-> \R$ is a multi-armed bandit algorithm.

Here, $\delta_1$ is the Dirac measure centered at $1$ for all actions $x\in X$, which is the only possible conditional probability distribution over a one-element set, and $\bigoplus$ denotes the disjoint union of measurable spaces.

An algorithm, then, assigns every dataset of arbitrary size to a probability measure describing what arms should be picked with what probability.
Each dataset consists of $(x, y)$ pairs where $x$ are the locations chosen by the algorithm, and $y$ are the noisy observed values—recall $\omega$ is the randomness used by the noise.
Some algorithms maximize $f$ faster than others—we consider this next.

For a given multi-armed bandit, let $f(x^*)$ be the global maxima of $f$, which we assume to exist.
Define the regret of an algorithm $\pi$ at time $T$ to be

\[
R(\omega,T) = \sum_{t=1}^T f(x^*) - f(x_t(\omega))
\]

where $x_t \~ \pi(x_1,y_1,..,x_{t-1},y_{t-1})$, $y_t(\omega) = f(x_t) + \eps_t(\omega)$, and $\eps_t \~ \eps(\.,x_t)$.

[Multi-armed bandits]Here, we illustrate a simple three-armed bandit. Each of the three arms has its own reward distribution, shown on the left. These are unknown to the agent, which only sees the rewards obtained by actions it takes, shown in the center, and must use the obtained information to decide which arm to pick. Each arm's regret is given by its expected decrease in reward compared to the optimal arm, shown on the right.

Regret thus counts the total reward lost by virtue of not knowing the optimal arm in advance—this is illustrated in <Ref>.
Algorithms which achieve small regret therefore learn the optimal rewards effectively, and avoid getting stuck pulling bad arms.

Regret behavior in multi-armed bandit problems is strongly dependent on properties of the underlying function $f$, and in particular its domain $X$.
It is clear by considering for instance $X = \R$ that if $X$ is too large or too unstructured, no algorithm achieves better than linear regret.
We therefore begin studying the simplest non-trivial domain class.

We say that a multi-armed bandit defined over a finite set with cardinality $K=|X|$ is a $K$-armed bandit.
Moreover, if $y(\.,x)$ is a Bernoulli random variable for all $x$, we say it is a Bernoulli bandit.

For such bandits, regret can be decomposed on a per-arm basis in the manner given as follows.

Regret satisfies

\[
R(\omega,T) = \sum_{x=1}^K \Delta(x) n_T(\omega,x)
\]

where $n_T(\omega,x) = \sum_{t=1}^T \1_{x_t(\omega) = x}$ and $\Delta(x) = f(x^*) - f(x)$.

The claim follows directly from the definitions of $R$ and $n_T$ by grouping terms in the sum.

What kind of performance is possible on such a problem?
One can ask and answer this question as follows.

For any algorithm defined over a class of $K$-armed bandits there is an $f$ such that

\[
\E(R(\.,T)) \geq \Omega(\sqrt{KT})
\]

[Theorem 2.11]slivkins19.

A similar lower bound, but where the right-hand-side depends on an extra term controlling the difficult of the problem, was originally proved by lai85.
Such bounds are called instance-dependent: in that case, the rate obtained is logarithmic, owing to the presence of the extra term.
In our analysis, we focus on instance-independent bounds such as the one presented above.

This result tells us that regret is necessarily incurred as consequence of not knowing the expected rewards of each arm given by $f$.
The next step, then, is to ask: is there an algorithm which achieves this rate?
We first consider simply evaluating each arm once, and then making choices according to the empirical averages.

For a Bernoulli bandit, the algorithm which chooses actions by first trying each arm out once, and then selecting arms according to the maximum empirical average

\[
x_{t+1} &= \argmax_{x\in X} \mu_t(x)
\mu_t(x) &= \frac{\sum_{t=1}^T \1_{y_t = 1} \1_{x_t = x}}{\sum_{t=1}^T \1_{x_t = x}}
\]

of the random observed data achieves linear regret.

We first show that it suffices to exhibit a bandit along with an event $S$, which we call the stuck event, which occurs with constant probability and induces linear regret.
To see this, for an event $S$, write

\[
\E(R(\.,T)) &= \E(R(\.,T) \given S) \P(S) + \E(R(\.,T) \given S^c) \P(S^c)
\\
&\geq \E(R(\.,T) \given S) \P(S)
\]

Now, take $S$ to be the event that during the initial first arm pulls, the optimal arm yields zero reward, while some other arm yields a reward: clearly $\P(S) > 0$.
On the other hand, $\E(R(\.,T) \given S) = \c{O}(T)$, because the optimal arm conditional on $S$ has empirical mean zero, which is always smaller than some alternative, and will never be selected again.
The claim follows.

This algorithm fails to resolve the explore-exploit tradeoff, and gets stuck with suboptimal choices.
While the setup is obviously contrived—for instance, giving each arm some non-zero selection probability even if it previously gave bad rewards would seemingly break the lower bound—it also illustrates an important principle: the algorithm must explore.
The question, then, is how should it explore?
As an alternative, consider a simple uncertainty-based decision rule for selecting arms.

Upper confidence bound

[Upper confidence bound algorithm]The idea behind the upper confidence bound algorithm is to use the observed data to construct a set of error bars that reflect what has been learned about the mean rewards.
Here, we construct these using the posterior distribution under a Gaussian prior and likelihood.
The next arm is chosen as the maximum of the quantile-based error bars.
This process is repeated iteratively as additional data is obtained, with the quantile becoming more strict over time, ensuring the bounds hold increasingly often.

Define the Hoeffding upper confidence bound algorithm for Bernoulli bandits, which tries each arm once, and then, for a total of $T$ rounds, selects actions by

\[
x_{t+1} &= \argmax_{x\in X} f^+_t(x)
f^+_t(x) &= \mu_t(x) + c_t\sigma_t(x)
\]

where $\mu_t$ are the empirical means, $\sigma_t(x) = \sqrt{\frac{1}{n_t(x)}}$, and $c_t = \sqrt{2\ln(T)}$.

Many algorithms within the upper confidence bound class have been proposed <cit.>: we focus on the variant analyzed by auer02, but illustrate an alternative in <Ref>.
We study the Hoeffding-based algorithm's regret, and present a simplified but less sharp analysis compared to the one given in that work.
It turns out the simple modification of adding a set of error bars, with a threshold growing in time and decreasing with data size, is enough to result in favorable regret behavior.

The Hoeffding upper confidence bound algorithm achieves an expected regret of

\[
\E(R(\.,T)) \leq \tl{\c{O}}(\sqrt{KT})
\]

uniformly for all $f$, where $\tl{\c{O}}$ denotes asymptotics up to logarithmic factors.

We adapt the argument presented by slivkins19.
First, let $B$ denote the event that the bound holds, namely $f^-_t(\.,x) \leq f(x) \leq f^+_t(\.,x)$ for all $x\in X$ and all $t \leq T$, where $f^-_t(\.,x) = \mu_t(\.,x) - c_t \sigma_t(\.,x)$, where all stochasticity is defined with respect to the data.
Under this event, the upper confidence bound is a true bound on the rewards.
Thus, we have

\[
\E(R(\.,T)) &= \E(R(\.,T) \given B) \P(B) + \E(R(\.,T) \given B^c) \P(B^c)
\\
&\leq \E(R(\.,T) \given B) + T \P(B^c)
\]

since $\P(B) \leq 1$, and a regret of $T$ is the maximum possible value as there are $T$ rounds and rewards are in $\{0,1\}$.
Our strategy will be to bound $\E(R(\.,T) \given B)$, while choosing the scaling factor $c_t$ in the width of the upper confidence bound to ensure that $\P(B^c)$ decays fast enough that the latter term is negligible.
We begin with the latter.
By the union bound, we have

\[
\P(B^c) &\leq \sum_{t=1}^T \sum_{x\in X} \P(p_x \leq f^-_t(\.,x)) + \P(p_x \geq f^+_t(\.,x))
\\
&\leq 2KT \max_{\substack{x\in X\\1 \leq t \leq T}} \P(p_x \geq f^+_t(\.,x))
\]

where $p_x$ is the Bernoulli parameter for arm $x$, and we have used symmetry of the error bars to combine terms.
For each $n > 0$, if $y$ is a sum of independent Bernoulli random variables, the probability we want to bound is

\[
\P\del{p_x \geq \frac{y(\.)}{n} + \sqrt{\frac{2\ln(T)}{n}}} \leq \exp\del{-2 \del{\sqrt{\frac{2\ln(T)}{n}}}^2 n} = \frac{1}{T^4}
\]

using Hoeffding's inequality.
Since holds uniformly in $n$, it follows that

\[
\P(p_x \geq f^+_t(\.,x)) = \P\del{p_x \geq \frac{\sum_{x_t = x} y_t(\.)}{n_t(\.,x)} + \sqrt{\frac{2\ln(T)}{n_t(\.,x)}}} \leq \frac{1}{T^4}
\]

We therefore conclude that

\[
T\P(B^c) \leq \frac{2KT}{T^4} = \c{O}(1)
\]

which shows that the upper confidence bounds are exceeded sufficiently rarely that this possibility incurs no regret in the asymptotic limit, and completes this part of the argument.
We now proceed to analyze the $\E(R(\.,T) \given B)$ term—assume henceforth that $B$ holds.
We have that

\[
\Delta(x_t(\omega)) &= f(x^*) - f(x_t(\omega))
\\
&\leq f^+_t(\omega,x^*) - f^-_t(\omega,x_t(\omega))
\\
&\leq f^+_t(\omega,x_t(\omega)) - f^-_t(\omega,x_t(\omega))
\\
&= 2\sigma_t(\omega,x_t(\omega))
\]

using the definition of the event $B$.
Now, if we consider the rescaled error bar width $c_t\sigma_t(\omega,x_t(\omega))$ of the arm chosen at time $t$, this is

\[
\Delta(x_t(\omega)) \leq \c{O}(c_t\sigma_t(\omega,x_t(\omega))) = \tl{\c{O}}\del{\frac{1}{\sqrt{n_t}}}
\]

By <Ref>—which we note still holds conditional on $B$ via the exact same argument—we have

\[
\E(R(\.,T)\given B) &= \sum_{x=1}^K \Delta(x) n_T(\.,x)
\\
&\leq \sum_{x=1}^K \tl{\c{O}}\del{\sqrt{n_T(x)}}
\\
&\leq \tl{\c{O}}\del{\sqrt{K} \sqrt{\sum_{x\in X} n_T(x)}}
\\
&= \tl{\c{O}}(\sqrt{KT})
\]

using that by definition, every $x$ for which $n_T(\omega, x) > 0$ is equal to $x_t(\omega)$ for some $t$, enabling us to apply the preceding inequality.
The claim follows.

This shows that the proposed algorithm, which uncertainty built via the error bars, effectively balances exploration and exploitation in this setting.
The above behavior is not unique to the given form of the upper confidence bound, nor to the upper confidence bound rule itself.
For example, we can consider a variation of the above algorithm, where the width of the error bars is constructed via a Bayesian model.

Define a Bayesian model via the likelihood $\gamma(x) \~[Ber](\mu(x))$ and prior $\mu(x) \~[Beta](a,b)$.
Define the beta–Bernoulli upper confidence bound algorithm which selects actions by maximizing the function

\[
x_{t+1} &= \argmax_{x\in X} f^+_t(x)
f^+_t(x) &= \mu_t(x) + c_t \sigma_t(x)
\]

where $(\mu_t, \sigma_t)$ are the mean and standard deviation of the posterior distribution $\mu\given\gamma(x_1) = y_1, .., \gamma(x_t) = y_t$, and $c_t$ is defined as previously.

One can show this algorithm, illustrated in <Ref>, achieves the same regret as its Hoeffding-based analogue: the proof is similar, but instead uses the empirical Bernstein inequality, and is more messy owing to the complicated form of the upper confidence bound.
Even more generally, one can select arms by optimizing a function $\alpha : X \-> \R$, called an acquisition function, built from a posterior distribution, confidence set, or other appropriate construction.
Many different acquisition functions have been proposed.

We have made no attempt to optimize the bound, and significant improvements are possible, particularly when considering variations of the upper confidence bound algorithm—see lattimore20 for state-of-the-art results.
Of these, algorithms built using confidence sets rather than posterior distributions, such as the variant presented previously, are particularly important in the $K$-armed bandit setting.
We focus on the Bayesian view because it generalizes well to more complex settings.

We now proceed to explore a more general setting which will enable us to employ the ideas developed so far to develop efficient black-box optimization algorithms.
This will enable us to use Bayesian methods to solve a broad class of decision problems of practical interest.

§.§ Bayesian optimization

We now describe the formalism of Bayesian optimization for using ideas built on Bayesian learning and multi-armed bandits to design global optimization algorithms.
We use frazier18 as our main reference.
Our goal now is to minimize a black-box function

\[
f : X \-> \R
\]

which is assumed continuous, and defined on a compact set $X \subseteq \R^d$.
Such a function is automatically bounded.
Our goal is to minimize $f$ with as few evaluations as possible.
To measure performance, we again introduce a notion of regret.
One choice is to define

\[
R(T) = \sum_{t=1}^T f(x_t) - f(x^*)
\]

where we have used the opposite sign convention compared to bandits and reinforcement learning, because our goal is to minimize $f$ rather than maximizing rewards.
Note that unlike before, for a deterministic algorithm this is now a purely deterministic quantity, and not a random variable.

As before, we can approach this problem by building a Bayesian model.
For an arbitrary sequence of points $x_1,..,x_t$, define the likelihood

\[
y_t(\omega) &= f(x_t) + \eps_t(\omega)
\eps_t &\~[N](0,\sigma^2)
\]

Here, the space of observed values is $\R$, and our quantity of interest is the actual function $f$, which we at least implicitly view as an element of a space of functions, such as for instance the Banach space of continuous functions $C^0(X;\R)$.

$t = 1$

$t = 2$

$t = 3$

$t = 4$

[Bayesian optimization using Thompson sampling]Bayesian optimization using Thompson sampling.
Here, we sample a random function from the posterior distribution over possible functions, and choose the next evaluation point to be the minima of the sampled random function.
As additional data is obtained, the posterior distribution concentrates around the true function.
Unlike the upper confidence bound acquisition function considered previously, Thompson sampling makes use of a full probability distribution, rather than just a set of error bars.

It is now clear why we bothered with setting up a dense, abstract formalism for Bayesian learning in general measure spaces: this formalism is rich and powerful enough to enable us to properly define priors on spaces like this—we explore these in the sequel.
Suppose for the moment that this is possible.
Then, we can calculate the posterior distribution

\[
f \given y_1,..,y_t
\]

which is now a probability measure supported on $C^0(X;\R)$.
To determine the next point to query, we introduce and optimize an acquisition function.
Typical acquisition functions include the upper confidence bound <cit.> acquisition function considered previously, as well as Thompson sampling <cit.>, first proposed half a century before Bayesian optimization was otherwise popularized.
This is defined as a random acquisition function given by

\[
x_{t+1}(\omega) &= \argmin_{x\in X} \alpha_t(\omega,x)
\alpha_t&\~ f \given y_1,..,y_t
\]

We show an example of Bayesian optimization using Thompson sampling in <Ref>.
Another typical choice is given by expected improvement <cit.>, which is defined as

\[
x_{t+1} = \argmax_{x\in X} \E \max(0,f_t(\.,x) - f(x^*_t))
\]

where $f_t = f \given y_1,..,y_t$ is the posterior, and $x^*_t = \argmin_{t \in \{1,..,t\}} f(x_t)$ is location with the smallest value observed so far.
Many other effective acquisition functions, including probability of improvement <cit.>, predictive entropy search <cit.> and information directed sampling <cit.>, are also possible.

Regret analysis for all of these choices can be performed, and will in general depend on the detailed properties of the Bayesian model, acquisition function, and unknown function $f$ <cit.>.
In particular, regularity and smoothness properties may play a role, as well as structure present in the domain $X$.
In settings where the function $f$ is unknown, it's also possible to analyze

\[
BR(T) = \E \sum_{t=1}^T f(\omega,x_t) - f(\omega,x^*(\omega))
\]

where regret is now considered in expectation with respect to the prior—this is called the Bayesian regret.
Some acquisition functions, such as Thompson sampling, admit particularly simple analyses in this setting <cit.>.

This concludes our development and showcasing of decision-making algorithms.
We now proceed to develop the final piece of the puzzle not yet studied in detail: how to place priors on function spaces, in order to build Bayesian models for settings such as Bayesian optimization.

§ GAUSSIAN PROCESSES

In the preceding section, considerations arising from Bayesian decision-making algorithms motivated us to find a way to place priors on spaces of functions $f : X \-> \R$.
Gaussian processes are a broad class of random functions capable of this.
In <Ref>, to perform Bayesian optimization using Thompson sampling, we used a Gaussian process model to define the posterior distribution over random functions, as needed by the decision rule.

To develop Gaussian processes, we will need to start with simpler notions and gradually work towards increasing levels of generality.
Gaussian processes are defined by the key property that no matter what angle one views them from, they yield Gaussian marginal distributions.
We therefore begin by studying properties of Gaussian random variables and Gaussian random vectors, which will be used as the basic building blocks for the general case.

The notion of a Gaussian process as a random function will turn out to be too restrictive for all of our settings of interest: it is not possible to view certain random variables, which do deserve to be called Gaussian processes, in this way.
The obstructions can be geometric or analytic in nature.
For example, a vector field on a manifold is not a vector-valued continuous function, it is a section of a vector bundle: what, then, should the term Gaussian actually mean?
Difficulties also occur when considering white noise processes.

To handle these issues, we work with the notion of a Gaussian process $f : \Omega \-> V$ where $V$ is a real vector space equipped with additional structure arising from the setting at hand.
Gaussian process are fundamentally linear objects which reflect this structure.
The simplest possible settings one can study arise when $V$ is smallest.

Choosing $V = \{0\}$ to be the trivial vector space, there is exactly one $V$-valued random variable, whose distribution is the Dirac measure centered at $0$, which can trivially be called Gaussian.
This random variable is not very interesting, so we do not consider it further.
Choosing $V = \R$ to coincide with the underlying scalar field yields the setting of Gaussian random variables.
This is the next-simplest setting and the first one we explore in detail.
Choosing $V = \R^d$ yields the setting of Gaussian random vectors, whose components are multivariate Gaussian—or, equivalently, whose dot products are all scalar-valued Gaussian.
Choosing $V$ to be a space of functions $f : X \-> \R$ yields the setting of Gaussian processes, whose finite-dimensional marginals are multivariate Gaussian.
This setting is well-studied when $X$ is itself a Euclidean space: in <Ref>, we will examine cases where $X$ instead possesses various kinds of geometric structure.
Finally, choosing $V$ to be a possibly infinite-dimensional topological vector space yields the most general setting we examine.
This level of generality will be important for two reasons: (i) to provide a formalism for studying stochastic partial differential equations, and (ii) to develop a coordinate-free notion of Gaussian random vectors that will be useful in differential-geometric settings.

In what follows, our goal will be to lay the groundwork for subsequent development.
Thus, we will not discuss Bayesian learning with Gaussian processes, which will instead be presented in <Ref>.
We use lifshits12 as our standard reference.
We first study the scalar case, before generalizing to other settings.

§.§ Review of vector spaces

Since Gaussian processes are closely connected to vector spaces and related mathematical structures, here we review ideas from functional analysis that are useful in understanding them.
We use lang12 as our reference on these topics.

Topological vector space
A set equipped with a notion of addition and multiplication by scalars taking values in an underlying field is called a vector space.
We work with vector spaces where the underlying field is the real numbers.
A topological vector space is a vector space equipped with a topology under which the vector operations are continuous—the topology gives rise to notions of convergence, continuity, compactness, and many others.
Many interesting spaces of functions and generalizations thereof are topological vector spaces.

A Banach space is a topological vector space equipped with a norm, which assigns a notion of size to all vectors, and induces a notion of distance given by a topologically complete metric, which in turn induces the space's topology.
A Hilbert space, similarly, is a topological vector space equipped with an inner product which induces a complete norm.
Many, but not all, useful spaces of functions are Banach or Hilbert spaces.
Similarly, many, but not all, useful notions of convergence come from norms or inner products.

Every vector space admits a Hamel basis, which allows arbitrary vectors to be decomposed into linear combinations of linearly independent vectors in the basis.
Bases are generally highly non-unique.
The dimension of a vector space is the cardinality of such a basis, which does not depend on which basis is chosen.
Hamel bases of infinite-dimensional vector spaces are usually poorly-behaved, and alternative notions, such as that of a Schauder basis of a Banach space, are often used instead.

Given two vector spaces $V$ and $W$, we can form the direct sum vector space $V \oplus W$ by taking the Cartesian product of both spaces, and defining vector operations for $V \oplus W$ componentwise using the vector operations on $V$ and $W$.
This yields a vector space whose dimension is the sum of the dimensions of $V$ and $W$.
If $V$ and $W$ are topological vector spaces, equipping $V \oplus W$ with the product topology makes it into a topological vector space.
In this work, we only consider finite direct sums.

Linear operator
Maps between vector spaces are called linear operators if they preserve the vector operations.
By linearity, continuous maps between Banach spaces are automatically bounded with respect to their inputs, and vice versa, so in this context we refer to boundedness and continuity interchangeably.
Given two respective bases, linear maps between finite-dimensional vector spaces can be viewed as matrices via their action on a vector's basis coefficients.
In doing so, composition of linear maps becomes matrix multiplication.

Dual space
For a topological vector space $V$, the space $V^*$ consisting of continuous linear maps $\phi : V \-> \R$ into the underlying scalar field is called its dual space.
The dual space can be made into a vector space: for many topological vector spaces, including Banach and Hilbert spaces, it can canonically be assigned a topology making it into a topological vector space.
We generally call elements of $V^*$ linear functionals.

Space of bounded linear operators
The spaces of linear operators between two vector spaces can also be given the structure of a vector space.
Using this, if $V$ and $W$ are Banach spaces, one can define the space of $L(V;W)$ of bounded linear operators.
This space can be made into a Banach space via the operator norm, defined as the largest possible size of a unit-norm vector in $V$ after it is mapped into $W$ by the given operator.
The adjoint $\c{A}^* : H \-> G$ of a bounded linear operator $\c{A} : G \-> H$ between Hilbert spaces is defined by $\innerprod{\c{A}g}{h}_H = \innerprod{g}{\c{A}^* h}_G$.

Space of continuous functions
We now describe some particularly important function spaces.
The space of continuous functions is denoted $C^0(X;\R)$: if $X$ is a compact set, this space can be equipped with the supremum norm to obtain a Banach space.
This is not the only useful topology on this space: in many settings, for instance, the topology induced by convergence of all bounded linear functionals is also used.
The space of infinitely differentiable functions is denoted by $C^\infty(X;\R)$: this space, in contrast, is generally not a Banach space.

For a measure space $X$, the Lebesgue space $L^p(X;\R)$, defined for $1 \leq p < \infty$, is the Banach space of equivalence classes of measurable functions $f : X \-> \R$, identified up to almost sure equality, whose norm, given by integration of absolute $p$th powers, is finite.
Of these spaces, only $L^2(X;\R)$ is a Hilbert space.
The space $L^\infty(X;\R)$ is the Banach space of essentially bounded equivalence classes of functions equipped with the essential supremum norm.

Space of distributions
Many useful topological vector spaces are not spaces of functions.
The space of distributions $D'(X)$ is defined as the dual of the space $D(X)$ of infinitely differentiable functions with compact support, equipped with the topology of uniform convergence of the function and all derivatives on compact sets.
Locally integrable functions embed into $D'(X)$, so we can view it as a space of generalized functions, containing both functions and other, less regular elements such as the Dirac delta function, which is not a classical function.

§.§ Gaussian random variables

A Gaussian random variable is a map which takes in an abstract random number, and returns a real scalar.
The basic object from which other Gaussians will be constructed is the standard scalar Gaussian, defined as follows.

A random variable $z : \Omega\->\R$ is called standard Gaussian if it admits the Lebesgue density

\[
f(z) = \frac{1}{\sqrt{2\pi}} \exp\del{-\frac{z^2}{2}}
\]

From this, we define general scalar Gaussians.

A random variable $y : \Omega\->\R$ is called Gaussian if there are scalars $\mu, \sigma\in\R$, and a standard Gaussian $z$, such that

\[
y = \sigma z + \mu
\]

Note that we do not require $\sigma \geq 0$: hence, every Gaussian random variable is determined uniquely in
distribution by the pair $(\mu,\sigma^2)$, which we call its mean and variance, respectively.
We write $y \~[N](\mu,\sigma^2)$.
True to these parameter names, we have

\[
\E(y) &= \mu
\E\del{(y - \mu)^2} &= \sigma^2
\]

A Gaussian random variable is called centered if $\mu = 0$.
For a given variance $\sigma^2$ and standard Gaussian $z$, the expressions $\sigma z$ and $-\sigma z$ define two different centered Gaussians with the same distribution.
At this stage, pointing these distinctions may appear needlessly pedantic: they will become more pronounced and important once we consider more general objects.
The density of a Gaussian random variable, if it exists, takes on a form analogous to that of a standard Gaussian, namely

\[
f(y) = \frac{1}{\sqrt{2\pi}\sigma} \exp\del{-\frac{(y-\mu)^2}{2\sigma^2}}
\]

This density is visualized in <Ref>.
Note that the density will not exist if $\sigma^2 = 0$: the distributions of such Gaussians are Dirac measures centered at $\mu$.
By examining this density, one sees that Gaussian random variables respect the additive and multiplicative structures of the reals.

[Gaussian distributions]Probability density functions for a set of Gaussian random variables, each with different mean and standard deviation parameters. All of these can be obtained by shifting and rescaling any of the others via affine transformations.

Let $y \~[N](\mu,\sigma^2)$.
Then for $a,b \in \R$ we have that

\[
a y + b \~[N](a\mu+b, a^2\sigma^2)
\]

Immediate by definition.

This compatibility with linear structure will be true at all levels of generality we consider.
We now lift this definition to construct multivariate analogs.

§.§ Gaussian random vectors

A multivariate Gaussian random vector is a random variable taking values in $\R^d$.
We write vectors defined in $\R^d$ in bold italics to emphasize this distinction, and similarly distinguish matrices from linear maps by writing the former in bold upface letters.
As before, we begin by defining a standard Gaussian.

A random variable $\v{z} : \Omega\->\R^d$ is called standard multivariate Gaussian if its distribution is the product measure of the distributions of $d$ standard Gaussians.

Once the notion of a standard Gaussian is available, we can again define multivariate Gaussians as transformations of standard Gaussians.

A random variable $\v{y} : \Omega\->\R^d$ is called multivariate Gaussian if there is a vector $\v\mu\in\R^d$, matrix $\m{L}\in\R^{d\x d}$, and standard multivariate Gaussian $\v{z}$, such that

\[
\v{y} = \m{L}\v{z} + \v\mu
\]

$\rho = 0.9$

$\rho = -0.6$

[Multivariate Gaussian densities]Two multivariate Gaussian densities, in dimension two, with unit variances and different correlation coefficients $\rho$. As $|\rho| \to 1$, the individual components of the multivariate Gaussian become more and more dependent.

Every multivariate Gaussian is determined by its mean vector $\v\mu$ and positive semi-definite covariance matrix $\m\Sigma = \m{L}\m{L}^T$, and, as before, is called centered if $\v\mu = \v{0}$.
We write $\v{y}\~[N](\v\mu,\m\Sigma)$.
We can obtain a centered multivariate Gaussian with a given distribution by calculating a matrix square root of $\m\Sigma$, multiplying it with a standard Gaussian.
The mean and covariance are

\[
\E(\v{y}) &= \v\mu
\Cov(\v{y}) &= \E\del{(\v{y}-\v\mu)(\v{y}-\v\mu)^T} = \m\Sigma
\]

which mirror the previous situation.
If the determinant $|\m\Sigma|$ is non-zero, the density, displayed in two different ways in <Ref>, is

\[
f(\v{y}) = \frac{1}{\sqrt{(2\pi)^d|\m\Sigma|}} \exp\del{-\frac{1}{2}(\v{y}-\v\mu)^T\m\Sigma^{-1}(\v{y}-\v\mu)}
\]

which now might not exist even if the distribution of $\v{y}$ is not Dirac.
In such cases, one can see that at least some eigenvalues of $\m\Sigma$ must be zero, and so Gaussians which do not admit densities must, when viewed in an appropriate basis, be products of Dirac measures with Gaussians which do admit densities.
Already in the multivariate case, then, we see the technical power of densities weakening: this will become more pronounced as we consider more general settings.
Affine maps, however, behave as before.

Let $\v{y}\~[N](\v\mu,\m\Sigma)$. Then for $\m{A}\in\R^{d\x d}$ and $\v{b}\in\R^d$ we have

\[
\m{A} \v{y} + \v{b} \~[N](\m{A}\v\mu + \v{b}, \m{A}\m\Sigma\m{A}^T)
\]

Immediate by definition.

$\rho = 0.9$

$\rho = -0.6$

[Multivariate Gaussian quantile ellipsoids]Quantile ellipsoids for two multivariate Gaussian densities in dimension two with unit variances, different correlation coefficients $\rho$, and quantile levels equal to $0.8$, $0.95$, and $0.99$.

We now pause and reflect.
First, we distinguish $\R^d$ from a generic $d$-dimensional vector space: the former comes with a product structure $\R^d = \R \x .. \x \R$, including projection maps onto each coordinate, which in turn induce a canonical choice of inner product given by the Euclidean dot product.
A generic finite-dimensional vector space lacks this structure: it admits many different inner products, and provides for no canonical choice.

With this in mind, we observe that we have not actually used the product structure of $\R^d$: suppose that $V$ is a finite-dimensional inner product space, and define

\[
y(\omega)  &= \sum_{i=1}^d z_i(\omega)  e_i
z_i \~[N](0,1)
\]

where $e_i$ is any orthonormal basis.
By noting that the choice of basis $e_i$ induces a Borel isomorphism $V \isom \R^d$, it is easy to see that this definition is basis-independent, since linear maps associated with changes of orthonormal bases are represented by orthogonal matrices.
This expression therefore defines a standard Gaussian with respect to the given inner product.

Subtle difference such as the ones considered become increasingly important in the sequel.
Therefore, we introduce an alternative definition to help build intuition for later.

A random variable $\v{y} : \Omega \-> \R^d$ is called multivariate Gaussian if, for any vector $\v\phi \in \R^d$, the Euclidean dot product

\[
\v\phi\.\v{y} : \Omega \-> \R
\]

is univariate Gaussian.

This definition turns out to be equivalent to the original one.

The notions of multivariate Gaussians in the sense of transformations and in the sense of duality coincide.

Since the dot product is a linear map, it is clear that multivariate Gaussians in the sense of transformations are also Gaussian in the sense of duality.
To see the other direction, consider unit vectors $\v\phi$ where all coordinates except one are zero.
By using dot products with such vectors to reassemble the mean vector and covariance matrix, one sees that the claim follows from eigenvalue factorization of positive semi-definite matrices.

This definition allows one to begin imagining what a substantially more general notion of Gaussianity might look like: dot products become linear functionals, and covariance matrices become bilinear forms.
The preceding argument even suggests that such random vectors can be studied using spectral theory.
Of course, in infinite-dimensional settings, topological and analytic considerations come into play, making theory more difficult.
We develop these ideas subsequently, but first consider a simpler setting.

§.§ Gaussian random functions

We now consider Gaussian random functions, which are the first notion of a Gaussian random variable that is generally called a Gaussian process.
Here, we will adopt a bottom-up view which, from a technical perspective, departs slightly from the notions introduced so far.

Recall that for a set $X$, an $\R$-valued stochastic process is a map $f : \Omega \x X \-> \R$ measurable in its first argument.
If we have a set of points $x_1,..,x_n \in X$, we can plug them into $f$ to obtain a map $(f(\.,x_1),..,f(\.,x_n)) : \Omega \-> \R^n$, which, by virtue of being a product of measurable maps, is a random variable.
We call its distribution a finite-dimensional marginal distribution.
Using this notion, we define Gaussian processes.

Squared exponential

[Gaussian processes]Random functions generated from Gaussian processes with four different covariance kernels, which differ in particular according to their regularity, ranging from nowhere-differentiable in the exponential case to infinitely-differentiable in the squared exponential case.

Let $X$ be a set.
A random process $f: \Omega \x X \-> \R$ is called a Gaussian process if, for any finite set of points $x_1,..,x_n \in X$, the random variable $(f(\.,x_1),..,f(\.,x_n)) : \Omega \-> \R^n$ is multivariate Gaussian.

Immediately upon writing this definition, one is left to wonder: does it actually make sense?
In particular, do Gaussian processes in the sense given here exist?

Under general conditions on the finite-dimensional marginals, Kolmogorov's Consistency Theorem states that there exists a unique probability measure on a suitable $\sigma$-algebra, called the cylindrical $\sigma$-algebra and defined below, whose finite-dimensional projections coincide with the distributions of the random variables written above.
Thus, Gaussian processes exist so long as a family of multivariate Gaussians satisfying the conditions of Kolmogorov's Consistency Theorem can be found.
We visualize such processes in <Ref>.

The same results also allow us to reinterpret Gaussian processes as function-valued random variables, which, in many ways, are a more natural point of view to take.
Define $\R^X = \{f : X \-> \R\}$ and equip it with the cylindrical $\sigma$-algebra, namely the smallest $\sigma$-algebra containing all sets of the form $\{f \in \R^X : f(x_1),..,f(x_n) \in A\}$, where $A$ is a measurable set on $\R^n$, for all $n$, to make it into a measurable space.
Let $V \subseteq \R^X$, and equip it with the subset $\sigma$-algebra.
We can then reinterpret Gaussian processes as follows.

Let $X$ be a set.
A random variable $f: \Omega \-> V \subseteq \R^X$ is called a Gaussian process if, for any finite set of points $x_1,..,x_n \in X$, the random variable $(f(\.)(x_1),..,f(\.)(x_n)) : \Omega \-> \R^n$ is multivariate Gaussian.

It is clear that every Gaussian process in the random process sense induces a Gaussian process in the random variable sense, and vice versa.
Thus, these two notions are simply two different ways of viewing the same object.

We can then ask, what properties of multivariate Gaussians are still true in this setting?
Since a Gaussian process is uniquely determined by its finite-dimensional marginals, we need to find functions that generate a family of mean vectors and covariance matrices which are positive semi-definite.
The former is straightforward: every function $\mu : X \-> \R$ can be evaluated at a finite set of points $x_1,..,x_n$ to obtain a mean vector $\v\mu = \mu(x_1,..,x_n)$.
The latter requires only slightly more thinking.

A symmetric function $k : X \x X \-> \R$ is called a positive semi-definite kernel if, for any finite set of points $x_1,..,x_n\in X$, the kernel matrix

\[
\begin{bmatrix}
k(x_1,x_1) & \dots &k(x_1,x_n)
\\
\vdots & \ddots & \vdots
\\
k(x_n,x_1) & \dots & k(x_n,x_n)
\end{bmatrix}
\]

is positive semi-definite.

We can recover such a kernel from a given Gaussian process.

Define the covariance kernel of a Gaussian process to be

\[
k(x,x') = \Cov(f(\.,x),f(\.,x'))
\]

It is then clear that the pair $(\mu,k)$ uniquely define a Gaussian process, and so we write $f\~[GP](\mu,k)$.
Defining a centered Gaussian process then amounts to defining a positive semi-definite kernel.

In the Euclidean case, this can be done straightforwardly by noting that (i) the linear kernel $k(x,x') = \innerprod{x}{x'}$ is positive semi-definite by non-degeneracy of the inner product, (ii) that sums, powers, and limits of positive semi-definite kernels are positive semi-definite.
By this technique, we see that the widely-used exponential and squared exponential kernels are positive semi-definite.

For many spaces of functions $V \subseteq \R^X$, we can define addition and scalar multiplication, thereby making $V$ into a vector space.
If we do so, then, as before, affine maps preserve Gaussianity.

Let $f\~[GP](\mu, k)$, and equip $V$ and $W \subseteq \R^X$ with the topology of pointwise convergence.
Suppose that $f$ can be written as an infinite sum of deterministic basis functions as $f(\omega,x) = \sum_{n=0}^\infty w_i(\omega) f_i(x)$ where $w_i$ are independent Gaussian random variables.
If $\c{A} : V \-> W$ is a continuous linear map, and $b \in W$, then $\c{A} f + b$ is a Gaussian process.

It is clear by considering a set of finite-dimensional marginal distributions that adding $b$ preserves Gaussianity.

\[
(\c{A}f)(x) = \c{A} \del{\sum_{i=0}^\infty w_i(\omega) f_i(\.)}(x) = \sum_{i=0}^\infty w_i(\omega) (\c{A}f_i)(x)
\]

gives the statement for all finite-dimensional marginals, since sums and limits of Gaussian random variables are Gaussian.
The claim follows.

This is a substantially weaker assertion than previously: though we can conclude that affine maps of Gaussian processes are Gaussian, for a generic affine map we cannot immediately determine the form of the resulting kernel.
Moreover, the conditions are unnatural: by supposing them, we have more-or-less assumed the resulting kernel to exist.
The problem here is that the notion of a kernel is too rigid to permit a broad statement—an analogous property will hold if this is replaced with a different notion of covariance.

The situation for the other properties considered previously is even worse.
In particular, it is clear that, as an infinite-dimensional object, $f$ does not admit a probability density analogous to the ones considered previously, because an infinite-dimensional Lebesgue measure, in the sense of a locally finite translation invariant measure, does not exist.
It is also not clear what a standard Gaussian process, nor what the analog of a matrix square root of a positive semi-definite kernel, might be.

The loss of these technical tools has consequences on what can be said about Gaussian processes.
In <Ref>, we will study Gaussian processes whose domain $X$ is a Riemannian manifold.
In that setting, one cannot begin by proving positive semi-definiteness of linear kernels, because there is simply no useful analog of this concept.
Defining positive semi-definite kernels there turns out to be non-trivial, and the kernels we study do not admit simple expressions like they do in the Euclidean case.

We therefore proceed to adopt a function-analytic perspective which is more technical, but significantly more powerful, and will allow us to recover the previous set of tools in a much more pure and abstract form.

§.§ Gaussian processes in the sense of duality

In the preceding section, we established a notion of a Gaussian process, generalizing the notion of Gaussianity to random functions.
In developing this viewpoint, some of the appeal of Gaussianity was lost: in particular, unlike in the finite-dimensional setting, no simple covariance kernel transformation rule was available.
Here, we consider an alternative function-analytic view which avoids these limitations at cost of a higher degree of abstraction.

Let $V$ and $W$ be a pair of measurable real vector spaces equipped with a jointly measurable non-degenerate bilinear form $\dualprod{\.}{\.} : W \x V \-> \R$.
We say that a random variable $f : \Omega \-> V$ is a Gaussian process in the sense of duality if, for any $\phi \in W$, the scalar-valued random variable $\dualprod{\phi}{f}$ is Gaussian.

This is an exceedingly broad definition, which is more useful as an organizing framework connecting different concepts than as a technical tool in its own right.
At first glance, it is not clear what, if anything, this notion has to do with the Gaussian processes we have considered previously.
To better understand this, an illustrative example is in order.

Let $f \~[GP](\mu,k)$ be a Gaussian process defined on $[0,1]$ whose sample paths are almost surely continuous, and note by compactness that they are automatically bounded.
Our Gaussian process can therefore be viewed as a random variable $f : \Omega \-> C^0([0,1];\R)$.
We endow the latter space with the supremum norm, making it into a Banach space.
If we take another function $\phi \in C^\infty([0,1],\R)$, called the test function, we can define the pairing

\[
\dualprod{\phi}{f} = \int_0^1 \phi(x)f(x) \d x
\]

which by boundedness is almost surely finite.
Now, note that since sums of Gaussians are Gaussian, Riemann sums of Gaussian processes are Gaussian.
Since limits of Gaussians are Gaussian, the quantity $\dualprod{\phi}{f}$ will be a Gaussian scalar.
Therefore, a Gaussian in the sense of duality can loosely be thought of as a Gaussian process whose integral against arbitrary test functions is always Gaussian—an adaptation of the dot product notion encountered previously to the infinite-dimensional setting.

Of course, the above only describes how one can intuitively think about Gaussians in the sense of duality.
For general Gaussians, there is no defined notion of an integral—only of a dual pairing.
Such a Gaussian also need not be a random real-valued function: it's possible to consider random distributions and other objects generalizing the usual notion of a function.
It is therefore clear this view is very general.

This generality comes with a price: the resulting random vectors become more abstract, and it is much more difficult to understand whether or not they actually exist or say anything useful about them.
This difficulty is handled by specializing to less-general settings: for example, bogachev98 studies the setting where $V$ is a locally convex topological vector space and $W$ is its dual, and hairer09,lifshits12 study cases where $V$ is a Banach or Hilbert space.
In those cases, a number of results are available.

Provided we are considering a Gaussian $f$ which does exist, what objects play the role of a mean and covariance, and uniquely characterize $f$?
Since the canonical pairing $\dualprod{\.}{\.}$ is linear and non-degenerate, we know that $f = g$ holds if and only if $\dualprod{\phi}{f} = \dualprod{\phi}{g}$ for all $\phi\in W$, with both equalities used in the same sense.
It is therefore clear that the mean of a Gaussian process is simply a vector $\mu\in V$, since such a vector uniquely determines all expectations $\E\dualprod{\phi}{f}$.
The covariance is only slightly more subtle.

We say that a symmetric positive semi-definite bilinear form $k : W \x W \-> \R$ is a covariance form.

As before, we can construct a covariance form from a given Gaussian process.

The covariance form of a Gaussian in the sense of duality is defined by

\[
k(\phi,\psi) = \Cov(\dualprod{\phi}{f},\dualprod{\psi}{f})
\]

It is clear that two Gaussian processes $f$ and $g$ are equal in distribution if and only if they have the same mean and covariance form.
This can be seen by noting these requirements force $\dualprod{\phi}{f} = \dualprod{\phi}{g}$ to hold for all $\phi\in W$.

We now ask: what relationship does the covariance form have with the covariance kernel defined previously?
Consider again our Gaussian process $f$ defined on the space of continuous functions.
For any pair of test functions $\phi,\psi \in C^\infty([0,1];\R)$, the symmetric positive semi-definite bilinear form

\[
(\phi,\psi) \|> \int_0^1\int_0^1 \phi(x)k(x,x')\psi(x') \d x \d x'
\]

is the covariance form of $f$.

So far, this perspective mirrors the preceding ones, albeit with a more abstract presentation.
We have only introduced definitions.
These definitions, however, suffice to recover a notion of affine maps.

Let $f \~[GP](\mu, k)$. For $\c{A} : V \-> V'$ and $b\in V'$ we have

\[
\c{A} f + b \~[GP](\c{A}\mu + b, k(\c{A}^*(\.),\c{A}^*(\.)))
\]

where $V'$ and $W'$ are a dual pair, and $\c{A}^* : W' \-> W$ is the adjoint operator defined by $\dualprod{\c{A}^*\phi}{v} = \dualprod{\phi}{\c{A}v}$.

Immediate by definition.

Suppose now that $V$ is a locally convex topological vector space, and $W = V^*$ its topological dual.
Since the covariance form is a map $k : V^* \x V^* \-> \R$, it gives rise to an operator $\c{K} : V^* \-> V^{**}$, called the covariance operator, which in the case that $V$ is reflexive becomes an operator $\c{K} : V^* \-> V$ via composition with the canonical isomorphism given by reflexivity.
For our running example, this operator is given by

\[
\phi \|> \int_0^1 \phi(x)k(x,\.) \d x
\]

One can easily see that the transformation rule for covariance operators under affine maps is given by

\[
\c{K} \|> \c{A} \c{K}\c{A}^*
\]

If it is often unclear whether or not a Gaussian process with a given covariance form exists, it can be even less clear whether or not a Gaussian process with a given covariance operator exists.
Still, this notion is important, as these operators can potentially be studied using spectral theory.
We now move to the final concept we consider at this stage: that of a standard Gaussian.

Let $W \subseteq H$ where $H$ is a Hilbert space.
We say that $\c{W} : \Omega \-> V$ is a Gaussian white noise process if its covariance form coincides with the inner product in all cases where the former is defined.

In cases where $\c{K}$ admits a sufficiently rich spectral theory, this can potentially provide a suitable notion of an operator square root which relates Gaussian white noise to other Gaussian processes.

With these definitions in hand, we have recovered the two technical notions lost when transitioning from multivariate Gaussians to Gaussian processes in the sense of random functions.
The introduced machinery gives a way of constructing Gaussian processes that does not rely on defining a kernel: (i) construct a Gaussian white noise process, and (ii) define the Gaussian process of interest as an affine map of white noise.

We employ a variation of this strategy in <Ref> to construct Gaussian processes on Riemannian manifolds, where defining a kernel directly leads to difficulties.
Loosely speaking, we do so by taking $V$ to be a space of distributions, and $\c{A}$ to be the inverse of a differential operator.
This is chosen to ensure that $\c{A}$ smooths its inputs, so that its output is regular enough to be understood as a random function.

This line of thought leads one to the theory of stochastic partial differential equations driven by Gaussian white noise.
The main point swept here under the rug is that there is a convenient way to sidestep much of the analysis needed to carry out the above calculations using the theory of reproducing kernel Hilbert spaces.
Roughly, this amounts to instead working with certain vector spaces that uniquely determine the Gaussian processes of interest.

We also use the theory of Gaussians in the sense of duality for a second purpose: to construct a coordinate-free notion of Gaussianity as a suitable building block for constructing Gaussian vector fields on Riemannian manifolds.
The issue here is that the Gaussian process cannot be understood as a real-valued random function due to topological obstructions—it is instead a random section, which we define in <Ref>.
Finite-dimensional Gaussians in the sense of duality end up being the right tool for this setting.

For this, we prove a general existence theorem on Gaussians in the sense of duality in finite-dimensional settings.
Aside from giving an intrinsic reinterpretation of the multivariate Gaussians described previously, this ensures our framework is suited for its purpose in the coordinate-free setting.

Let $V$ be a finite-dimensional real topological vector space, and let $W = V^*$ be its topological dual.
Then for any vector $\mu \in V$ and covariance form $k : V^* \x V^* \-> \R$, there exists a unique-in-distribution random vector $y \~[N](\mu, k)$.

To prove this, choose a basis $e_i$ on $V$, and let $e^i$ be the dual basis.
Let $\c{E} : V \-> \R^d$ be the continuous linear isomorphism induced by the basis, and define $y = \c{E}^{-1} \v{y}$ with $\v{y}\~[N](\v\mu,\m{K})$ defined by

\[
\v\mu &= \begin{bmatrix}
\dualprod[0]{e^1}{\mu}
\\
\vdots
\\
\dualprod[0]{e^d}{\mu}
\end{bmatrix}
\m{K} &= \begin{bmatrix}
k(e^1, e^1) & \hdots & k(e^1, e^d)
\\
\vdots & \ddots & \vdots
\\
k(e^d, e^1) & \hdots & k(e^d, e^d)
\end{bmatrix}
\]

It is clear by direct calculation using Gaussians on $\R^d$ that the resulting vector is Gaussian with the right mean and covariance form.
The claim follows by noting that the assumed non-degeneracy of the dual pairing forces the distribution of every Gaussian in the sense of duality to be uniquely determined by its mean and covariance form.

Here, we see the key difference between the coordinate-free view and the matrix-vector view considered previously: the Gaussian random vector can be viewed as a real-valued multivariate Gaussian in any basis, but itself is defined on $V$ independent of this choice.
The value of considering this distinction in the first place will become clear in <Ref>.
Here, finite-dimensionality suffices to ensure existence: the story in infinite-dimensional settings is completely different and requires case-by-case analysis.

To conclude, we reflect on the introduced ideas.
We began by studying Gaussian random variables and random vectors, before generalizing these to Gaussian random functions, which offered a concrete framework where certain aspects of Gaussianity were seemingly lost.
Adopting a function-analytic view restored these aspects, at cost of increased abstraction.
This also gave a coordinate-free way to reinterpret the preceding multivariate constructions.

While much of the technical power of the duality framework is extraneous for our purposes, studying it nonetheless helps provide a unified conceptual perspective from which to interpret our developments.
By considering these notions, it becomes much clearer how one should understand the constructions encountered later, which might otherwise appear as if they arise out of thin air.
This completes our study of Gaussianity for its own sake, independent of Bayesian learning and other machine-learning-related considerations.

§ DISCUSSION

The preceding sections paint a rich and detailed picture of what a mathematical theory of decision-making under uncertainty looks like.
We now recap the steps taken so far, and reflect on them, before proceeding to describe contributions to be presented.

We began with the concept of probability, constructed in the language of measure theory.
We used this language to formalize the concept of learning via the notion of conditional probability, thereby obtaining the theory of Bayesian learning.
By working in an abstract measure-theoretic setting, we obtained a formalism suitable for learning about very general unknown quantities of interest.

We then took a step back, examining how to formalize the notion of an agent selecting actions in an unknown environment on basis of interactions, obtaining the key concept of a Markov decision process.
We then immediately restricted to the simpler setting of multi-armed bandits.
We saw that model-based algorithms built atop Bayesian learning yielded decision systems that perform effectively.
With these notions, we described how to efficiently solve global optimization problems using Bayesian optimization.

To transform the preceding ideas into a workable class of methods, we proceeded to study Gaussian process models in depth.
We developed these models in sequence, starting from the simplest settings, and ending with the highest generality.
These ideas provide us with key tools for understanding Gaussian processes, so that we can use them as building blocks of Bayesian models and high-performance decision systems atop those models.

It is worth pausing to reflect on the merits of the setting chosen, within a broader context of artificial intelligence.
In choosing to work with Bayesian learning, we opted to represent uncertainty using probability—a powerful but computationally limiting choice.
This choice was counterbalanced by working with simple models in bandit-like settings, and is most effective when the decisions of interest must be made in a data-efficient manner that only algorithms with near-asymptotically-optimal regret can achieve.

Not all settings fit these criteria well.
In many reinforcement learning problems of interest in robotics, the complexity of the dynamics—which, for multi-armed bandits, are totally absent from the problem—is a key difficulty.
Gaussian processes are largely not expressive enough to represent multi-object collision dynamics and related phenomena.
We have also not addressed partial observations—another key difficulty in that setting.

On the other hand, no other currently known theoretical framework comes close to understanding decision to the degree of command one can obtain from the notions described.
In the absence of a probabilistic framework, it is unclear how to assess, represent, and propagate uncertainty to resolve explore-exploit tradeoffs and minimize regret in non-trivial settings.
Thus, when probabilistic methods can be considered, they absolutely should be.

As a step towards building increasingly sophisticated decision systems, it seems fruitful to expand probabilistic approaches built via Bayesian learning to more general settings.
Improved understanding of these phenomena may yield lessons of broad interest to the application of decision systems.
Contributions presented here include development of pathwise conditioning methods for making Gaussian process models easier to work with, and a variety of non-Euclidean Gaussian processes, both described next.

CHAPTER: PATHWISE CONDITIONING

Gaussian processes admit analytic conditional distributions, making them a key model class for Bayesian learning.
The standard way of viewing these conditional distributions mirrors the general measure-theoretic setup common to all Bayesian models, and has strongly influenced how people think about Gaussian processes.

In the early 1970s, an alternative view emerged in the geostatistics community.
Miraculously, in the Gaussian case it is also possible to develop conditioning in a manner not purely based on distributions, but on random variables directly.
This view turns out to lift from the multivariate to the Gaussian process setting, yielding pathwise representations of posterior Gaussian processes, which have largely been overlooked in machine learning until now.

The pathwise perspective turns out to be a powerful point of view with wide-ranging consequences.
We will show how to use it to resolve a long-standing difficulty in Bayesian optimization: constructing a posterior approximation whose computational cost is linear both at training time and at test time, with excellent approximation properties and error control.

One of the key ingredients used within the construction will be basis function expansions of prior Gaussian processes.
We will thus examine a number of methods for constructing such expansions for different classes of priors.
We will also reinterpret sparse approximations in a function-based manner simpler than the typical viewpoint.
We conclude by benchmarking Bayesian optimization using pathwise sampling.
We proceed to these developments.

§ CONDITIONING MULTIVARIATE GAUSSIANS

We now describe conditioning of multivariate Gaussians.
Recall that using Bayes' Rule, a prior and likelihood combine into a joint distribution, which factorizes into the marginal distribution of the data and the posterior.
The posterior is the conditional distribution of the parameters given the data, which is unique almost everywhere with respect to the marginal distribution.
We now study how to represent this distribution for the case of interest.

§.§ Distributional conditioning

The most obvious way to represent a Gaussian conditional distribution is to simply calculate its distribution as a closed-form analytic expression.
This is given below.

\[
\begin{bmatrix}
\v\theta
\\
\v{y}
\end{bmatrix}
\~[N]\del{
\begin{bmatrix}
\v\mu_{\v\theta}
\\
\v\mu_{\v{y}}
\end{bmatrix}
\begin{bmatrix}
\m\Sigma_{\v\theta\v\theta} & \m\Sigma_{\v\theta\v{y}}
\\
\m\Sigma_{\v{y}\v\theta} & \m\Sigma_{\v{y}\v{y}}
\end{bmatrix}
\]

be non-singular.
Then we have that

\[
(\v\theta\given\v{y})(\.,\v\gamma) \~[N]\del{\v\mu_{\v\theta} + \m\Sigma_{\v\theta\v{y}}\m\Sigma_{\v{y}\v{y}}^{-1}(\v\gamma-\v\mu_{\v{y}}), \m\Sigma_{\v\theta\v\theta} - \m\Sigma_{\v\theta\v{y}}\m\Sigma_{\v{y}\v{y}}^{-1}\m\Sigma_{\v{y}\v\theta}}
\]

By non-singularity, $(\v\theta,\v{y})$ admits a Lebesgue density, and the claim follows by direct calculation via applying Bayes' Rule for densities.

Calculate conditional

Draw samples

[Distributional conditioning of multivariate Gaussians]Illustration of distributional conditioning of a bivariate Gaussian. Here, we form the joint distribution, calculate the conditional distribution, and then draw samples from it. Note that all steps except the very last one are distributional in nature and do not involve the use of random variables, which only appear at the very end of the process.

See [Appendix A]rasmussen06, for a full derivation.
The non-singularity requirement here is more-or-less necessary: otherwise, the marginal distribution of $\v{y}$ may admit too many null sets, rendering the desired conditional distribution non-unique, except in regions where one can apply a linear map to recover a suitable Lebesgue density and apply the above argument on the obtained subspace.

If one wants to work with this expression numerically, then it's possible to simply compute the desired conditional mean and covariance.
Then, one can generate conditional samples via the expression

\[
(\v\theta\given\v{y})(\omega,\v\gamma) &= \m{L}_{\v\theta\given\v{y}}\v{z}(\omega) + \v\mu_{\v\theta\given\v{y}}
\v{z} &\~[N](\v{0},\m{I})
\]

where $\v\mu_{\v\theta\given\v{y}}$ is the conditional mean, and $\m{L}_{\v\theta\given\v{y}}$ is a Cholesky factor of the conditional covariance, or any other matrix square root obtained numerically.
We illustrate this procedure in <Ref>.
This enables one to calculate any quantity of interest depending on the conditional distribution numerically via the Monte Carlo method.
The computational costs will in general be cubic in the dimension of both $\v\theta$ and $\v{y}$, owing to the need to compute $\m{L}_{\v\theta\given\v{y}}$ and invert $\m\Sigma_{\v{y}\v{y}}$, respectively.

§.§ Pathwise conditioning

The preceding considerations gave closed-form analytic expressions for Gaussian conditionals in terms of matrix-vector expressions that can be computed numerically.
From this, one might be tempted to conclude that there is nothing more to say about conditioning multivariate Gaussians—this, however, would miss an alternative view: Gaussian conditionals, which in general are a purely distributional notion, can also be described in a pathwise manner using random variables.

[Matheron's update rule]theoremthmmvnpw

For $\v\theta,\v{y}$ defined in <Ref>, we have that

\[
(\v\theta\given\v{y})(\omega,\v\gamma) = \v\theta(\omega) + \m\Sigma_{\v\theta\v{y}}\m\Sigma_{\v{y}\v{y}}^{-1}(\v\gamma - \v{y}(\omega))
\]

By direct calculation,

\[
\E(\v\theta + \m\Sigma_{\v\theta\v{y}}\m\Sigma_{\v{y}\v{y}}^{-1}(\v\gamma - \v{y})) = \v\mu_{\v\theta} + \m\Sigma_{\v\theta\v{y}}\m\Sigma_{\v{y}\v{y}}^{-1}(\v\gamma - \v\mu_{\v{y}}) = \E(\v\theta\given\v{y}=\v\gamma)
\]

\[
\Cov(\v\theta + \m\Sigma_{\v\theta\v{y}}\m\Sigma_{\v{y}\v{y}}^{-1}(\v\gamma - \v{y})) &= \m\Sigma_{\v\theta\v\theta} + \m\Sigma_{\v\theta\v{y}}\m\Sigma_{\v{y}\v{y}}^{-1}  \m\Sigma_{\v{y}\v\theta} - 2\m\Sigma_{\v\theta\v{y}}\m\Sigma_{\v{y}\v{y}}^{-1} \m\Sigma_{\v{y}\v\theta}
\\
&= \m\Sigma_{\v\theta\v\theta} - \m\Sigma_{\v\theta\v{y}}\m\Sigma_{\v{y}\v{y}}^{-1}  \m\Sigma_{\v{y}\v\theta} = \Cov(\v\theta\given\v{y}=\v\gamma)
\]

where we have cancelled a factor of $\m\Sigma_{\v{y}\v{y}}\m\Sigma_{\v{y}\v{y}}^{-1}$ in the middle term.

This is illustrated in <Ref>.
The above argument affirms the claim, but gives few hints on where this expression originates or how to obtain it from first principles.
To better understand this, we now prove <Ref> in a different way.
To do so, we first prove a plug-in property of conditioning.

Sample jointly

Transform into conditional

[Pathwise conditioning of multivariate Gaussians]Illustration of pathwise conditioning of a bivariate Gaussian. Here, we first sample a random vector from the joint distribution, then transform it into a sample from the conditional distribution.
These steps are called pathwise because they are defined directly using random variables, rather than indirectly through probability distributions.

Consider three random vectors $\v{a} : \Omega \-> \R^m$, $\v{b} : \Omega \-> \R^n$, and $\v{c} : \Omega \-> \R^m$ such that

\[
\v{a} = f(\v{b}) + \v{c}
\]

where $f : \R^n \-> \R^m$ is a measurable function, and where the random variables $\v{b}$ and $\v{c}$ are independent.
Then we have

\[
(\v{a} \given \v{b} = \v\beta) = f(\v\beta) + \v{c}
\]

This follows by direct calculation by writing

\[
\int_{A_{\v{b}}} \pi_{\v{a}\given\v{b}}(A_{\v{a}}\given\v\beta) \d\pi_{\v{b}}(\v\beta) &= \P(\v{a} \in A_{\v{a}}, \v{b} \in A_{\v{b}})
\\
&= \P(f(\v{b}) + \v{c} \in A_{\v{a}}, \v{b} \in A_{\v{b}})
\\
&= \int_{\R^m\x\R^n} \1_{f(\v\beta) + \v\varsigma \in A_{\v{b}}, \v\beta\in A_{\v{b}}} \d(\pi_{\v{c}}\ox\pi_{\v{b}})(\v\varsigma,\v\beta)
\\
&= \int_{A_{\v{b}}} \int_{\R^m} \1_{f(\v\beta) + \v\varsigma \in A_{\v{b}}} \d\pi_{\v{c}}(\v\varsigma) \d\pi_{\v{b}}(\v\beta)
\\
&= \int_{A_{\v{b}}} \P(f(\v\beta) + \v{c} \in A_a) \d\pi_{\v{b}}(\v\beta)
\]

where we have used independence to represent the probability as an integral over a product measure, followed by Tonelli's Theorem.
By the Disintegration Theorem, $\pi_{\v{a}\given\v{b}}$ is $\pi_{\v{b}}\ae[-]$ unique, implying that the conditional distributions of interest are equal, and the claim follows.

The key idea behind our argument will be to choose the conditional expectation for our function $f$, or more precisely, the map $\v{y} \|> \E(\v\theta\given\v{y})$.
We now recall this notion and some of its key properties.

Conditional expectation is defined as the orthogonal projection from the Lebesgue space $L^2(\Omega,\c{F},\P;\R^n)$ onto the subspace $L^2(\Omega,\sigma(\v{y}),\P;\R^n)$ where $\sigma(\v{y})$ is the smallest $\sigma$-algebra containing all preimages $\v{y}^{-1}(A_{\v{y}})$ where $A_{\v{y}}\in\c{B}(\R^n)$.
Recall that the preimage is a map $\v{y}^{-1} : \c{B}(\R^n) \-> \c{F}$ between $\sigma$-algebras.
This definition is reasonably intuitive: we can think of this as projecting onto the subspace induced by all collections of random numbers in $\Omega$ which play a role in determining what $\v{y}$ does.

Recall that $L^2(\Omega,\c{F},\P;\R^n)$ is the Hilbert space of equivalence classes of random variables with inner product given by $\innerprod{\v{a}}{\v{b}} = \E(\v{a}\cdot\v{b})$.
Using this, it follows from the Projection Theorem for Hilbert spaces that the terms $\E(\v\theta\given\v{y})$ and $(\v\theta - \E(\v\theta\given\v{y}))$ are uncorrelated—and, in the Gaussian case, that they are independent.
This gives us the candidate random variables to use for $\v{b}$ and $\v{c}$, if we choose conditional expectation for $f$.

Finally, recall that for multivariate Gaussians, the conditional expectation is given by $\E(\v\theta\given\v{y}) = \v\mu_{\v\theta} + \m\Sigma_{\v\theta\v{y}}\m\Sigma_{\v{y}\v{y}}^{-1}(\v{y} - \v\mu_{\v{y}})$, where we note in our setting that the inverse always exists by non-singularity of $\v{y}$.