Exact GPR Method

An instance of response y from a Gaussian process regression (GPR) model can be modeled as

$$P\left(y_i|f\left(x_i\right),x_i\right) ~N\left(y_i|h{\left(x_i\right)}^T\beta +f\left(x_i\right),{\sigma }^2\right)$$

Hence, making predictions for new data from a GPR model requires:

• Knowledge of the coefficient vector, $$\beta$$, of fixed basis functions

• Ability to evaluate the covariance function $$k\left(x,x^{\prime }|\theta \right)$$ for arbitrary $$x$$ and $$x^{\prime }$$, given the kernel parameters or hyperparameters, $$\theta$$.

• Knowledge of the noise variance $${\sigma }^2$$ that appears in the density $$P\left(y_i|f\left(x_i\right),x_i\right)$$

That is, one needs first to estimate $$\beta$$, $$\theta$$, and $${\sigma }^2$$ from the data $$\left(X,y\right)$$.

Parameter Estimation

One approach for estimating the parameters $$\beta$$, $$\theta$$, and $${\sigma }^2$$ of a GPR model is by maximizing the likelihood $$P\left(y|X\right)$$ as a function of $$\beta$$, $$\theta$$, and $${\sigma }^2$$[1]. That is, if $$\widehat{\beta }$$, $$\widehat{\theta }$$, and $${\widehat{\sigma }}^2$$ are the estimates of $$\beta$$, $$\theta$$, and $${\sigma }^2$$, respectively, then:

$$\widehat{\beta },\widehat{\theta },{\widehat{\sigma }}^2=\underset{\beta ,\theta ,{\sigma }^2}{\text{arg max}}\log P\left(y|X,\beta ,\theta ,{\sigma }^2\right).$$

Because

$$P\left(y|X\right)=P\left(y|X,\beta ,\theta ,{\sigma }^2\right)=\mathcal{N}\left(y|H\beta ,K\left(X,X|\theta \right)+{\sigma }^2{I}_n\right),$$

the marginal log likelihood function is as follows:

$$\begin{array}{ll}\log P\left(y|X,\beta ,\theta ,{\sigma }^2\right)=\hfill & -\frac{1}{2}{\left(y-H\beta \right)}^T{\left[K\left(X,X|\theta \right)+{\sigma }^2{I}_n\right]}^{-1}\left(y-H\beta \right)\hfill \\ \hfill & -\frac{n}{2}\log 2\pi -\frac{1}{2}\log \left|K\left(X,X|\theta \right)+{\sigma }^2{I}_n\right|.\hfill \end{array}$$

where $$H$$ is the vector of explicit basis functions, and $$K\left(X,X|\theta \right)$$ is the covariance function matrix (for more information, see Gaussian Process Regression Models).

To estimate the parameters, the software first computes $$\widehat{\beta }\left(\theta ,{\sigma }^2\right)$$, which maximizes the log likelihood function with respect to $$\beta$$ for given $$\theta$$ and $${\sigma }^2$$. It then uses this estimate to compute the $$\beta$$-profiled likelihood:

$$\log \left\{P\left(y|X,\widehat{\beta }\left(\theta ,{\sigma }^2\right),\theta ,{\sigma }^2\right)\right\}.$$

The estimate of $$\beta$$ for given $$\theta$$, and $${\sigma }^2$$ is

$$\widehat{\beta }\left(\theta ,{\sigma }^2\right)={\left[H^T{\left[K\left(X,X|\theta \right)+{\sigma }^2{I}_n\right]}^{-1}H\right]}^{-1}{H}^T{\left[K\left(X,X|\theta \right)+{\sigma }^2{I}_n\right]}^{-1}y.$$

Then, the $$\beta$$-profiled log likelihood is given by

$$\begin{array}{ll}\log P\left(y|X,\widehat{\beta }(\theta ,{\sigma }^2),\theta ,{\sigma }^2\right)=\hfill & -\frac{1}{2}{\left(y-H\widehat{\beta }\left(\theta ,{\sigma }^2\right)\right)}^T{\left[K\left(X,X|\theta \right)+{\sigma }^2{I}_n\right]}^{-1}\left(y-H\widehat{\beta }\left(\theta ,{\sigma }^2\right)\right)\hfill \\ \hfill & -\frac{n}{2}\log 2\pi -\frac{1}{2}\log \left|K\left(X,X|\theta \right)+{\sigma }^2{I}_n\right|\hfill \end{array}$$

