## 2 Theoretical Background  
### Gaussian Process Regression  
A Gaussian process regression can be used to identify input-output relationships among observed data. The aim is to obtain a function \(f(\cdot)\) that approximates the given data points fully characterized by the mean and variance.  
Consider the dataset \(\mathcal{D}=\left\{\mathbf{z}_{k},y_{k}\right\}_{k=1}^{N}\) containing \(N\) pairs of input-output noisy data. The regression problem is of the form [3]:  
\[y_{k}=f(\mathbf{z}_{k})+e_{k},\quad e_{k}\sim\mathcal{N}(0,\sigma_{n}^{2}) \tag{1}\]  
where  
\[f(\mathbf{z})\sim\mathcal{GP}\big{(}m(\mathbf{z}),k(\mathbf{z},\mathbf{z}^{ \prime})\big{)} \tag{2}\]is a stochastic process with mean \(m(\mathbf{z})\) and covariance \(k(\mathbf{z},\mathbf{z}^{\prime})\):  
\[m(\mathbf{z}) =\mathbb{E}[f(\mathbf{z})] \tag{3}\] \[k(\mathbf{z},\mathbf{z}^{\prime}) =\mathbb{E}[(f(\mathbf{z})-m(\mathbf{z}))\big{(}f(\mathbf{z}^{ \prime})-m(\mathbf{z}^{\prime})\big{)}] \tag{4}\]  
corrupted by Gaussian noises with variance \(\sigma_{n}^{2}\).  
The mean and covariance functions are selected a priori. Usually, a function \(m(\mathbf{z})=0\) is considered [5]. For Equation (4), also called kernel, a squared exponential covariance function is usually employed:  
\[k(\mathbf{z},\mathbf{z}^{\prime})=\sigma_{f}^{2}\text{exp}\bigg{(}-\frac{\| \mathbf{z}-\mathbf{z}^{\prime}\|^{2}}{2l^{2}}\bigg{)}+\sigma_{n}^{2}\mathbb{I} \tag{5}\]  
with \(\sigma_{f}^{2}\) and \(l\) defining the magnitude and characteristic length scales of the regressor functions, respectively. Together with the variance \(\sigma_{n}^{2}\), they are known as hyperparameters of the GP. Let \(\mathbf{Z}=\{\mathbf{z}_{1},\dots,\mathbf{z}_{N}\}\) be a training dataset, with \(\mathbf{z}\in\mathbb{R}^{d}\), \(i,j=1,\dots,N\), and \(\mathbf{Z}_{*}=\{\mathbf{z}_{*1},\dots,\mathbf{z}_{*M}\}\) a test dataset. The goal is to find the predictive distribution of the test dataset based upon the set of \(N\) training input-output data pairs. Under the assumption of independence between noises and function values, the joint distribution of the observed values and the unknown function values can be written as  
\[\begin{pmatrix}\mathbf{y}\\ f(\mathbf{Z}_{*})\end{pmatrix}\sim\mathcal{N}\bigg{(}\begin{pmatrix}m(\mathbf{ Z})\\ m(\mathbf{Z}_{*})\end{pmatrix},\begin{pmatrix}(\mathbf{K}+\sigma_{n}^{2}\mathbb{I})& \mathbf{k}_{*}^{T}\\ \mathbf{k}_{*}&\mathbf{k}_{**}\end{pmatrix}\bigg{)} \tag{6}\]  
where \(\mathbf{y}\) is the vector of training targets, \(\mathbf{K}\) is the covariance matrix denoted by elements \(\mathbf{K}_{ij}=k(\mathbf{z}_{i},\mathbf{z}_{j})\), \(\mathbf{k}_{*}=k(\mathbf{Z}_{*},\mathbf{Z})\) is the matrix of covariances between training data and test data, and \(\mathbf{k}_{**}=k(\mathbf{Z}_{*},\mathbf{Z}_{*})\) is the autocovariance matrix of the test data.  
The posterior conditional distribution of the function values can be obtained by conditioning this joint Gaussian distribution on the observations \(\mathbf{y}\):  
\[f(\mathbf{Z}_{*})|\mathbf{Z},\mathbf{y},\mathbf{Z}_{*}\sim\mathcal{N}\big{(} \mathbf{\tilde{f}}_{*},\text{cov}(\mathbf{f}_{*})\big{)} \tag{7}\]  
where  
\[\mathbf{\tilde{f}}_{*} =\mathbf{k}_{*}[\mathbf{K}+\sigma_{n}^{2}\mathbb{I}]^{-1}\mathbf{ y} \tag{8}\] \[\text{Cov}(\mathbf{f}_{*}) =\mathbf{k}_{**}-\mathbf{k}_{*}[\mathbf{K}+\sigma_{n}^{2}\mathbb{ I}]^{-1}\mathbf{k}_{*}^{T} \tag{9}\]  
are the mean and the covariance of the normal distribution, respectively.