Intro to Bayesian Optimization 4. Gaussian Processes

This is the fourth part in my Intro to Bayesian Optimization

1. Optimization Intuitions
2. The Bayesian Optimization Framework
3. Intro to Bayesian Statistics
4. Gaussian Processes

So far we have described the Bayesian Optimization framework, and have a basic understanding of Bayesian statistics. In BO, we have no knowledge about the function, so the model we use needs to be able to describe arbitrary continuous functions.

Gaussian Process regression is a nonparametric method to model arbitrary continuous functions given limited data. Sounds perfect for our situation. Most literature and implementations of BO use GPs for this reason. However the usefulness of GPs is wide-ranging, and the study of GPs is a large research area.

A change of perspective

Let’s say we have 5 data points. Each with an $x$ and a $y$ value. Let’s ignore the $x$ values for now and just think of our vector of $y$ values:

$$\begin{array}{c}[y_1,y_2,y_3,y_4,y_5]\end{array}$$

Now, we could think of it as 5 points in 1-dimensional space, but if we tilt our head, we could also interpret this as 1 point in 5-dimensional space.

What if all the data points in our infinite sample space are just a single point in infinite-dimensional space, With every possible point corresponding to its own dimension?

The Gaussian Process

Gaussian Process regressions take this perspective, then define a Gaussian distribution over this infinite-dimensional space, which we call our function space. A point in function space now corresponds to a realization of the data (a function) in our original data space. We pick a Gaussian Process as the prior distribution over the function space.

The specific definition is that a Gaussian Process is any (possibly infinite) collection of Random Variables, such that any finite collection of them have a joint Gaussian distribution. So here we define a GP over infinite-dimensional space, and any finite subset (e.g. our data points) is also Gaussian in its joint distribution.

Note that we are not saying that our $y$ values are gaussian with respect to the $x$ values. We are saying that the $y$ values themselves fall out of a Gaussian distribution with respect to each other. This will become clearer when we visualize the prior.

Kernels

Now think about what it means for a function to be continuous, it means an arbitrarily small change in our $x$ variables would result in a bounded change in our $y$ variable, i.e. a small $x$ change is a small $y$ change. We can encode this in our new space by saying that if our change in $x$ between two points is small, the correlation in the corresponding dimensions of the function space is large because the $y$ values should be similar. We call this mapping of distance of two points to covariance a kernel $k(x,x^{\prime })$ where $x$ and $x^{\prime }$ are two points.

A common kernel is the radial basis function (RBF) aka the exponentiated quadratic kernel:

$$\begin{array}{c}k(x,x^{\prime })=e^{-\frac{1}{2}\Vert x-x^{\prime }{\Vert }^2}\text{where }\Vert x-x^{\prime }\Vert \text{ is the euclidean distance between x and x′}\end{array}$$

When $x$ and $x^{\prime }$ are close, $k(x,x^{\prime })$ goes to 1, as they go far away, $k(x,x^{\prime })$ goes to 0. Note this works for any number of dimensions of explanatory variables. We actually include two hyperparameters into the RBF for tuning, ${\sigma }_f$ controls the amplitude and $l$ controls how fast covariance drops off:

$$\begin{array}{c}k(x,x^{\prime };{\sigma }_f,l)={\sigma }_f^2{e}^{-\frac{1}{2l^2}\Vert x-x^{\prime }{\Vert }^2}\end{array}$$

There are many other reasonable choices of kernel you can use depending on what you know about the function you are approximating (e.g. if your function is periodic and you want to encode that). Read more

So if we some response variable vector $\mathbf{y}=[y_1,y_2,\dots ,y_n]$ and some explanatory variables $\mathbf{X}=[x_1^T,x_2^T,\dots ,x_n^T]$, then our prior distribution for $\mathbf{y}$ is the Gaussian Process:

$$\begin{array}{c}\mathbf{y}=\left[\begin{array}{c}f(x_1)\\ f(x_2)\\ ⋮\\ f(x_n)\end{array}\right]\sim \mathcal{N}(0,K(\mathbf{X},\mathbf{X}))\text{ where }K(\mathbf{X},\mathbf{X})=\left[\begin{array}{cccc}k(x_1,x_1)& k(x_1,x_2)& \dots & k(x_1,x_n)\\ k(x_2,x_1)& k(x_2,x_2)& & k(x_2,x_n)\\ ⋮& & \ddots & ⋮\\ k(x_n,x_1)& k(x_n,x_2)& \dots & k(x_n,x_n)\end{array}\right]\end{array}$$

