#todo

Resources

• Rasmussen, Carl Edward, and Christopher K. I. Williams. 2008. Gaussian Processes for Machine Learning. 3. print. Adaptive Computation and Machine Learning. MIT Press.
• Krause, Andreas, and Jonas Hübotter. 2025. “Probabilistic Artificial Intelligence.” arXiv:2502.05244. Preprint, arXiv, February 7. https://doi.org/10.48550/arXiv.2502.05244.

Related

Main Idea

Gaussian Process

Recall that a Stochastic Process is the infinite dimension analogue of a random vector.
A Gaussian process (GP) imposes a Gaussian structure. Suppose $x$ is the domain of this function (similar to the vector index). For any points ${\mathbf{x}}_1,\dots ,{\mathbf{x}}_n$, a GP $f$ with Kernel $k(\mathbf{x},{\mathbf{x}}^{\prime })$ has function evaluations distributed according to

where $K_{ij}=k({\mathbf{x}}_i,{\mathbf{x}}_j)$ (Gram matrix), and ${\mu }_i=\mu ({\mathbf{x}}_i)$. According to prior information, this kernel can be further specified. If we assume $\mu =0$, $k$ fully defines the behavior of the process $f$. More specifically, a covariance kernel is Symmetric Positive Semidefinite.

For instance if $k$ is stationary, it can be written as $k({\mathbf{x}}_i,{\mathbf{x}}_j)=k^{\prime }({\mathbf{x}}_i-{\mathbf{x}}_j)$, and the process $f$ is stationary in both strict and wide sense because the first two moments are finite.
More restrictively, if $k({\mathbf{x}}_i,{\mathbf{x}}_j)=k^{\prime }(||{\mathbf{x}}_i-{\mathbf{x}}_j||)$, then the covariance is isotropic and a Radial Basis Function. The structure of $k$ also encodes information about the smoothness of $f$ realizations, for instance relating to Square Integrable functions ($L^2$) depending also on the spatial domain.

Weight-space view

Let us take a step back to the discrete case before continuing with the functional version afterwards.

Let the independent training data be $\mathcal{D}=\{({\mathbf{x}}_i,y_i)|i=1,\dots ,n\}\subset {\mathbb{R}}^D\times \mathbb{R}$, or the matrix version with design matrix $\mathbf{X}\in {\mathbb{R}}^{D\times n}$ and target vector $y\in {\mathbb{R}}^n$. We are interested in modeling the distribution of $y$, given some $\mathbf{x}$ (e.g. the conditional distribution $p(y|\mathbf{x})$, not a generative distribution $p(\mathbf{x})$ or $p(\mathbf{x}|y)$).

The standard Bayesian Linear Regression paradigm poses model $f$ with parameters $\mathbf{w}\in {\mathbb{R}}^d$ and a Gaussian noise model on the observations:

The model and noise assumption define the likelihood, a probability distribution over the data observations, given the model. This is

In summary, this gives the likelihood as a Multivariate Gaussian:

The Bayesian framework also requires a prior over the parameter values, which is assumed to also be Gaussian, centered at $0$ without loss of generality (we could normalize the data instead), which for covariance ${\mathrm{\Sigma }}_p$ is

Deriving parameter posterior

The posterior (the distribution over the parameters, given the data and observations) is computed according to Bayes Theorem as

The marginal likelihood in the denominator is independent of the parameters $\mathbf{w}$, and is computed according to the (usually intractable) integral

While in this case, one could conceivably compute this, there is no need, as we determine the shape of the posterior as Gaussian and can normalize afterwards because of this nice distribution. Thus, we continue by noting

We have the argument of $\exp$ as

which we want to write in the Gaussian form by finding some unknown mean $\overline{\mathbf{w}}$ such that

Setting these two equal, we can drop the ${\overline{\mathbf{w}}}^T\mathbf{A}\overline{\mathbf{w}}$ term, which is constant, and find

Thus, the posterior is

Expanding the parameter mean (and mode / Maximum a Posteriori) gives

Recall, the normal equations from (linear) Least Squares give the parameters as

or the version with Tikhonov Regularization or Ridge Regression magnitude $\lambda$:

We can relate this Bayesian version by noting that $\lambda \mathbf{I}={\mathrm{\Sigma }}_p^{-1}$ gives the same coefficient prediction (ignoring $\sigma$); in other words, the Gaussian prior regularizes the coefficients in exactly the same way as Tikhonov Regularization (we can generalize to non-diagonal regularization if needed). Alternatively, if we have ${\mathrm{\Sigma }}_p^{-1}=\mathbf{0}$, e.g. by infinite variances for the prior, then we recover the least squares case.

Predictive posterior

Future predictions for new $x_{\ast }$ include pushing the posterior distribution through the forward model with that $x_{\ast }$ (i.e. classic Forward UQ). In other words, denote the predictive distribution as $f_{\ast }$, which has true value $f(x_{\ast })$. The distribution is computed as

That is, sample $\mathbf{w}$ from the posterior. Then with that value, compute $f_{\ast }$ by also using $x_{\ast }$, which is a deterministic map in this case. Because $f$ is linear, this output is also Gaussian, with mean and covariance adjusted according to the respective transform:

To predict the distribution of the observations, rather than the underlying function, we would also add in our Gaussian noise model (e.g. add ${\sigma }^2$ in the variance). This would give

The variance over these predicted observations can be broken down into epistemic and aleatoric components. The epistemic uncertainty corresponds to ${\mathbf{x}}_{\ast }^T{\mathbf{A}}^{-1}{\mathbf{x}}_{\ast }$, which can be reduced with more data. It is uncertainty about the model itself. On the other hand, aleatoric uncertainty encapsulates the ${\sigma }^2\mathbf{I}$ term, which is noise inherent in the observations themselves, and is irreducible.

