## 2 Gaussian Processes  
### Gaussian Process Regression  
By definition, a Gaussian process (GP) is a continuous collection of random variables, any finite subset of which is normally distributed as a multivariate distribution. A more comprehensive treatment can be found in [1]. Here, we give a brief introduction with focus on practical implications. Let \(\mathcal{D}=\{X,\mathbf{y}\}\) be a dataset of \(n\) observations, where \(X=\left[\mathbf{x}_{1},\mathbf{x}_{2},...,\mathbf{x}_{n}\right]^{T}\) is an \(n\times d\) matrix of \(n\) input vectors of dimension \(d\) and \(\mathbf{y}=\left[y_{1},y_{2},...,y_{n}\right]^{T}\) is a vector of continuous-valued scalar outputs. These data points are also called training points. Regression aims to find a mapping \(f:\mathbb{R}^{d}\rightarrow\mathbb{R}\),  
\[y=f(\mathbf{x})+\epsilon,\quad\epsilon\sim\mathcal{N}(0,\sigma_{\epsilon}^{2}) \tag{1}\]  
with \(\epsilon\) being identically distributed observation noise. In stochastic numerics, this mapping is often implemented by a Gaussian process, which is fully defined by its mean \(m(\mathbf{x})\) and _covariance function_\(k(\mathbf{x},\mathbf{x}^{\prime})\). It is generally denoted as \(f(\mathbf{x})\sim\mathcal{GP}(m(\mathbf{x}),k(\mathbf{x},\mathbf{x}^{\prime}))\). Typically, the covariance function, also referred to as a _kernel_ or _kernel function_, is parametrized by a vector of hyperparameters \(\mathbf{\theta}\). By definition, a GP yields a distribution over a collection of functions that have a joint normal distribution [1],  
\[ \begin{bmatrix}\mathbf{f}\\ \mathbf{f}_{*}\end{bmatrix} \sim\mathcal{N}\left( \begin{bmatrix}\mathbf{m}_{X}\\ \mathbf{m}_{X_{*}}\end{bmatrix} , \begin{bmatrix}K_{X,X}&K_{X,X_{*}}\\ K_{X_{*},X}&K_{X_{*},X_{*}}\end{bmatrix} \right) \tag{2}\]  
where \(X\) and \(X_{*}\) are the input vectors of the \(n\) observed training points and \(n_{*}\) the unobserved test points, respectively. The vectors of mean values for \(X\) and \(X_{*}\) are given by \(\mathbf{m}_{X}\) and \(\mathbf{m}_{X_{*}}\). The covariance matrices \(K_{X,X}\), \(K_{X_{*},X_{*}}\), \(K_{X_{*},X}\) and \(K_{X,X_{*}}\) are constructed by evaluating \(k\) at their respective pairs of points. In practice, we don't have access to the latent function values directly, but are dependent on the noisy observations \(\mathbf{y}\). Putting it all together, the conditional (predictive posterior) distribution of the GP can be written as:  
\[\mathbf{f}_{*}|X,X_{*},\mathbf{y},\mathbf{\theta},\sigma_{\epsilon}^{2}\sim \mathcal{N}\left(\mathbb{E}(\mathbf{f}_{*}),\mathbb{V}(\mathbf{f}_{*})\right) \tag{3}\]  
\[\mathbb{E}(\mathbf{f}_{*})=\mathbf{m}_{X_{*}}+K_{X_{*},X}\left[K_{X,X}+\sigma_{ \epsilon}^{2}I\right]^{-1}\mathbf{f} \tag{4}\]  
\[\mathbb{V}(\mathbf{f}_{*})=K_{X_{*},X_{*}}-K_{X_{*},X}\left[K_{X,X}+\sigma_{ \epsilon}^{2}I\right]^{-1}K_{X,X_{*}} \tag{5}\]  
Without loss of generality, the mean function \(m\) is often set to zero. This does not imply a zero mean posterior \(\mathbb{E}(\mathbf{f}_{*})\), which in fact does not depend on the mean of the training points at all. However, recent work by [31] shows that the rate of convergence for Bayesian optimization (BO) can depend sensitively on the choice for a non-zero mean function. An extensive overview of kernels and how to implement and combine them can be found in [32]. Hyperparameters \(\mathbf{\theta}\) are usually learned by using a gradient based optimization algorithm to maximize the log marginal likelihood,  
\[\log p(\mathbf{y}|\theta,X)\propto-\frac{1}{2}\left[\mathbf{y}^{T}\left[K_{X,X}+\sigma_{\epsilon}^{2}I\right]\mathbf{y}+\log|K_{X,X}+\sigma_{\epsilon}^{2 }I)|\right] \tag{6}\]which is a combination of a data fit term and complexity penalty and thus automatically incorporates Occam's Razor [33]. This guards the GP against overfitting. A comparison of common optimizers can be found in [34]. In our experiments, we use BFGS, a quasi-Newton method described in [35].