Gaussian process

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:[9]

• Constant : ${\displaystyle K_{\operatorname {C} }(x,x')=C}$
• Linear: ${\displaystyle K_{\operatorname {L} }(x,x')=x^{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 {\Big (}-{\frac {|d|^{2}}{2\ell ^{2}}}{\Big )}}$
• Ornstein–Uhlenbeck: ${\displaystyle K_{\operatorname {OU} }(x,x')=\exp \left(-{\frac {|d|}{\ell }}\right)}$
• Matérn: ${\displaystyle K_{\operatorname {Matern} }(x,x')={\frac {2^{1-\nu }}{\Gamma (\nu )}}{\Big (}{\frac {{\sqrt {2\nu }}|d|}{\ell }}{\Big )}^{\nu }K_{\nu }{\Big (}{\frac {{\sqrt {2\nu }}|d|}{\ell }}{\Big )}}$
• Periodic: ${\displaystyle K_{\operatorname {P} }(x,x')=\exp \left(-{\frac {2\sin ^{2}\left({\frac {d}{2}}\right)}{\ell ^{2}}}\right)}$
• Rational quadratic: ${\displaystyle K_{\operatorname {RQ} }(x,x')=(1+|d|^{2})^{-\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.

Clearly, 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}$.[9] This approach is also known as maximum likelihood II, evidence maximization, or empirical Bayes.[11]

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 (\mathbb {e} ^{-x^{2}})\,dx,}$

these two integrals being equal according to integration by substitution ${\displaystyle h=\mathbb {e} ^{-x^{2}},}$ ${\displaystyle \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 ${\displaystyle \sigma (\mathbb {e} ^{-x^{2}})={\tfrac {1}{x^{a}}}}$ for large ${\displaystyle x,}$ that is, ${\displaystyle \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

${\displaystyle (*)}$   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.}$

Theorem 1.   Let ${\displaystyle \sigma }$ be continuous and satisfy ${\displaystyle (*).}$ Then the condition ${\displaystyle I(\sigma )<\infty }$ is necessary and sufficient for sample continuity of ${\displaystyle X.}$

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

There exist sample continuous processes ${\displaystyle X}$ such that ${\displaystyle I(\sigma )=\infty$ ;} they violate condition ${\displaystyle (*).}$ An example found by Marcus and Shepp [20]^: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 ${\displaystyle \textstyle \sum _{n}c_{n}<\infty .}$ The latter relation implies ${\displaystyle \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 ${\displaystyle \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 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 θ.[9] As such the log marginal likelihood is:

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

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)^{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 ${\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 prediction is cubic in the number of points |x| and as such can become unfeasible for larger data sets.[8] 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.[28][29]

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). It allows predictions from Bayesian neural networks to be more efficiently evaluated, and provides an analytic tool to understand deep learning models.

• 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)
• 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)

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