Feature space

We can generalize to a different class of models by instead "engineering" features depending on $\mathbf{x}$, via the vector-valued $\varphi (\mathbf{x})$, i.e. $\varphi :{\mathbb{R}}^D\to {\mathbb{R}}^N$. The matrix $\mathrm{\Phi }(\mathbf{X})\in {\mathbb{R}}^{N\times n}$ comes from applying $\varphi$ to each column of $\mathbf{X}$. We now seek $\mathbf{w}\in {\mathbb{R}}^N$, with one parameter for each feature. The new model is nonlinear in $\mathbf{x}$:

This is fundamentally just changing what inputs we build a linear model for, substituting $\mathrm{\Phi }$ for $\mathbf{X}$ in the analysis. The predictive posterior distribution changes the most (and is the most relevant, as it is more connected with ${\mathbf{x}}_{\ast }$), giving

with $\mathbf{A}={\sigma }^{-2}\mathrm{\Phi }{\mathrm{\Phi }}^T+{\mathrm{\Sigma }}_p^{-1}$. In this case, predictions, even deterministic ones, hinge on inverting the $N\times N$ matrix $\mathbf{A}$, which may be difficult for large $N$.

Kernel Trick

Let $\mathbf{K}={\mathrm{\Phi }}^T{\mathrm{\Sigma }}_p\mathrm{\Phi }$, noting $\mathbf{K}\in {\mathbb{R}}^{n\times n}$. We will derive a new form of the predictive posterior:

Now, left multiply by ${\mathbf{A}}^{-1}$ and right multiply by $(\mathbf{K}+{\sigma }^2\mathbf{I}{)}^{-1}$ to give

Thus, we can rewrite the mean as

The covariance matrix can be rewritten using the Sherman-Morrison inversion formula for ${\mathbf{A}}^{-1}=({\mathrm{\Sigma }}_p^{-1}+{\sigma }^{-2}\mathrm{\Phi }{\mathrm{\Phi }}^T{)}^{-1}$, giving

Then the entire covariance term just plugs this in to give

The main point of this tedious derivation is to have the inverse on $(\mathbf{K}+{\sigma }^2\mathbf{I}{)}^{-1}$, which is $n\times n$, rather than the $N\times N$ previous version ${\mathbf{A}}^{-1}$. Thus, the cost of computing the posterior predictive (propagating parameter uncertainties) scales independently of the feature / parameter dimension $N$, instead scaling based on the data size $n$. This is uncommon for uncertainty propagation, where we expect cost to scale based on the parameter dimension. Also with many $\varphi$, we expect $N$ to scale quite poorly. For instance if $\varphi$ is a polynomial expansion of degree $p$, then $N=\mathcal{O}(d^p)$.

In the final forms for the predictive mean and covariance, the feature space always comes in as ${\mathrm{\Phi }}^T{\mathrm{\Sigma }}_p\mathrm{\Phi }$ or $\varphi ({\mathbf{x}}_{\ast }{)}^T{\mathrm{\Sigma }}_p\varphi ({\mathbf{x}}_{\ast })$ (or a mix between the two). This motivates the definition of a covariance function or kernel $k(\mathbf{x},{\mathbf{x}}^{\prime })=\varphi (\mathbf{x}{)}^T{\mathrm{\Sigma }}_p\varphi ({\mathbf{x}}^{\prime })$. Recalling that the prior covariance matrix on the parameters, ${\mathrm{\Sigma }}_p$, is positive definite, we can also define $\psi (\mathbf{x})={\mathrm{\Sigma }}_p^{1/2}\varphi (\mathbf{x})$ and write the kernel as an Inner Product (or dot product): $k(\mathbf{x},{\mathbf{x}}^{\prime })=\psi (\mathbf{x}{)}^T\psi (\mathbf{x})=\psi (\mathbf{x})\cdot \psi (\mathbf{x})$. Fundamentally, the kernel trick will rely on this kernel to replace inner products. As we have derived in this section, this substitution avoids explicitly calculating the feature vectors themselves $\varphi (\mathbf{x})$, even allowing for $\varphi$ to conceptually have an infinite-dimensional output.

Function-space view

Previously we had the model

where we treated $\mathbf{w}$ as a random vector (in a Bayesian sense). This then means that the function $f$ is a random function. This is like a standard basis function expansion from coefficients ($\mathbf{w}$) to an actual function.

Through the kernel trick, we already eliminated the direct influence of $\varphi$ and $\mathrm{\Sigma }$, replacing it with the kernel. If we just model $f(\mathbf{x})$, we can still presumably do the same things as before (generate a predictive posterior). Because $\mathbf{w}$ is Gaussian (in both prior and posterior form) and $f$ is a linear combination of $\mathbf{w}$, $f$ is also Gaussian. Hence, we define it directly as a Gaussian process (in a function space).

To match the above case, we can define $f(\mathbf{x})=\varphi (\mathbf{x}{)}^T\mathbf{w}$ for prior $\mathbf{w}\sim \mathcal{N}(0,{\mathrm{\Sigma }}_p)$, giving the Gaussian process mean and covariance

This converts from the weight-space view to the function-space view.

We can alternatively start directly in the function-space view, defining $f$ as a Gaussian process, which can induce a certain meaning to $\varphi (\mathbf{x})$ and $\mathbf{w}$. For instance, the squared exponential kernel

corresponds implicitly to infinitely many basis functions $(N=\mathrm{\infty })$. Via Mercer's Theorem, any positive definite covariance kernel induces a (possibly infinite) dimensional basis function expansion $\varphi$. The choice of the kernel (including the length scale $\ell$) specifies the distribution over functions.