Usual covariance functions

The effect of choosing different kernels on the prior function distribution of the Gaussian process. Left is a squared exponential kernel. Middle is Brownian. Right is quadratic.

There are a number of common covariance functions:^[7]

• Constant : ${\displaystyle K_{\operatorname {C} }(x,x')=C}$
• Linear: ${\displaystyle K_{\operatorname {L} }(x,x')=x^{\mathsf {T}}x'}$
• white Gaussian noise: ${\displaystyle K_{\operatorname {GN} }(x,x')=\sigma ^{2}\delta _{x,x'}}$
• Squared exponential: ${\displaystyle K_{\operatorname {SE} }(x,x')=\exp \left(-{\tfrac {|d|^{2}}{2\ell ^{2}}}\right)}$
• Ornstein–Uhlenbeck: ${\displaystyle K_{\operatorname {OU} }(x,x')=\exp \left(-{\tfrac {|d|}{\ell }}\right)}$
• Matérn: ${\displaystyle K_{\operatorname {Matern} }(x,x')={\tfrac {2^{1-\nu }}{\Gamma (\nu )}}\left({\tfrac {{\sqrt {2\nu }}|d|}{\ell }}\right)^{\nu }K_{\nu }\left({\tfrac {{\sqrt {2\nu }}|d|}{\ell }}\right)}$
• Periodic: ${\displaystyle K_{\operatorname {P} }(x,x')=\exp \left(-{\tfrac {2}{\ell ^{2}}}\sin ^{2}(d/2)\right)}$
• Rational quadratic: ${\displaystyle K_{\operatorname {RQ} }(x,x')=\left(1+|d|^{2}\right)^{-\alpha },\quad \alpha \geq 0}$

Here ${\displaystyle d=x-x'}$. The parameter ${\displaystyle \ell }$ is the characteristic length-scale of the process (practically, "how close" two points ${\displaystyle x}$ and ${\displaystyle x'}$ have to be to influence each other significantly), ${\displaystyle \delta }$ is the Kronecker delta and ${\displaystyle \sigma }$ the standard deviation of the noise fluctuations. Moreover, ${\displaystyle K_{\nu }}$ is the modified Bessel function of order ${\displaystyle \nu }$ and ${\displaystyle \Gamma (\nu )}$ is the gamma function evaluated at ${\displaystyle \nu }$. Importantly, a complicated covariance function can be defined as a linear combination of other simpler covariance functions in order to incorporate different insights about the data-set at hand.

The inferential results are dependent on the values of the hyperparameters ${\displaystyle \theta }$ (e.g. ${\displaystyle \ell }$ and ${\displaystyle \sigma }$) defining the model's behaviour. A popular choice for ${\displaystyle \theta }$ is to provide maximum a posteriori (MAP) estimates of it with some chosen prior. If the prior is very near uniform, this is the same as maximizing the marginal likelihood of the process; the marginalization being done over the observed process values ${\displaystyle y}$.^[7] This approach is also known as maximum likelihood II, evidence maximization, or empirical Bayes.^[9]

Stationary case

For a stationary Gaussian process ${\displaystyle X=(X_{t})_{t\in \mathbb {R} },}$ some conditions on its spectrum are sufficient for sample continuity, but fail to be necessary. A necessary and sufficient condition, sometimes called Dudley–Fernique theorem, involves the function ${\displaystyle \sigma }$ defined by

$${\displaystyle \sigma (h)={\sqrt {{\mathbb {E} }{\big [}X(t+h)-X(t){\big ]}^{2}}}}$$

(the right-hand side does not depend on ${\displaystyle t}$ due to stationarity). Continuity of ${\displaystyle X}$ in probability is equivalent to continuity of ${\displaystyle \sigma }$ at ${\displaystyle 0.}$ When convergence of ${\displaystyle \sigma (h)}$ to ${\displaystyle 0}$ (as ${\displaystyle h\to 0}$) is too slow, sample continuity of ${\displaystyle X}$ may fail. Convergence of the following integrals matters:

$${\displaystyle I(\sigma )=\int _{0}^{1}{\frac {\sigma (h)}{h{\sqrt {\log(1/h)}}}}\,dh=\int _{0}^{\infty }2\sigma (e^{-x^{2}})\,dx,}$$

these two integrals being equal according to integration by substitution ${\textstyle h=e^{-x^{2}},}$ ${\textstyle x={\sqrt {\log(1/h)}}.}$ The first integrand need not be bounded as ${\displaystyle h\to 0+,}$ thus the integral may converge (${\displaystyle I(\sigma )<\infty }$) or diverge (${\displaystyle I(\sigma )=\infty }$). Taking for example ${\textstyle \sigma (e^{-x^{2}})={\tfrac {1}{x^{a}}}}$ for large ${\displaystyle x,}$ that is, ${\textstyle \sigma (h)=(\log(1/h))^{-a/2}}$ for small ${\displaystyle h,}$ one obtains ${\displaystyle I(\sigma )<\infty }$ when ${\displaystyle a>1,}$ and ${\displaystyle I(\sigma )=\infty }$ when ${\displaystyle 0<a\leq 1.}$ In these two cases the function ${\displaystyle \sigma }$ is increasing on ${\displaystyle [0,\infty ),}$ but generally it is not. Moreover, the condition

(*)   there exists ${\displaystyle \varepsilon >0}$ such that ${\displaystyle \sigma }$ is monotone on ${\displaystyle [0,\varepsilon ]}$

does not follow from continuity of ${\displaystyle \sigma }$ and the evident relations ${\displaystyle \sigma (h)\geq 0}$ (for all ${\displaystyle h}$) and ${\displaystyle \sigma (0)=0.}$

Some history.^[17]: 424 Sufficiency was announced by Xavier Fernique in 1964, but the first proof was published by Richard M. Dudley in 1967.^{[16]: Theorem 7.1} Necessity was proved by Michael B. Marcus and Lawrence Shepp in 1970.^[18]: 380

There exist sample continuous processes ${\displaystyle X}$ such that ${\displaystyle I(\sigma )=\infty$ ;} they violate condition (*). An example found by Marcus and Shepp ^[18]: 387 is a random lacunary Fourier series

$${\displaystyle X_{t}=\sum _{n=1}^{\infty }c_{n}(\xi _{n}\cos \lambda _{n}t+\eta _{n}\sin \lambda _{n}t),}$$

where ${\displaystyle \xi _{1},\eta _{1},\xi _{2},\eta _{2},\dots }$ are independent random variables with standard normal distribution; frequencies ${\displaystyle 0<\lambda _{1}<\lambda _{2}<\dots }$ are a fast growing sequence; and coefficients ${\displaystyle c_{n}>0}$ satisfy ${\textstyle \sum _{n}c_{n}<\infty .}$ The latter relation implies

${\textstyle {\mathbb {E} }\sum _{n}c_{n}(|\xi _{n}|+|\eta _{n}|)=\sum _{n}c_{n}{\mathbb {E} }[|\xi _{n}|+|\eta _{n}|]={\text{const}}\cdot \sum _{n}c_{n}<\infty ,}$

whence ${\textstyle \sum _{n}c_{n}(|\xi _{n}|+|\eta _{n}|)<\infty }$ almost surely, which ensures uniform convergence of the Fourier series almost surely, and sample continuity of ${\displaystyle X.}$

Autocorrelation of a random lacunary Fourier series

Its autocovariation function

$${\displaystyle {\mathbb {E} }[X_{t}X_{t+h}]=\sum _{n=1}^{\infty }c_{n}^{2}\cos \lambda _{n}h}$$

is nowhere monotone (see the picture), as well as the corresponding function ${\displaystyle \sigma ,}$

$${\displaystyle \sigma (h)={\sqrt {2{\mathbb {E} }[X_{t}X_{t}]-2{\mathbb {E} }[X_{t}X_{t+h}]}}=2{\sqrt {\sum _{n=1}^{\infty }c_{n}^{2}\sin ^{2}{\frac {\lambda _{n}h}{2}}}}.}$$

Gaussian process prediction, or Kriging

Gaussian Process Regression (prediction) with a squared exponential kernel. Left plot are draws from the prior function distribution. Middle are draws from the posterior. Right is mean prediction with one standard deviation shaded.

When concerned with a general Gaussian process regression problem (Kriging), it is assumed that for a Gaussian process ${\displaystyle f}$ observed at coordinates ${\displaystyle x}$, the vector of values ${\displaystyle f(x)}$ is just one sample from a multivariate Gaussian distribution of dimension equal to number of observed coordinates ${\displaystyle n}$. Therefore, under the assumption of a zero-mean distribution, ${\displaystyle f(x')\sim N(0,K(\theta ,x,x'))}$, where ${\displaystyle K(\theta ,x,x')}$ is the covariance matrix between all possible pairs ${\displaystyle (x,x')}$ for a given set of hyperparameters θ.^[7] As such the log marginal likelihood is:

$${\displaystyle \log p(f(x')\mid \theta ,x)=-{\frac {1}{2}}\left(f(x)^{\mathsf {T}}K(\theta ,x,x')^{-1}f(x')+\log \det(K(\theta ,x,x'))+n\log 2\pi \right)}$$

and maximizing this marginal likelihood towards θ provides the complete specification of the Gaussian process f. One can briefly note at this point that the first term corresponds to a penalty term for a model's failure to fit observed values and the second term to a penalty term that increases proportionally to a model's complexity. Having specified θ, making predictions about unobserved values ${\displaystyle f(x^{*})}$ at coordinates x* is then only a matter of drawing samples from the predictive distribution ${\displaystyle p(y^{*}\mid x^{*},f(x),x)=N(y^{*}\mid A,B)}$ where the posterior mean estimate A is defined as

$${\displaystyle A=K(\theta ,x^{*},x)K(\theta ,x,x')^{-1}f(x)}$$

and the posterior variance estimate B is defined as:

$${\displaystyle B=K(\theta ,x^{*},x^{*})-K(\theta ,x^{*},x)K(\theta ,x,x')^{-1}K(\theta ,x^{*},x)^{\mathsf {T}}}$$

where ${\displaystyle K(\theta ,x^{*},x)}$ is the covariance between the new coordinate of estimation x* and all other observed coordinates x for a given hyperparameter vector θ, ${\displaystyle K(\theta ,x,x')}$ and ${\displaystyle f(x)}$ are defined as before and ${\displaystyle K(\theta ,x^{*},x^{*})}$ is the variance at point x* as dictated by θ. It is important to note that practically the posterior mean estimate of ${\displaystyle f(x^{*})}$ (the "point estimate") is just a linear combination of the observations ${\displaystyle f(x)}$; in a similar manner the variance of ${\displaystyle f(x^{*})}$ is actually independent of the observations ${\displaystyle f(x)}$. A known bottleneck in Gaussian process prediction is that the computational complexity of inference and likelihood evaluation is cubic in the number of points |x|, and as such can become unfeasible for larger data sets.^[6] Works on sparse Gaussian processes, that usually are based on the idea of building a representative set for the given process f, try to circumvent this issue.^[29][30] The kriging method can be used in the latent level of a nonlinear mixed-effects model for a spatial functional prediction: this technique is called the latent kriging.^[31]

Often, the covariance has the form ${\textstyle K(\theta ,x,x')={\frac {1}{\sigma ^{2}}}{\tilde {K}}(\theta ,x,x')}$, where ${\displaystyle \sigma ^{2}}$ is a scaling parameter. Examples are the Matérn class covariance functions. If this scaling parameter ${\displaystyle \sigma ^{2}}$ is either known or unknown (i.e. must be marginalized), then the posterior probability, ${\displaystyle p(\theta \mid D)}$, i.e. the probability for the hyperparameters ${\displaystyle \theta }$ given a set of data pairs ${\displaystyle D}$ of observations of ${\displaystyle x}$ and ${\displaystyle f(x)}$, admits an analytical expression.^[32]

Bayesian neural networks as Gaussian processes

Bayesian neural networks are a particular type of Bayesian network that results from treating deep learning and artificial neural network models probabilistically, and assigning a prior distribution to their parameters. Computation in artificial neural networks is usually organized into sequential layers of artificial neurons. The number of neurons in a layer is called the layer width. As layer width grows large, many Bayesian neural networks reduce to a Gaussian process with a closed form compositional kernel. This Gaussian process is called the Neural Network Gaussian Process (NNGP).^[7][33][34] It allows predictions from Bayesian neural networks to be more efficiently evaluated, and provides an analytic tool to understand deep learning models.

Literature

• The Gaussian Processes Web Site, including the text of Rasmussen and Williams' Gaussian Processes for Machine Learning
• Ebden, Mark (2015). "Gaussian Processes: A Quick Introduction". arXiv:1505.02965 [math.ST].
• A Review of Gaussian Random Fields and Correlation Functions
• Efficient Reinforcement Learning using Gaussian Processes

Software

• GPML: A comprehensive Matlab toolbox for GP regression and classification
• STK: a Small (Matlab/Octave) Toolbox for Kriging and GP modeling
• Kriging module in UQLab framework (Matlab)
• CODES Toolbox: implementations of Kriging, variational kriging and multi-fidelity models (Matlab)
• Matlab/Octave function for stationary Gaussian fields
• Yelp MOE – A black box optimization engine using Gaussian process learning
• ooDACE – A flexible object-oriented Kriging Matlab toolbox.
• GPstuff – Gaussian process toolbox for Matlab and Octave
• GPy – A Gaussian processes framework in Python
• GSTools - A geostatistical toolbox, including Gaussian process regression, written in Python
• Interactive Gaussian process regression demo
• Basic Gaussian process library written in C++11
• scikit-learn – A machine learning library for Python which includes Gaussian process regression and classification
• - The Kriging toolKit (KriKit) is developed at the Institute of Bio- and Geosciences 1 (IBG-1) of Forschungszentrum Jülich (FZJ)

Video tutorials

• Gaussian Process Basics by David MacKay
• Learning with Gaussian Processes by Carl Edward Rasmussen
• Bayesian inference and Gaussian processes by Carl Edward Rasmussen