Gaussian Process Regression

Summary

Regression — inferring a function $f:\mathbb{X}\to \mathbb{R}$ from finitely many (noisy) evaluations — is the low-level operation inside every numerical method and the canonical PN latent object. This note develops it two ways: parametric (Gaussian prior on weights of feature functions) and nonparametric (a Gaussian process $\mathcal{G}\mathcal{P}(m,k)$, the infinite-dimensional limit). Both give closed-form Gaussian posterior mean and covariance over function values. The GP posterior mean coincides with kernel ridge regression / least-squares, and its posterior variance equals the RKHS worst-case error — the bridge PN uses to reinterpret classical methods.

Overview

In numerical applications the objects of interest are usually real-valued functions: the integrand for quadrature, the objective for optimisation, the vector field for an ODE. Applying the Gaussian inference machinery to functions is the central tool of PN. Regression is “an internal operation of all numerical methods and arguably a low-level numerical algorithm in itself.” From the ML/statistics viewpoint regression is not computation but learning — highlighting once more the fundamental similarity of computation and learning (Computation as Probabilistic Inference).

Two equivalent constructions: start finite (weights over $F$ features) then take $F\to \infty$ to obtain a Gaussian process, a Gaussian measure over an infinite-dimensional function space in which only finite marginals — governed by a mean function and a covariance kernel — are ever manipulated.

Main Content

Parametric (weight-space) Gaussian regression

Linear feature model

Assume $f$ is a weighted sum of $F$ feature functions ${\varphi }_i:\mathbb{X}\to \mathbb{R}$:

$$f(x)=\sum _{i=1}^F{\varphi }_i(x)w_i=:{\Phi }_x^{\top }w,w\in {\mathbb{R}}^F.$$

For a data set $X=[x_1,\dots ,x_N]$, write ${\Phi }_X\in {\mathbb{R}}^{F\times N}$ (columns are feature vectors), and $f_X=[f(x_1),\dots ,f(x_N)]$. The model is linear in the weights $w$ (though generally nonlinear in $x$) — a linear regression model. (Eq. 4.1.)

Place a Gaussian prior $p(w)=\mathcal{N}(w;\mu ,\Sigma )$ and observe $y=[y_1,\dots ,y_N]$ corrupted by Gaussian noise, $p(y\mid f)=\mathcal{N}(y;f_X,\Lambda )$ with $\Lambda \in {\mathbb{R}}^{N\times N}$ (the case $\Lambda \to 0$ gives exact, machine-precision observations $p(y\mid f)=\delta (y-f_X)$). By the conditioning formula (Eq. 3.13), the posterior over weights is

$$p(w\mid y)=\mathcal{N}(w;\tilde{\mu },\tilde{\Sigma }),\tilde{\Sigma }=({\Sigma }^{-1}+{\Phi }_X{\Lambda }^{-1}{\Phi }_X^{\top }{)}^{-1},\tilde{\mu }=\tilde{\Sigma }({\Sigma }^{-1}\mu +{\Phi }_X{\Lambda }^{-1}y).$$

Inference in weight space costs $\mathcal{O}(N{F}^2+F^3)$. Via the matrix inversion lemma this can be rewritten with an $N\times N$ inverse. The induced posterior over function values $f_x$ at test points is (Eq. 4.3):

$$p(f_x)=\mathcal{N}(f_x;{\Phi }_x^{\top }\tilde{\mu },{\Phi }_x^{\top }\tilde{\Sigma }{\Phi }_x),$$ $${\Phi }_x^{\top }\tilde{\mu }={\Phi }_x^{\top }\mu +{\Phi }_x^{\top }\Sigma {\Phi }_X({\Phi }_X^{\top }\Sigma {\Phi }_X+\Lambda {)}^{-1}(y-{\Phi }_X\mu ).$$

Crucially, function-space expressions contain the features only through inner products

$$k_{ab}:={\Phi }_a^{\top }\Sigma {\Phi }_b=\sum _{ij=1}^F{\varphi }_i(a){\varphi }_j(b){\Sigma }_{ij},m_a:={\Phi }_a^{\top }\mu ,$$