The software then maximizes the $$\beta$$-profiled log-likelihood over $$\theta$$, and $${\sigma }^2$$ to find their estimates.

Prediction

Making probabilistic predictions from a GPR model with known parameters requires the density $$P\left(y_{new}|y,X,x_{new}\right)$$. Using the definition of conditional probabilities, one can write:

$$P\left(y_{new}|y,X,x_{new}\right)=\frac{P\left(y_{new},y|X,x_{new}\right)}{P\left(y|X,x_{new}\right)}.$$

To find the joint density in the numerator, it is necessary to introduce the latent variables $$f_{new}$$ and $$f$$ corresponding to $$y_{new}$$, and $$y$$, respectively. Then, it is possible to use the joint distribution for $$y_{new}$$, $$y$$, $$f_{new}$$, and $$f$$ to compute $$P\left(y_{new},y|X,x_{new}\right)$$:

$$\begin{array}{l}\begin{array}{ll}P\left(y_{new},y|X,x_{new}\right)\hfill & ={\displaystyle \int {\displaystyle \int P}}\left(y_{new},y,f_{new},f|X,x_{new}\right)dfd{f}_{new}\hfill \\ \hfill & ={\displaystyle \int {\displaystyle \int P}}\left(y_{new},y|f_{new},f,X,x_{new}\right)P\left(f_{new},f|X,x_{new}\right)dfd{f}_{new}.\hfill \end{array}\\ \end{array}$$

Gaussian process models assume that each response $$y_i$$ only depends on the corresponding latent variable $$f_i$$ and the feature vector $$x_i$$. Writing $$P\left(y_{new},y|f_{new},f,X,x_{new}\right)$$ as a product of conditional densities and based on this assumption produces:

$$\begin{array}{l}P\left(y_{new},y|f_{new},f,X,x_{new}\right)=P\left(y_{new}|f_{new},x_{new}\right){\displaystyle \prod _{i=1}^{n}P}\left(y_i|f\left(x_i\right),x_i\right)\hfill \end{array}.$$

After integrating with respect to $$y_{new}$$, the result only depends on $$f$$ and $$X$$:

$$\begin{array}{l}P\left(y|f,X\right)={\displaystyle \prod _{i=1}^{n}P}\left(y_i|f_i,x_i\right)={\displaystyle \prod _{i=1}^n\mathcal{N}}\left(y_i{|h\left(x_i\right)}^T\beta +f_i,{\sigma }^2\right)\hfill \end{array}.$$

Hence,

$$P\left(y_{new},y|f_{new},f,X,x_{new}\right)=P\left(y_{new}|f_{new},x_{new}\right)P\left(y|f,X\right).$$

Again using the definition of conditional probabilities,

$$P\left(f_{new},f|X,x_{new}\right)=P\left(f_{new}|f,X,x_{new}\right)*P\left(f|X,x_{new}\right),$$

it is possible to write $$P\left(y_{new},y|X,x_{new}\right)$$ as follows:

$$P\left(y_{new},y|X,x_{new}\right)={\displaystyle \int {\displaystyle \int P}}\left(y_{new}|f_{new},x_{new}\right)P\left(y|f,X\right)P\left(f_{new}|f,X,x_{new}\right)P\left(f|X,x_{new}\right)dfd{f}_{new}.$$

Using the facts that

$$P\left(f|X,x_{new}\right)=P\left(f|X\right)$$

and

$$P\left(y|f,X\right)P\left(f|X\right)=P\left(y,f|X\right)=P\left(f|y,X\right)P\left(y|X\right),$$

one can rewrite $$P\left(y_{new},y|X,x_{new}\right)$$ as follows:

$$P\left(y_{new},y|X,x_{new}\right)=P\left(y|X\right){\displaystyle \int {\displaystyle \int P}}\left(y_{new}|f_{new},x_{new}\right)P\left(f|y,X\right)P\left(f_{new}|f,X,x_{new}\right)dfd{f}_{new}.$$

It is also possible to show that

$$P\left(y|X,x_{new}\right)=P\left(y|X\right).$$

Hence, the required density $$P\left(y_{new}|y,X,x_{new}\right)$$ is:

$$\begin{array}{ll}P\left(y_{new}|y,X,x_{new}\right)\hfill & =\frac{P\left(y_{new},y|X,x_{new}\right)}{P\left(y|X,x_{new}\right)}=\frac{P\left(y_{new},y|X,x_{new}\right)}{P\left(y|X\right)}\hfill \\ \hfill & ={\displaystyle \int {\displaystyle \int \underset{\left(1\right)}{\underbrace{P\left(y_{new}|f_{new},x_{new}\right)}}\underset{\left(2\right)}{\underbrace{P\left(f|y,X\right)}}\underset{\left(3\right)}{\underbrace{P\left(f_{new}|f,X,x_{new}\right)}}dfd{f}_{new}.}}\hfill \end{array}$$