We may choose the mean not to be 0, but some function $\mu (x)$ in practice. But knowing nothing else about the function, f(x), a mean of 0 is a reasonable choice.

Sampling the prior

So now we have a distribution that describes our function space. Remember, our reasoning here was to find a pretty general way to model arbitrary continuous functions. If we were going to sample a point, what would the function look like?

The Gaussian part of Gaussian Processes refers to these wiggly curves. Note they are not shaped as Gaussians; we are not saying that the function has a Gaussian shape, we are saying that the points of the function are a priori one sample from an infinite-dimensional Gaussian distribution with a covariance described by the kernel. Some wiggly curves are more likely than others, note how all of these ones centre around 0, that is because of our choice of mean of the GP. The probability of the curves is Gaussian, and that’s what gives them their arbitrary continuous (wiggly) shape.

Consider the earlier example of Bayesian linear regression. That was in a sense Gaussian too, as we said our two parameters came from a joint normal distribution. Same thing here, except instead of parameters, we consider all infinitely many points on the real number line.

Visualizing the prior

We can take the limit of the number of samples and overlay them to visualize the probability distribution of the prior at each point $x$.

Posterior

We start with a very useful theorem. If we have a Gaussian Process, then by definition the joint distribution of any two subsets of the GP are also Gaussian. And it turns out the conditional distribution of Gaussians is also Gaussian.

Specifically, if you consider two subsets of random variables of a Gaussian Process, ${\mathbf{y}}_1$ and ${\mathbf{y}}_2$

$$\begin{array}{c}\left[\begin{array}{c}{\mathbf{y}}_1\\ {\mathbf{y}}_2\end{array}\right]\sim \mathcal{N}\left(\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\end{array}\right],\left[\begin{array}{cc}{\Sigma }_{11}& {\Sigma }_{12}\\ {\Sigma }_{21}& {\Sigma }_{22}\end{array}\right]\right)\end{array}$$

then

$$\begin{array}{c}p({\mathbf{y}}_2|{\mathbf{y}}_1)=\mathcal{N}\left({\mu }_{2|1},{\Sigma }_{2|1}\right)\end{array}$$

where

$$\begin{array}{c}{\mu }_{2|1}={\mu }_2+{\Sigma }_{21}{\Sigma }_{11}^{-1}({\mathbf{y}}_1-{\mu }_1)\text{and}{\Sigma }_{2|1}={\Sigma }_{22}-{\Sigma }_{21}{\Sigma }_{11}^{-1}{\Sigma }_{12}\end{array}$$

So if we have observed data ${\mathbf{y}}_{data},{\mathbf{X}}_{data}$, and want to calculate the posterior for prediction points ${\mathbf{y}}_{\ast },{\mathbf{X}}_{\ast }$, we have prior:

$$\begin{array}{c}\left[\begin{array}{c}{\mathbf{y}}_{data}\\ {\mathbf{y}}_{\ast }\end{array}\right]\sim \mathcal{N}\left(\left[\begin{array}{c}0\\ 0\end{array}\right],\left[\begin{array}{cc}K& K_{\ast }\\ K_{\ast }^T& K_{\ast \ast }\end{array}\right]\right)\end{array}$$

where

$$\begin{array}{c}K=K({\mathbf{X}}_{data},{\mathbf{X}}_{data}),K_{\ast }=K({\mathbf{X}}_{data},{\mathbf{X}}_{\ast }),K_{\ast \ast }=K({\mathbf{X}}_{\ast },{\mathbf{X}}_{\ast })\end{array}$$

then we have the posterior:

$$\begin{array}{c}p({\mathbf{y}}_{\ast }|{\mathbf{y}}_{data})=\mathcal{N}\left(K_{\ast }^T{K}^{-1}{\mathbf{y}}_{data},K_{\ast \ast }-K_{\ast }^T{K}^{-1}{K}_{\ast }\right)\end{array}$$

Note that ${\mathbf{y}}_{\ast }$ doesn’t need to come from a small or even finite set, it can be the entire infinite-dimensional space which is equivalent to every point on the number line.

Sampling the posterior

Sampling a few times

Sampling many times

