定义：若对于任意时刻ti(i=1,2,...,n)，随机过程的任意n维随机变量Xi=X(ti)(i=1,2,...,n)服从高斯分布，则称X(t)为高斯随机过程或正太过程。

The key idea behind GPs is that a function can be modeled using an infinite dimensional multivariate Gaussian distribution. In other words, every point in the input space is associated with a random variable and the joint distribution of these is modeled as a multivariate Gaussian.

Ok, so what does that mean and what does it actually look like? Well lets start with a simpler case: a unit 2D Gaussian. How can we start to view this as a distribution over functions? Here’s what we have:

Normally this is visualised as a 3D bell curve with the probability density represented as height. But what if, instead of visualising the whole distribution, we just sample from the distribution. Then we will have two values which we can plot on an graph. Let’s do this 10 times, putting the first value at x=0x=0 and the second at x=1x=1 and then drawing a line between them.

Looking at all these lines on a graph, it starts to look like we’ve just sampled 10 linear functions… What if we now use a 20-dimensional Gaussian, joining each of the sampled points in order?

These definitely look like functions of sorts but they look far too noisy to be useful for our purposes. Let’s think a bit more about what we want from these samples and how we might change the distribution to get better looking samples…

The multivariate Gaussian has two parameters, its mean and covariance matrix. If we changed the mean then we would only change the overall trend (i.e. if the mean was ascending integers, e.g. np.arange(D) then the samples would have an overall positive linear trend) but there would still be that jagged noisy shape. For this reason we tend to leave the GP mean as zero - they actually turn out to be powerful enough to model many functions without changing this.

Instead we want some notion of smoothness: i.e. if two input points are close to each other then we expect the value of the function at those points to be similar. In terms of our model: random variables corresponding to nearby points should have similar values when sampled under their joint distribution (i.e. high covariance).

The covariance of these points is defined in the covariance matrix of the Gaussian. Suppose we have an NN dimensional Gaussian modeling y0,…,yNy0,…,yN, then the covariance matrix ΣΣ is N×NN×N and its (i,j)(i,j)-th element is Σij=cov(yi,yj)Σij=cov(yi,yj). In other words ΣΣ is symmetric and stores the pairwise covariances of all the jointly modeled random variables.

So how should we define our covariance function? This is where the vast literature on kernels comes in handy. For our purposes we will choose a squared exponential kernel which (in its simplest form) is defined by:

This function (which we plot in a moment) is 1 when x=x′x=x′ and tends to zero as its arguments drift apart.

We can plot this kernel to show how it’s maximised when x=x′x=x′ and then smoothly falls off as the two inputs start to differ.

So, to get the sort of smoothness we want we will consider two random variables yiyi and yjyjplotted at xixi and xjxj to have covariance cov(yi,yj)=κ(xi,xj)cov(yi,yj)=κ(xi,xj) – the closer they are together the higher their covariance.

Using the kernel function from above we can get this matrix with k(xs, xs). Now lets try plotting another 10 samples from the 20D Gaussian, but this time with the new covariance matrix. When we do this we get:

Now we have something thats starting to look like a distribution over (useful) functions! And we can see how as the number of dimensions tends to infinity we don’t have to connect points any more because we will have a point for every possible choice of input.

Let’s use more dimensions and see what it looks like across a bigger range of inputs:

Now that we have a distribution over functions, how can we use training data to model a hidden function so that we can make predictions?

First of all we need some training data! And to get that we are going to create our secret function ff.

For this intro we’ll use a 5th order polynomial:

I chose this because it has a nice wiggly graph but we could have chosen anything.

Now we get to the heart of GPs. Theres a bit more maths required but it only consists of consolidating what we have so far and using one trick to condition our joint distribution on observed data.

So far we have a way to model p(y|x)p(y|x) using a multivariate normal:

where K=κ(x,x)K=κ(x,x) and m(x)=0m(x)=0.

This is a prior distribution representing the kind out outputs yy that we expect to see over some inputs xx before we observe any data.

So we have some training data with inputs xx, and outputs y=f(x)y=f(x). Now lets say we have some new points x∗x∗ where we want to predict y∗=f(x∗)y∗=f(x∗).

Now recalling the definition of a GP, we will model the joint distribution of all of yy and y∗y∗as:

where K=κ(x,x)K=κ(x,x), K∗=κ(x,x∗)K∗=κ(x,x∗) and K∗∗=κ(x∗,x∗)K∗∗=κ(x∗,x∗). As before we are going to stick with a zero mean.

However this is modeling p(y,y∗|x,x∗)p(y,y∗|x,x∗) and we only want a distribution over y∗y∗!

Rather than deriving it from scratch we can just make use of this standard result. If we have a joint distribution over yy and y∗y∗ as above, and we want to condition on the data we have for yy then we have the following:

Now we have a posterior distribution over y∗y∗ using a prior distribution and some observations!

NB: The code below would not be used in practice since KK can often be poorly conditioned, so its inverse might be inaccurate. A better approach is covered in part II of this guide!

And that’s it! We can now use these two parameters to draw samples from the conditional distribution. Here we plot them against the true function f(x)f(x) (the dashed black line). Since we are using a GP we also have uncertainty information in the form of the variance of each random variable. We know the variance of the ii-th random variable will be Σ∗iiΣ∗ii - in other words the variances are just the diagonal elements of Σ∗Σ∗. Here we plot the samples with an uncertainty of ±2±2 standard deviations.

In practice we need to do a bit more work to get good predictions. You may have noticed that the kernel implementation contained two parameters σσ and ll. If you try changing those when sampling from the prior then you can see how σσ changes the vertical variation and ll changes the horizontal scale. So we would need to change these to reflect our prior belief about the hidden function ff. For instance if we expect ff to have a much bigger range of outputs (for the domain we are interested in) then we would need to scale up σσaccordingly (try scaling the return value of ff by 100 to see what happens, then setsigma=100). In fact, as with anything that uses kernels, we might change our kernel entirely if we expect a different kind of function (e.g. a periodic function).

Picking the kernel is up to a human expert but choosing the parameters can be done automatically by minimising a loss term. This is the realm of Gaussian process regression.

Finally we should consider how to handle noisy data - i.e. when we can’t get perfect samples of the hidden function ff. In this case we need to factor this uncertainty into the model to get better generalisation.

These two topics will be the focus of Introduction to Gaussian Processes - Part II.