the covariance function $k:\mathbb{X}\times \mathbb{X}\to \mathbb{R}$ and mean function $m:\mathbb{X}\to \mathbb{R}$. The features are “encapsulated” — motivating the kernel trick and the nonparametric limit. (The maximum-likelihood/OLS estimate $w_{ML}=({\Phi }_X{\Phi }_X^{\top }{)}^{-1}{\Phi }_{X}y$ is the flat-prior special case, Exercise 4.1.)

Positive-definite kernels

Kernel (positive semi-definite)

A function $k:\mathbb{X}\times \mathbb{X}\to \mathbb{R}$ is a positive (semi-)definite kernel if for every finite set $X=[x_1,\dots ,x_N]$ the Gram matrix $[k_{XX}{]}_{ij}=k(x_i,x_j)$ is positive (semi-)definite. Such $k$ are also called Mercer kernels / covariance functions. (Def. 4.2.)

Kernels form a semi-ring: if $k,h$ are kernels and $\alpha \in {\mathbb{R}}_{+}$, then $\alpha k$, $k+h$, the Hadamard product $(k\odot h)(a,b)=k(a,b)h(a,b)$ (Schur’s theorem), and $k(\varphi (y),\varphi (y^{\prime }))$ for any map $\varphi$ are all kernels. Standard examples over $\mathbb{R}$:

$$k_{\text{SE}}(a,b)={\theta }^2\exp (-\frac{(a-b{)}^2}{2{\lambda }^2})\text{(square-exponential)},k_{\text{Wiener}}(a,b)=\min (a,b),k_{\text{OU}}(a,b)=\exp (-|a-b|).$$

The square-exponential kernel arises as the $F\to \infty$ limit of a Gaussian-feature model with feature prior variance scaled as $\propto 1/F$ (Eq. 4.4): infinitely many degrees of freedom, finite total prior mass.

Nonparametric limit: the Gaussian process

Gaussian process

Given a mean function $\mu :\mathbb{X}\to \mathbb{R}$ and a positive-definite kernel $k$, the Gaussian process $p(f)=\mathcal{G}\mathcal{P}(f;\mu ,k)$ is the measure such that, for any finite subset $X=[x_1,\dots ,x_N]$, the function values $f_X$ are jointly Gaussian, $p(f_X)=\mathcal{N}(f_X;{\mu }_X,k_{XX})$. (Def. 4.4.) It extends regression from finite weight vectors to real-valued functions.

GP regression posterior

Given a GP prior $p(f)=\mathcal{G}\mathcal{P}(f;\mu ,k)$ and Gaussian likelihood $p(y\mid f)=\mathcal{N}(y;f_X,\Lambda )$, the posterior is a GP, $p(f\mid y)=\mathcal{G}\mathcal{P}(f;m,\mathbb{V})$, with

$$m_x={\mu }_x+k_{xX}(k_{XX}+\Lambda {)}^{-1}(y-{\mu }_X),\text{\hspace{1em}(4.6)}$$ $$\mathbb{V}(x,x^{\prime })=k_{x{x}^{\prime }}-k_{xX}(k_{XX}+\Lambda {)}^{-1}{k}_{X{x}^{\prime }}.\text{\hspace{1em}(4.7)}$$

Here $k_{xX}$ is the row vector of covariances between test point $x$ and data $X$. The dominant cost is the $\mathcal{O}(N^3)$ Cholesky factorisation of the Gram matrix $G=k_{XX}+\Lambda$; the mean prediction at a new point is then $\mathcal{O}(N)$ and the marginal variance $\mathcal{O}(N^2)$.

Connection to least-squares / kernel ridge regression / RKHS

Each kernel is associated (Moore–Aronszajn) with a reproducing kernel Hilbert space (RKHS) ${\mathcal{H}}_k$: functions $f$ with $k(x,\cdot )\in {\mathcal{H}}_k$ and the reproducing property $f(x)=⟨f,k(\cdot ,x)⟩$. The GP posterior mean is exactly the kernel ridge regression estimate — the minimiser of a regularised empirical risk over ${\mathcal{H}}_k$:

$$m_x=\arg \underset{f\in {\mathcal{H}}_k}{\min }{\sigma }^2\Vert f{\Vert }_{{\mathcal{H}}_k}^2+\frac{1}{N}\sum _{i=1}^N(y_i-f(x_i){)}^2.\text{\hspace{1em}(4.8)}$$