As we have a closed form solution for the posterior (one of the nice things about GPs), we can visualize it without sampling. Just specify your level of credibility.

Noisy observations

Say instead we have $y=f(x)+\epsilon$, where measurements have error

$$\begin{array}{c}\epsilon \sim \mathcal{N}(0,{\sigma }_m^2)\end{array}$$

Then the prior covariance matrix of the data points is

$$\begin{array}{c}K=K({\mathbf{X}}_{data},{\mathbf{X}}_{data})+{\sigma }_m^{2}I\end{array}$$

and we can just replace this in the earlier process.

Note that this gives the posterior distribution, i.e. the distribution for the noiseless function, not the predictive distribution for future measurements which will need to have noise added.

$$\begin{array}{c}p\left(f\left({\mathbf{X}}_{\ast }\right)|{\mathbf{y}}_{data}\right)\end{array}$$

In the noiseless case, we called $p\left({\mathbf{y}}_{\ast }|{\mathbf{y}}_{data}\right)$ the posterior, but that was just because of the special case where $y=f(x)$

Prior on the functions? Where are the parameters?

In the linear regression case, we defined a prior of the parameters, not the functions that they create. Here, we set the prior on the functions directly:

$$\begin{array}{c}p(f(x)|y)=\frac{p(y|f(x))\cdot p(f(x))}{p(y)}\end{array}$$

That is because Gaussian Process regression is a nonparametric method. This is a catch-all term for any model that doesn’t use specified models with parameters that can be optimized — including histograms, decision trees, and some support vector machines. However, often when we use this term, what we are actually considering is infinitely many parameters.

With a GP, the complexity of the model (i.e. its number of dimensions) increases with the number of data points. Contrast this with linear regression where complexity increases with the number of explanatory variables. Scaling complexity with the number of data points has a huge advantage in low data models. Consider a linear regression, where our fit can be characterized just by the coefficients. No matter how many data points we have, in a one-dimensional case, we can compress the information of the model down into an intercept and coefficient term. However, with GPs, we are never compressing/losing information, in a sense, we are maximally leveraging the data we have available.

The flip side of this is that GPs are computationally constrained to low data regimes. We have won statistical efficiency at the expense of computational efficiency. Practically speaking, we must calculate the inversion of the Kernel matrix of the observed data $K^{-1}$, which is an expensive $O(n^3)$ operation where $n$ is the number of data points. Compare this with a linear regression where we must calculate $({\mathbf{X}}^T\mathbf{X}{)}^{-1}$ which is an $O(p^3)$ operation, where $p$ is the number of parameters. With GPs we are constrained to working at about 10,000 data points given modern tooling. However, this is fine for Bayesian Optimization where we don’t have so many data points. Also, we do not require this model to be computationally performant, as most time should be spent in function evaluation, i.e. performing experiments. It’s ok for modelling to take a few minutes if the time between experiments is a few hours.

Hyperparameter selection

We have introduced 3 hyperparameters. The kernel has terms $l$ and ${\sigma }_f$, and with noise we introduced ${\sigma }_m$. A common method used to pick these is to maximize the marginal likelihood $\arg {\max }_{l,{\sigma }_f,{\sigma }_m}p(\mathbf{y}|\mathbf{X},l,{\sigma }_f,{\sigma }_m)$ which can be done through gradient methods.

However, there are arguments to be made in the context of Bayesian Optimization that we don’t want to find optimal model hyperparameters, as we may want to err on the side of overestimating the kernel variance so we don’t get stuck in local minima through the iterative process.

Summary

We want to find the combination of $x$ values that gives us an optimal $y$ value, where $y=f(x)+\epsilon$. $f$ is a black box we can sample, but we cannot sample its derivatives. Sampling is expensive so we want to minimize the number of samples we take.

The approach we take is to create a probabilistic (Bayesian) model of the function with the data we have, then find the right balance of exploration and exploitation to pick the next point to sample (maximizing the acquisition function).

We use a Gaussian Process for this probabilistic model, because it can model arbitrary continuous functions, and is powerful in low data regimes.

1. Start with some initial sample measurements
2. Fit a Gaussian Process to the data
3. Compute the acquisition function based on the GP
4. Find $x^{\ast }$ that maximizes the acquisition function
5. Sample $f(x^{\ast })$ and add to your data
6. If time/resources are remaining, go to step 2