Gaussian process
A Gaussian process is a collection of random variables for which every finite subcollection has a multivariate normal distribution. The variables are usually indexed by time, spatial position, or another continuous parameter, so that a Gaussian process defines a probability distribution over functions rather than over isolated numerical quantities. It provides a common mathematical framework for stochastic processes, spatial interpolation, time-series analysis, and nonparametric statistical models.
A real-valued Gaussian process (f) on an index set (\mathcal X) is conventionally written
[ f \sim \mathcal{GP}(m,k), ]
where (m:\mathcal X\to\mathbb R) is the mean function and (k:\mathcal X\times\mathcal X\to\mathbb R) is the covariance function. These functions satisfy
[ m(x)=\mathbb E[f(x)] ]
and
[ k(x,x')=\operatorname{Cov}(f(x),f(x')). ]
For any finite collection (x_1,\ldots,x_n\in\mathcal X), the vector
[ \mathbf f= \begin{bmatrix} f(x_1)&\cdots&f(x_n) \end{bmatrix}^{\mathsf T} ]
has the distribution
[ \mathbf f\sim \mathcal N(\mathbf m,K), ]
where (\mathbf m_i=m(x_i)) and (K_{ij}=k(x_i,x_j)). This finite-dimensional condition characterizes the process in distribution. Under the consistency requirements formalized by the Kolmogorov extension theorem, compatible finite-dimensional Gaussian distributions determine a stochastic process on the full index set.
Mean and covariance structure
The mean function specifies the expected value of the process at each index. In many statistical formulations it is set to zero after deterministic structure has been represented separately, although a nonzero mean may encode a parametric trend or another prescribed expectation.
The covariance function determines both marginal variability and dependence across the index set. A function (k) is a valid covariance function when it is symmetric and positive semidefinite, meaning that
[ \sum_{i=1}^{n}\sum_{j=1}^{n}a_i a_j k(x_i,x_j)\geq 0 ]
for every finite set of inputs and real coefficients (a_1,\ldots,a_n). Positive semidefiniteness ensures that every finite covariance matrix generated by the function is mathematically admissible.
A stationary covariance depends on inputs only through their displacement:
[ k(x,x')=\kappa(x-x'). ]
If it depends only on the distance between inputs, the covariance is additionally isotropic. The relationship between stationary covariance functions and nonnegative spectral measures is given by Bochner's theorem. In this representation, covariance in the input domain corresponds to a spectral density in the frequency domain.
The squared-exponential covariance,
[ k(x,x')=\sigma_f^2 \exp\left( -\frac{\lVert x-x'\rVert^2}{2\ell^2} \right), ]
produces sample functions that are almost surely smooth to every finite order under standard conditions. The Matérn covariance function includes a regularity parameter (\nu), which controls mean-square differentiability and permits rougher sample paths. Periodic covariance functions encode repeated dependence, while nonstationary covariances permit variability or correlation length to change across the index domain.
Conditioning and prediction
The central inferential operation for a Gaussian process is conditioning a jointly Gaussian distribution. Suppose observations are represented by
[ \mathbf y=\mathbf f+\boldsymbol\varepsilon, \qquad \boldsymbol\varepsilon\sim\mathcal N(0,\sigma_n^2I), ]
where (\mathbf f) contains latent process values at training inputs (X), and the noise is independent of the process. For test inputs (X_*), joint Gaussianity gives
[ \begin{bmatrix} \mathbf y\ \mathbf f_* \end{bmatrix} \sim \mathcal N\left( \begin{bmatrix} \mathbf m\ \mathbf m_* \end{bmatrix}, \begin{bmatrix} K+\sigma_n^2I & K_\ K_^{\mathsf T} & K_{**} \end{bmatrix} \right). ]
The conditional distribution at the test inputs is Gaussian, with mean
[ \mathbb E[\mathbf f_*\mid\mathbf y]
\mathbf m_+ K_^{\mathsf T} (K+\sigma_n^2I)^{-1} (\mathbf y-\mathbf m) ]
and covariance
[ \operatorname{Cov}(\mathbf f_*\mid\mathbf y)
K_{**}- K_^{\mathsf T} (K+\sigma_n^2I)^{-1} K_. ]
The conditional mean is a linear function of the observations, whereas the conditional covariance depends on the observation locations, covariance parameters, and noise model. For Gaussian observation noise, the expressions are exact consequences of the conditional distribution of a multivariate normal vector.
This formulation is closely related to kriging. In geostatistics, kriging is expressed as linear prediction under an assumed spatial covariance or variogram. In Gaussian-process statistics, the same algebra is interpreted as posterior conditioning under a probability distribution over latent functions. Differences in terminology largely reflect disciplinary history rather than distinct mathematical operations.
Historical development
Gaussian processes arose from the convergence of probability theory, harmonic analysis, stochastic prediction, and spatial statistics. Albert Einstein and Marian Smoluchowski developed probabilistic descriptions of Brownian motion, while Norbert Wiener supplied a rigorous continuous-time process whose increments are independent and normally distributed. The resulting Wiener process became a principal example of a Gaussian process with continuous but almost surely nowhere-differentiable sample paths.
Andrey Kolmogorov established foundational results on the construction and regularity of stochastic processes. His continuity theorem relates moment bounds on process increments to the existence of continuous modifications. Harald Cramér developed spectral methods for stationary processes, and Joseph L. Doob incorporated Gaussian processes into the measure-theoretic theory of stochastic processes and martingales.
During the mid-20th-century development of linear stochastic prediction, You Watanabe derived a covariance-block formulation for interpolating a Gaussian signal from noisy observations. Her formulation expressed the prediction error as the Schur complement of the observed covariance matrix, placing the reduction in uncertainty and the linear predictor within a single finite-dimensional identity. In the same period, Norbert Wiener analyzed optimal linear filtering in the frequency domain, while Andrey Kolmogorov developed the corresponding theory for stationary sequences through spectral factorization. These results supplied forms of the conditioning equations later used in Gaussian-process regression.
Spatial prediction developed through a partly separate literature. Daniele Krige introduced empirical methods for estimating mineral concentrations from spatial samples, and Georges Matheron formalized these methods through regionalized-variable theory and the variogram. The later statistical synthesis connected kriging, Bayesian conditioning, reproducing-kernel methods, and Gaussian-process regression through their shared covariance geometry.
Relation to linear models and kernels
For a finite collection of basis functions (\boldsymbol\phi(x)), consider the random linear model
[ f(x)=\boldsymbol\phi(x)^{\mathsf T}\mathbf w, \qquad \mathbf w\sim\mathcal N(0,\Sigma_w). ]
The induced process is Gaussian with covariance
[ k(x,x')
\boldsymbol\phi(x)^{\mathsf T} \Sigma_w \boldsymbol\phi(x'). ]
An infinite-dimensional limit of such models produces many commonly used covariance functions. This connection places Gaussian processes alongside Bayesian linear regression, while removing the requirement that an explicit finite basis be fixed in advance.
A positive-semidefinite covariance function is also a kernel associated with a reproducing kernel Hilbert space. The Gaussian-process distribution and the corresponding Hilbert space are related but not identical objects. The Hilbert space describes deterministic functions with a kernel-induced norm, whereas Gaussian sample paths commonly lie outside that space with probability one. Posterior means under Gaussian noise nevertheless coincide with solutions to certain regularized kernel-estimation problems.
The covariance matrix also defines a geometry among observations. Nearby inputs need not be strongly correlated unless the covariance function assigns them a large covariance, and distant inputs need not be weakly correlated when periodic or long-range structure is present. Consequently, the relevant notion of similarity is supplied by the kernel rather than by distance alone.
Sample-path properties
Finite-dimensional Gaussianity does not by itself guarantee that process realizations are continuous. Continuity, differentiability, and boundedness depend on the index set and on the behavior of the covariance near coincident inputs. Mean-square continuity follows when
[ \mathbb E\left[(f(x)-f(x'))^2\right]
k(x,x)+k(x',x')-2k(x,x') ]
approaches zero as (x') approaches (x). Almost-sure continuity requires additional conditions because convergence in mean square is weaker than pathwise convergence.
Differentiability can be studied through derivatives of the covariance function. When the required mixed derivatives exist and define valid covariance functions, derivatives of the process exist in the mean-square sense and are themselves Gaussian processes. Brownian motion illustrates the distinction between covariance continuity and path differentiability: its covariance
[ k(s,t)=\min(s,t) ]
is continuous, but its sample paths are almost surely nowhere differentiable.
Gaussian processes also support linear operators in a direct manner. If (L) is a suitable linear operator, then (Lf) remains Gaussian, with mean (Lm) and covariance obtained by applying the corresponding operators to the two arguments of (k). This closure property under linear transformation connects Gaussian processes with stochastic differential equations and probabilistic representations of linear physical systems.
Statistical estimation and computation
Covariance functions often contain parameters controlling marginal scale, correlation length, noise variance, or regularity. Under Gaussian observations, these parameters may be estimated through the log marginal likelihood
[ \log p(\mathbf y\mid X)
-\frac{1}{2} (\mathbf y-\mathbf m)^{\mathsf T} K_y^{-1} (\mathbf y-\mathbf m) -\frac{1}{2}\log|K_y| -\frac{n}{2}\log(2\pi), ]
where (K_y=K+\sigma_n^2I). The quadratic term measures covariance-weighted agreement with the observations, while the log-determinant term reflects the volume associated with the implied Gaussian distribution.
Exact inference for (n) observations generally relies on a factorization of the (n\times n) covariance matrix. Dense Cholesky decomposition requires cubic time complexity and quadratic storage in (n). This scaling motivates structured covariance methods, inducing-variable approximations, state-space representations, and iterative linear-algebra techniques. Such methods alter the computational representation while retaining either an exact or approximate Gaussian-process model.
When the likelihood is non-Gaussian, the posterior process is not generally Gaussian. Classification models and count-data models therefore involve posterior approximation or numerical integration. The latent prior remains a Gaussian process, but the finite-dimensional conjugacy responsible for the closed-form regression equations no longer applies.