Posterior variance = worst-case error (noise-free)

Let $\mathcal{H}$ be the RKHS of a positive-definite kernel $k$, and consider $f\in \mathcal{H}$ with $\Vert f\Vert \le 1$. For noise-free data ($\Lambda \to 0$) at $X$, the pointwise error between $f_x=f(x)$ and the posterior mean $m_x$ is tightly bounded above by the GP posterior marginal variance:

$$\underset{f\in \mathcal{H},\Vert f\Vert \le 1}{\sup }(m_x-f_x{)}^2=k_{xx}-k_{xX}{k}_{XX}^{-1}{k}_{Xx}.\text{\hspace{1em}(Thm. 4.8)}$$

Thus the Bayesian expected squared error equals the worst-case squared error over the unit RKHS ball — the linchpin identity that lets PN endow classical point estimates with calibrated error bars.

Inference on and from derivatives and integrals

Because linear maps preserve Gaussianity (Eq. 3.4), if $f\sim \mathcal{G}\mathcal{P}(m,k)$ with sufficiently differentiable/integrable $m,k$, then any partial derivatives ${\partial }^{\ell }f/\partial x^{\ell }$ and integrals ${\int }_a^{b}f\mathrm{d}\nu$ are jointly Gaussian with $f$. Their mean/covariance functions are obtained by applying the same linear operator to $m$ and to each argument of $k$, e.g.

$$\mathbb{E}\left(\frac{{\partial }^{\ell }f(x)}{\partial x^{\ell }}\right)=\frac{{\partial }^{\ell }m(x)}{\partial x^{\ell }},\mathrm{cov}(\frac{{\partial }^{\ell }f(x)}{\partial x^{\ell }},\frac{{\partial }^{m}f(x^{\prime })}{\partial x^{\prime m}})=\frac{{\partial }^{\ell }{\partial }^{m}k(x,x^{\prime })}{\partial x^{\ell }\partial x^{\prime m}}.$$

This is the engine of Bayesian Quadrature (observe $f$, infer $\int f$) and ODE filters (relate $f$ and $f^{\prime }$).

Examples

Two viewpoints on the same estimate

(1) Assign $p(f)=\mathcal{G}\mathcal{P}(f;0,k)$ with noise-free likelihood; the posterior is a GP with mean $m_x$ and variance ${\sigma }_x^2$. (2) Assume the true $f\in {\mathcal{H}}_k$ and take the regularised least-squares estimate (Eq. 4.8) — it equals $m_x$, and its distance to the truth is bounded by ${\sigma }_x$ times the (unknown) RKHS norm of $f$. The probabilistic and the approximation-theoretic stories give the same point estimate and error, but the prior is a generative object one can sample and criticise, whereas an empirical risk is not.

Integrated Wiener process posterior = cubic spline

Under the once-integrated Wiener process kernel (a Gauss–Markov prior, see Gauss-Markov Processes and SDEs) and noise-free observations, the GP posterior mean is the piecewise-cubic interpolant of the data — the classical cubic spline. This is a first concrete instance of “a classical numerical method is a Gaussian posterior mean.”

Connections

• Specialises Gaussian Distributions and Algebra to infinite dimension (marginalisation reads off finite sub-blocks).
• Univariate GPs with the Markov property are Gauss-Markov Processes and SDEs, enabling $\mathcal{O}(N)$ inference via Bayesian Filtering and Smoothing.
• Kernel/covariance choice and its hyperparameters are calibrated in Hierarchical Inference in Gaussian Models.
• The worst-case/expected-error identity underlies Convergence and Priors in Bayesian Quadrature and Uncertainty Calibration for Linear Solvers.

See Also

• Gaussian Distributions and Algebra — the finite-dimensional parent.
• Gauss-Markov Processes and SDEs — GPs with finite memory (linear-time inference).
• Bayesian Quadrature — GP regression used to infer integrals.
• Kernel Quadrature and Kernel Means — RKHS/kernel-mean view of integration.
• Classical Quadrature as Inference — trapezoid/Gauss rules as the posterior mean of this regression applied to the integrand.
• ODE Filters and Smoothers — GP regression in state-space form, conditioned on the ODE vector field.
• First- and Second-Order Optimisation Methods — BFGS as Gaussian regression on the Hessian.