# Gaussian Processes for Signal Processing and Representation in Control Engineering
## 3 Computational Methods of Gaussian Process
### Approximation of Gaussian Process  
GP computation is problematic, because of the earlier mentioned necessity of matrix inversion, which grows in size with the size of the training dataset. Nowadays, computers can easily handle reasonable sized problems with GP as solution. On the other hand, if we are dealing with large GP model and its computation becomes a burden, there are several methods of approximating GP. An example of such solution are sparse methods. Common to all different variations is that only a subset of the latent variables is treated exactly, and the remaining variables are given some approximate, which is computationally easier [16]. The authors of [16; 17] had described the whole idea in detail, but in general their approach is interpretable as exact inference with an approximate prior. To be more specific, if we infer from the GP model by putting the joint GP prior on training \(f\) and testing latent values \(f_{*}\) and combining it with the likelihood \(p(y|f)\) we achieve a joint posterior given by Equation (16). \[p(\boldsymbol{f}_{*}\mid\mathbf{y})=\frac{p(\boldsymbol{f}_{*}\boldsymbol{f}_{* })p(\mathbf{y}\mid\boldsymbol{f})}{p(\mathbf{y})}. \tag{16}\]  
Then there is a need to marginalize out the unwanted training set latent variables:  
\[p(\boldsymbol{f}_{*}\mid\mathbf{y})=\int p(\boldsymbol{f}_{*}\mid\mathbf{y}) \mathrm{d}\boldsymbol{f}=\frac{1}{p(\mathbf{y})}\int p(\mathbf{y}\mid \boldsymbol{f})p(\boldsymbol{f}_{*}\boldsymbol{f}_{*})\mathrm{d}\boldsymbol{f}. \tag{17}\]  
The posterior predictive is received:  
\[p(\boldsymbol{f}_{*}\boldsymbol{f}_{*})=\mathcal{N}\bigg{(}\mathbf{0},\left[ \begin{array}{cc}K_{\boldsymbol{f},\boldsymbol{f}}&K_{\boldsymbol{*}, \boldsymbol{f}}\\ K_{\boldsymbol{f},*}&K_{\boldsymbol{*},*}\end{array} \right]\bigg{)},\quad \text{and}\quad p(\mathbf{y}\mid\boldsymbol{f})=\mathcal{N}\Big{(}f,\sigma _{\text{noise}}^{2}\;I\Big{)}. \tag{18}\]Since both factors in the integral are Gaussian, the integral can be evaluated in closed form to give the Gaussian predictive distribution:  
\[p(\mathbf{f}_{*}\mid\mathbf{y})=\mathcal{N}\bigg{(}K_{*,\mathbf{f}}\Big{(}K_{\mathbf{f},\mathbf{ f}}+\sigma_{\text{noise}}^{2}\;I\Big{)}^{-1}\mathbf{y},K_{*,*}-K_{*,\mathbf{f}}\Big{(}K_{ \mathbf{f},\mathbf{f}}+\sigma_{\text{noise}}^{2}\;I\Big{)}^{-1}K_{\mathbf{f},*}\bigg{)}. \tag{19}\]  
The problem is that the equation requires the earlier mentioned inversion of matrix, the size of which depends on the number of training cases [16; 17]. The idea is to modify joint prior \(p(\mathbf{f}_{*},\mathbf{f})\) to reduce its computation requirements. In order to achieve this, inducing variables \(u\), which are a set of latent variables that are values of GP, corresponding to a set of input locations (which are called the inducing inputs). The consistency of GP allows us to recover the mentioned joint prior to \(p(\mathbf{f}_{*},\mathbf{f})\) by integrating (marginalizing) out \(\mathfrak{u}\) from joint GP prior \(p(\mathbf{f}_{*},\mathbf{f},u)\). \[p(\mathbf{f}_{*},\mathbf{f})=\int p(\mathbf{f}_{*},\mathbf{f},\mathbf{u})\mathrm{d}\mathbf{u}=\int p( \mathbf{f}_{*},\mathbf{f}\mid\mathbf{u})p(\mathbf{u})\mathrm{d}\mathbf{u},\quad\text{ where }\quad p(\mathbf{u})=\mathcal{N}(\mathbf{0},K_{\mathbf{u},\mathbf{u}}) \tag{20}\]  
Then joint prior is approximated by assuming that \(f_{*}\) and \(f\) are conditionally independent given \(u\). \[p(\mathbf{f}_{*},\mathbf{f})\simeq q(\mathbf{f}_{*},\mathbf{f})=\int q(\mathbf{f}_{*}\mid\mathbf{u}) q(\mathbf{f}\mid\mathbf{u})p(\mathbf{u})\mathrm{d}\mathbf{u}. \tag{21}\]  
The way of choosing the inducing variables as well as making addition assumptions to Equation (21) depends on an exact chosen sparse algorithm [16; 17]. Another approach to the approximation of GP is presented in [9] as Hilbert space approximate GP model for a low-rank representation of stationary GPs. Main principal consist of using basis function approximations based on approximation via Laplace eigenfunctions for stationary covariance functions proposed in [18].