Gaussian Processes

Gaussian processes can be viewed through a Bayesian lens of the function space: rather than sampling over individual data points, we are now sampling over entire functions. They extend the idea of bayesian linear regression by introducing an infinite number of feature functions for the input XXX.

In geostatistics, Gaussian processes are referred to as kriging regressions, and many other models, such as Kalman Filters or radial basis function networks, can be understood as special cases of Gaussian processes. In this framework, certain functions are more likely than others, and we aim to model this probability distribution.

Gaussian processes retain epistemic uncertainty, meaning we are still modeling what we don't know about the function, while remaining tractable. Essentially, they provide a normal distribution over functions.

If we index this infinite dimensional Gaussian, thanks to the properties of Gaussians with marginalization, we still have another Gaussian. We need a way to parametrize now. In this case we do it using a covariance function $k(x,x^{\prime })=\text{Cov}(f(x),f(x^{\prime }))$

This seems a quite nice resource on the topic.

Formal definition of Gaussian Process

A Gaussian process is an infinite set of random variables such that any finite number of them are jointly Gaussian.

The informal ideas is that we use the probability of the joints, given by our training dataset, to do some inference.

A Gaussian Process is a (infinite) set of random variables, indexed by some set $X$ (domain of the function) and for each $x\in X$ there is a random variable $f_x$ such that for all $A\subseteq X$ and $A=\left\{x_1,\dots ,x_m\right\}$ it holds that $f_A=[f_{x_1},\dots ,f_{x_m}]\sim \mathcal{N}({\mu }_A,K_{AA})$ , this is called noise-free prediction, to have the value of $y$ just add some Gaussians with a noise variance. You can see that GPs are very sensitive to adversarial data, it will be considered in the prediction and will condition the mean and Covariance of the prediction. Where

$$K_{AA}=\left(\begin{array}{c}k(x_1,x_1),\dots ,k(x_1,k_m)\\ ⋮\\ k(x_m,x_1),\dots ,k(x_m,k_m)\end{array}\right)\text{ and }{\mu }_A=\left(\begin{array}{c}\mu (x_1)\\ ⋮\\ \mu (x_m)\end{array}\right)$$

$k:X\times X\to \mathbb{R}$ is called the covariance (kernel) function and $\mu :X\to \mathbb{R}$ is called the mean function. We write $f\sim GP(\mu ,k)$ and define $f(x):=f_x$ See Kernel Methods for more about kernels.

Inference for Gaussian Processes

Inference is easy if you understand the conditional Gaussian distribution, see Gaussians. You already have the old points defined, you just need to get the conditioned Gaussian on the previous data points.

Here is the setting, and it's very similar to the Bayesian Linear Regression (indeed, GP can be seen as aa generalization of BLR) function view inference

$$\bar{f}=\left[\begin{array}{c}y\\ f\ast \end{array}\right],\bar{\Sigma }=\left[\begin{array}{cc}{\Sigma }_{AA}+{\sigma }_n^{2}I& k(x^{\ast },A{)}^T\\ k(x^{\ast },A)& k(x^{\ast },x^{\ast })\end{array}\right],\mu =\left[\begin{array}{c}{\mu }_A\\ {\mu }_{f^{\ast }}\end{array}\right]$$

We just use the conditional Gaussian to have that our wanted distribution $f^{\ast }\mid f$ has the following distribution:

