Starting to see these everywhere, so it may useful to write down a brief set of notes. These notes are based on Gaussian Processes for Machine Learning (Rasmussen and Williams, 2005).We’ll motivate Gaussian Processes (GPs) with linear regression. Consider a supervised learning task where we attempt to model the relationship in a dataset 𝐷={(𝑥𝑖,𝑦𝑖)}𝑁𝑖=1 with data points and their corresponding labels. We are essentially trying to find the joint relationship 𝑝(𝑥,𝑦). We can do this as either:𝑝(𝑥,𝑦)=𝑝(𝑦)𝑝(𝑥|𝑦)generative=𝑝(𝑥)𝑝(𝑦|𝑥)discriminativeWhy generative and discriminative? In the first instance, we’re modeling 𝑝(𝑥|𝑦), which takes the labels and tries to generate the data. Conversely, discriminative models will take the data and try to discriminate it into the possible classes (the labels).Briefly for the generative case, we can model each class-conditional probability as𝑝(𝑥|𝐶𝑐)=𝒩︀(𝜇𝑐,Σ𝑐)where 𝐶𝑐 is the label for a particular class 𝑐 with an appropriate prior for 𝑝(𝑦). The issue here is we are making strong assumptions that given the label, the data takes a normal density. Mixture of Gaussians and other such models take this approach and optimize for correct 𝜇𝑐,Σ𝑐.We instead look at the discriminative case. A simple approach is to take a regression output and transform it into a class probability with a “response function”. We “squash” the regressed value from the domain of (−∞,∞) into the probabilistic range of [0,1].𝑝(𝐶1|𝑥)=𝜎(𝑥⊤𝑤)where 𝑤 are the learned weights of our regressor. Remember that in linear regression, we place a Gaussian prior on the weights and then look at the “observations” and express the posterior distribution over the weights with Bayes’ rule.𝑤∼𝒩︀(0,Σ𝑝)𝑝(𝑤|𝑋,𝑦)=𝑝(𝑦|𝑋,𝑤)𝑝(𝑤)𝑝(𝑦|𝑋)where 𝑤∈ℝ𝑁 (with 𝑝 meaning prior) and 𝑋∈ℝ𝑛×𝑑 contains all the data points. This problem can be solved with MAP.But let’s take a closer look at our linear model 𝑓(𝑥)=𝜙(𝑥)⊤𝑤 with prior 𝑤∼𝒩︀(0,Σ𝑝). First, we modify our original form by adding this “feature” transformation of our original 𝑥. An example of 𝜙 could be (𝑥𝑥2), with possibly infinite entries.We can input any point 𝑥∈ℝ𝑑 in our linear model and retrieve a density for the regressed value. Computing the mean and covariance for this distribution:𝔼[𝑓(𝑥)]=𝔼[𝜙(𝑥)⊤𝑤]=𝜙(𝑥)⊤𝔼[𝑤]=0Cov(𝑓(𝑥),𝑓(𝑥′))=𝔼[𝑓(𝑥)𝑓(𝑥′)]−𝔼[𝑓(𝑥)]𝔼[𝑓(𝑥′)]=𝔼[(𝜙(𝑥)⊤𝑤)(𝑤⊤𝜙(𝑥′))]=𝜙(𝑥)⊤𝔼[𝑤𝑤⊤]𝜙(𝑥′)=𝜙(𝑥)⊤Σ𝑝𝜙(𝑥′)Effectively, since 𝑤 is a Gaussian random vector and 𝑓(𝑥) is a linear transformation of 𝑤, any 𝑓(𝑥) and 𝑓(𝑥′) are jointly Gaussian with mean 0 and covariance 𝜙(𝑥)⊤Σ𝑝𝜙(𝑥′). Interestingly, any finite function values 𝑓(𝑥1),𝑓(𝑥2),…,𝑓(𝑥𝑛) will be jointly Gaussian!This is a Gaussian process: a collection of random variables where any finite number of which have a joint Gaussian distribution.Here, we assume 𝜙(𝑥) and 𝑤 are both finite-dimensional, which means the joint covariance matrix will have finite dimensions with rank 𝑁 1 . But if we have more than 𝑁 data points (𝑛>𝑁), then this rank 𝑁 restriction means the joint Gaussian distribution becomes singular. The limited degrees of freedom means that we can’t fit all the points and the Gaussian process is degenerate 2 .Instead, we’ll extend 𝑁→∞ and place a Gaussian bump at every real point 𝑐∈ℝ𝑑, make 𝑤 a function over 𝑐 and integrate over ℝ𝑑:𝜙(𝑥;𝑐)=exp(−‖𝑥−𝑐‖222ℓ2)𝑓(𝑥)=∫ℝ𝑑𝑤(𝑐)𝜙(𝑥;𝑐)d𝑐𝑤 is now Gaussian white noise, defined by𝔼[𝑤(𝑐)]=0𝔼[𝑤(𝑐)𝑤(𝑐′)]=𝜎2𝛿(𝑐−𝑐′)where the Dirac delta 𝛿 enforces independence across different bumps.Similar to before, we can compute the mean and covariance of our new 𝑓(𝑥) 3 as 0 and a kernel function:𝑘(𝑥,𝑥′)=𝜎2𝑓exp(−‖𝑥−𝑥′‖222ℓ2)This is the “radial basis function” (RBF) or scaled exponential (SE) commonly used in the Bayesian linear regression model.This formulation gives the general Gaussian process:𝑓(𝑥)∼𝒢︀𝒫︀(𝑚(𝑥),𝑘(𝑥,𝑥′))where 𝑚(𝑥)=𝔼[𝑓(𝑥)] and 𝑘(𝑥,𝑥′)=𝔼[(𝑓(𝑥)−𝑚(𝑥))(𝑓(𝑥′)−𝑚(𝑥′)).To evaluate such a process on test points 𝑋∗, we simply condition the joint Gaussian prior distribution on the observations as 𝑓∗|𝑋∗,𝑋,𝑓, which is another Gaussian distribution that we can sample from.Interestingly, we no longer need to define feature spaces 𝜙 explicitly, turning the distributions over weights into distributions over functions.