It can be shown that

$$\left(1\right)P\left(y_{new}|f_{new},x_{new}\right)=\mathcal{N}\left(y_{new}|h{(x_{new})}^T\beta +f_{new},{\sigma }_{new}^2\right)$$

$$\left(2\right)P\left(f|y,X\right)=\mathcal{N}\left(f|\frac{1}{{\sigma }^2}{\left(\frac{I_n}{{\sigma }^2}+K{\left(X,X\right)}^{-1}\right)}^{-1}\left(y-H\beta \right),{\left(\frac{I_n}{{\sigma }^2}+K{(X,X)}^{-1}\right)}^{-1}\right)$$

$$\begin{array}{l}\left(3\right)\begin{array}{l}P\left(f_{new}|f,X,x_{new}\right)=\mathcal{N}\left(f_{new}|K\left(x_{new}^T,X\right)K{\left(X,X\right)}^{-1}f,\Delta \right)\hfill \end{array},\\ \text{where}\Delta =k\left(x_{new},x_{new}\right)-K\left(x_{new}^T,X\right)K{\left(X,X\right)}^{-1}K\left(X,x_{new}^T\right).\end{array}$$

After the integration and required algebra, the density of the new response $$y_{new}$$ at a new point $$x_{new}$$, given $$y$$, $$X$$ is found as

$$P\left(y_{new}|y,X,x_{new}\right)=\mathcal{N}\left(y_{new}|h{(x_{new})}^T\beta +\mu ,{\sigma }_{new}^2+\Sigma \right),$$

where

$$\mu =K\left(x_{new}^T,X\right)\underset{\alpha }{\underbrace{{\left(K(X,X)+{\sigma }^2{I}_n\right)}^{-1}\left(y-H\beta \right)}}$$

and

$$\Sigma =k\left(x_{new},x_{new}\right)-K\left(x_{new}^T,X\right){\left(K\left(X,X\right)+{\sigma }^2{I}_n\right)}^{-1}K\left(X,x_{new}^T\right).$$

The expected value of prediction $$y_{new}$$ at a new point $$x_{new}$$ given $$y$$, $$X$$, and parameters $$\beta$$, $$\theta$$, and $${\sigma }^2$$ is

$$\begin{array}{ll}E\left(y_{new}|y,X,x_{new},\beta ,\theta ,{\sigma }^2\right)\hfill & =h{\left(x_{new}\right)}^T\beta +K\left(x_{new}^T,X|\theta \right)\alpha \hfill \\ \hfill & =h{\left(x_{new}\right)}^T\beta +{\displaystyle \sum _{i=1}^n{\alpha }_i}k\left(x_{new},x_i|\theta \right),\hfill \end{array}$$

where

$$\alpha ={\left(K\left(X,X|\theta \right)+{\sigma }^2{I}_n\right)}^{-1}\left(y-H\beta \right).$$

Computational Complexity of Exact Parameter Estimation and Prediction

Training a GPR model with the exact method (when FitMethod is 'Exact') requires the inversion of an n-by-n kernel matrix $$K\left(X,X\right)$$. The memory requirement for this step scales as O(n²) since $$K\left(X,X\right)$$ must be stored in memory. One evaluation of $$\log P\left(y|X\right)$$ scales as O(n³). Therefore, the computational complexity is O(kn³), where k is the number of function evaluations needed for maximization and n is the number of observations.

Making predictions on new data involves the computation of $$\widehat{\alpha }$$. If prediction intervals are desired, this step could also involve the computation and storage of the Cholesky factor of $$\left(K\left(X,X\right)+{\sigma }^2{I}_n\right)$$ for later use. The computational complexity of this step using the direct computation of $$\widehat{\alpha }$$ is O(n³) and the memory requirement is O(n²).

Hence, for large n, estimation of parameters or computing predictions can be very expensive. The approximation methods usually involve rearranging the computation so as to avoid the inversion of an n-by-n matrix. For the available approximation methods, please see the related links at the bottom of the page.

References

[1] Rasmussen, C. E. and C. K. I. Williams. Gaussian Processes for Machine Learning. MIT Press. Cambridge, Massachusetts, 2006.