$$\{\begin{array}{l}\mu ={\mu }_{f^{\ast }}+k(x^{\ast },A{)}^T({\Sigma }_{AA}+{\sigma }_n^{2}I{)}^{-1}(A-{\mu }_A)\\ \Sigma =k(x^{\ast },x^{\ast })-k(x^{\ast },A{)}^T({\Sigma }_{AA}+{\sigma }_n^{2}I{)}^{-1}k(x^{\ast },A)\end{array}$$

And then you just sample with this, which we will solve in the next section. Where $k(x^{\ast },A)=\left[k(x^{\ast },x_1),\dots ,k(x^{\ast },x_n)\right]$

If we have only the GP prior then the mean ${\mu }_{f^{\ast }}=0$, this is often a convention because the prior mean is only used to compute the posterior mean, and can be easily added back in. One note: the posterior mean depends on the observations $y$ while the variance does not.

Sampling for Gaussian Processes

Standard Gaussian Sampling

One simple solution is just sampling using $X={\Sigma }^{1/2}\epsilon +\mu$, but taking the square root of a matrix takes $\mathcal{O}(n^3)$ time. This samples every needed point immediately. This is possible because by assumption we have that the GP is jointly Gaussian, meaning:

\begin{array} \\ f_{1:n} \sim \mathcal{N}(\mu_{f}, K_{f}) \\ \end{array}

Forward Sampling

Another approach, called forward sampling just samples everything one by one. Meaning

\begin{array} \\ f_{1} \sim p(f_{1}) \\ f_{2} \sim p(f_{2} \mid f_{1}) \\ \vdots \\ f_{n} \sim p(f_{n} \mid f_{1:n-1}) \end{array}

But we need to compute the matrix inverse in this case, which also takes $\mathcal{O}(n^3)$.

Learning Gaussians Processes

Kernel optimization

Once we have chosen a Kernel, we would like to find the best Hyper-parameters for our problem. Take a look at the Kernel Methods for possible Kernels. We want to optimize the following equation:

$$\begin{array}{rl}& \\ \widehat{{\theta }_{MLE}}& =\arg \underset{\theta }{\max }p(y_{\text{train}}\mid x_{\text{train}},\theta )\\ & =\arg \underset{\theta }{\max }\int p(y_{\text{train}}\mid f,x_{\text{train}},\theta )p(f\mid \theta )df\\ & =\arg \underset{\theta }{\max }\int p(y_{1:n}\mid f_{1:n},x_{1:n},\theta )p(f_{1:n}\mid \theta )df\\ & =\arg \underset{\theta }{\max }\int \prod _{i=1}^n\mathcal{N}(y_i\mid f_i;0,{\sigma }_n^2)\mathcal{N}(f_{1:n};0,K_f(\theta ))df\\ & =\arg \underset{\theta }{\max }\mathcal{N}(y_{1:n};0,K_y(\theta ))\\ & =\arg \underset{\theta }{\max }\frac{1}{\sqrt{(2\pi {)}^n|K_y(\theta )|}}\exp \left(-\frac{1}{2}{y}_{1:n}^T{K}_y^{-1}(\theta )y_{1:n}\right)\end{array}$$

Now let's take the negative log and we will get:

$$\begin{array}{rl}& \\ \widehat{{\theta }_{MLE}}& =\arg \underset{\theta }{\min }-\log \mathcal{N}(y_{1:n};0,K_y(\theta ))\\ & =\arg \underset{\theta }{\min }\frac{n}{2}\log 2\pi +\frac{1}{2}\log |K_y(\theta )|+\frac{1}{2}{y}_{1:n}^T{K}_y^{-1}(\theta )y_{1:n}\\ & =\arg \underset{\theta }{\min }\frac{1}{2}\log |K_y(\theta )|+\frac{1}{2}{y}_{1:n}^T{K}_y^{-1}(\theta )y_{1:n}\end{array}$$

Where

$$K_f(\theta )=K_{f_{1:n}}(\theta )=\left[\begin{array}{ccc}{k}_{\theta }(x_1,x_1),& \dots & ,k_{\theta }(x_1,x_n)\\ ⋮& & ⋮\\ k_{\theta }(x_n,x_1),& \dots & ,k_{\theta }(x_n,x_n)\end{array}\right]$$

And

$$K_y(\theta )=K_f(\theta )+{\sigma }_n^{2}I$$

One can also prove that the gradient of our likelihood function is

$$\frac{\partial \log p(y_{1:n}\mid x_{1:n},\theta )}{\partial {\theta }_j}=\frac{1}{2}\text{tr}\left((\alpha {\alpha }^T-K_{y,\theta }^{-1})\frac{\partial K_{y,\theta }}{\partial {\theta }_j}\right)$$

With $\alpha =K_{y,\theta }^{-1}y$, you just need some formulas from Multi Variable Derivatives. Now let's take some time to analyze what we good after all of this derivation: The quadratic form is interpretable as the likelihood, or the goodness of the fit. While the log of the Kernel matrix is the volume of the possible prediction, which is kinda the prior. Minimizing the above loss with gradient methods or similars can give us the best kernel parameters.

We can also obtain a Map Estimate if we place a prior on the kernel parameters $p(\theta )$.

Hyper-priors

Another way to optimize is using hyper-priors. Basically, we assume the priors have parameters too, and we try to optimize over those parameters. Usually a good way to choose kernels is posterior agreement. Given the data, we would like to make the probability of the given data as high as possible. This is very important but I have not understood this well

Approximating GPs

Leaning a Gaussian Process is usually quite slow. We will introduce some local methods for faster approximation, and a technique introduced in Rahimi, Recht NeurIPs 2007, it is a general algorithm for simplifying stationary kernels.

Local methods

One way is only use points such that satisfy $|k(x,x^{\prime })|\ge \tau$, the idea is that point that are far away usually have not such a great influence. This is still expensive if there are too many points close. Basically we can restrict ourselves to a subset of the Covariance matrix, as most other points will be set to $0$. This method essentially cuts off the tails of the distribution.

Kernel function approximation

This technique is quite hard to understand well. The idea is to reduce GPs to linear regression, then the cost would be $\mathcal{O}(n{m}^2+m^3)$ instead of $n^3$. One simple way to do it, is estimate the kernel function as some

$$k(x,x^{\prime })\approx \varphi (x{)}^T\varphi (x^{\prime })$$

Such that $\varphi (x)\in {\mathbb{R}}^m$, and $m$ is a much smaller dimensionality (easier to compute). Then it's just GP with a linear kernel, or a Bayesian linear regression.

Fourier Features

We make use of Bochner's theorem which states that a Kernel is a valid Mercer Kernel if its Fourier transform $p(\omega )$ is non-negative (which implies it is a valid probability distribution). When $p$ is a valid function, it is called spectral density.

Consider a shift invariant (stationary kernel), as they can be seen as function in one variable, these kernels have a Fourier transform, see Fourier Series.

$$k(x,x^{\prime })=k(x-x^{\prime })={\int }_{{\mathbb{R}}^d}p(\omega )e^{j{\omega }^T(x-x^{\prime })}d\omega ={\mathbb{E}}_{\omega ,b}[z_{\omega ,b}(x)\cdot z_{\omega ,b}(x^{\prime })]$$

Where $\omega \sim p(\omega ),b\sim U([0,2\pi ])$, and define $z_{\omega ,b}(x)=\sqrt{2}\cos ({\omega }^{T}x+b)$ where we are reinterpreting the integral as an expectation, one can show that these two are equivalent. The important thing to note here is that we have an unbiased approximator for some Kernel, given we have its spectral density.

Then we can approximate this expectation as follows:

$${\mathbb{E}}_{\omega ,b}[z_{\omega ,b}(x)\cdot z_{\omega ,b}(x^{\prime })]=\frac{1}{m}\sum _{i=1}^m{z}_{\omega ,b}(x)\cdot z_{\omega ,b}(x^{\prime })=z(x{)}^{T}z(x)$$

Where $z(x)=\frac{1}{\sqrt{m}}[z_{{\omega }_{1,b_1}}(x),\dots ,z_{{\omega }_{m,b_m}}(x)]$ Now we have an explicit random feature map, where $m$ works as a parameter for quality of the approximation.

One can also have some error analysis, telling you how close is the approximation to the original kernel. The theoretical thing is that this approximates the kernel function uniformly well, we use uniform convergence in this setting, which we don't know exactly what this is.

One nice thing to see is that one can prove BLR with fourier features is equivalent to GP with a stationary kernel.

Uniform Convergence of Fourier Features

The theorem states that if we have a compact subset of ${\mathbb{R}}^d$, called $\mathcal{M}$ of diameter $D_{\mathcal{M}}$, then for a stationary kernel $k$ and Fourier features $z$, it holds for any $\epsilon$ that:

$$\begin{array}{r}\mathbb{P}(|k(x,x^{\prime })-z(x{)}^{T}z(x^{\prime })|>\epsilon )\\ \le 2^8{\left(\frac{{\mathbb{E}}_{w\sim p}[w^{T}w]D_{\mathcal{M}}}{\epsilon }\right)}^2\exp \left(-\frac{m{\epsilon }^2}{8(d+2)}\right)\end{array}$$

Where we observe that the decay of the error is exponential in the number of Fourier features $m$

Proof: We first observe that $z(x{)}^{T}z(x)$ is shift invariant as the function that it attempts to approximate is. Then you can define some function $f(x,x^{\prime })=k(x,x^{\prime })-z(x{)}^{T}z(x^{\prime })$ and then we can define the function solely using the difference $\Delta =x-x^{\prime }$.

We use Hoeffding's Inequality to prove that:

$$\mathbb{P}(|f(\Delta )|>\epsilon )\le \exp \left(-m{\epsilon }^2\right)$$

Where we note that $0\le z(x{)}^{T}z(x)\le 2$

TODO: finish, this proof is hard...

Inducing points

Introduced in (Quinonero-Candela and Rasmussen, 2005), we want to summarize the points of the training set into the inducing points, and then use these points for inference.

Then the original Gaussian can be recovered in this manner:

$$p(f^{\ast },f)=\int p(f^{\ast },f,u)du=\int p(f^{\ast }\mid u)p(f\mid u)p(u)du$$

Recall that $f^{\ast }$ and $f$ are conditionally independen. The values $p(f^{\ast }\mid u)$ and $p(f\mid u)$ are respectively called testing conditional and training conditional.

Thus, the time complexity is cubic in the number of inducing points, but only linear in the number of data points.

Because we need to invert the matrix for the inducing point, but we don't need to invert the matrix for the training points.