Gaussian Processes & RKHS

1. The Infinite-Dimensional Gaussian

A Gaussian Process is a big pile of random variables indexed by a set $X$ and every finite linear combination of them drops out normal. The definition is short and clean but making it rigorous leans on the Kolmogorov Extension Theorem to stitch together a consistent measure on $\mathbb{R}^X$. The whole process sits on two pieces and that is all you need.

1. Mean function: $m(x) = \mathbb{E}[f(x)]$.
2. Covariance kernel: $k(x, x') = \mathbb{E}[(f(x)-m(x))(f(x')-m(x'))]$.

For the finite-dimensional distributions to hang together the kernel $k$ has to be positive semi-definite and for any set of points $\{x_i\}$ the Gram matrix $K_{ij} = k(x_i, x_j)$ has to satisfy $v^T K v \ge 0$ for every $v \in \mathbb{R}^n$. And that is the whole constraint.

2. The Reproducing Kernel Hilbert Space (RKHS)

Every PSD kernel $k$ carries with it a unique Hilbert space $\mathcal{H}_k$ of functions $f: X \to \mathbb{R}$ and this space is tiny. It holds the deterministic skeleton of the random process and wipes out the sample paths themselves with probability one.

You build it in a few steps and each one is short.

1. Define the pre-Hilbert space $\mathcal{H}_0 = \text{span} \{ k(\cdot, x) \mid x \in X \}$.
2. Define the inner product on basis elements: $$\langle k(\cdot, x), k(\cdot, y) \rangle_{\mathcal{H}_k} = k(x, y)$$
3. Extend by linearity and complete the space to get $\mathcal{H}_k$.

The Reproducing Property. For any $f \in \mathcal{H}_k$ you get

The evaluation functional $\delta_x: f \mapsto f(x)$ is a continuous linear operator and this is what makes point evaluation behave.

Mercer’s Theorem. Take $X$ compact and $k$ continuous and look at the integral operator $T_k: L^2(X) \to L^2(X)$.

This operator is compact and positive and self-adjoint and it comes with an orthonormal basis of eigenfunctions $\{\phi_i\}$ and eigenvalues $\lambda_i \ge 0$ and these let you write the kernel as a sum.

The RKHS norm then drops into a clean spectral form.

The RKHS norm punishes high-frequency pieces where $\lambda_i \to 0$ and the smaller the eigenvalue the more painful it is for a function to carry energy in that mode. Sample paths of the GP usually blow up the RKHS norm because $\sum \frac{|\xi_i|^2 \lambda_i}{\lambda_i} = \sum \xi_i^2$ diverges with $\xi_i \sim \mathcal{N}(0, 1)$ and this is the Driscoll Zero-One Law which says $P(f \in \mathcal{H}_k) = 0$. The sample paths actually sit in a slightly bigger space and that space is usually a fractional Sobolev space whose smoothness is set by how fast the spectrum of the kernel drops off.

3. The Geometry of the Posterior

You see noisy data $y_i = f(x_i) + \epsilon_i$ and you want the posterior $f_* | y$. On paper this boils down to Gaussian conditioning and geometrically it is an orthogonal projection sitting inside the $L^2(\Omega)$ space of random variables.

Look at the subspace spanned by the observations $V_n = \text{span}\{f(x_1), \dots, f(x_n)\}$. The variables are jointly Gaussian and so the conditional expectation $\mathbb{E}[f(x_*) | y]$ drops out as the orthogonal projection of $f(x_*)$ onto the subspace spanned by $y$.

This lines up with the Representer Theorem and the posterior mean sits inside the subspace spanned by the kernel sections $\{k(\cdot, x_i)\}$ and nowhere else.

The posterior variance is the squared norm of the residual and it picks out the piece of $f(x_*)$ that sits orthogonal to the observed subspace.

The geometry here is easy to see. Variance shrinks near data points because the angle between $f(x_*)$ and the observed variables $f(x_i)$ gets small and the projection grabs almost all of the norm and leaves almost nothing behind.

4. Spectral Kernels and Bochner’s Theorem

Theorem (Bochner). A continuous and shift-invariant function $k(\tau)$ on $\mathbb{R}^d$ is positive definite exactly when it is the Fourier transform of a finite non-negative measure $S(\omega)$.

Since $k$ is real the density $S(\omega)$ has to be symmetric and $S(\omega) = S(-\omega)$ and this $S(\omega)$ is what people call the spectral density.

Spectral Mixture Kernels (Wilson & Adams, 2013) come from a simple observation. Standard kernels match up with single-bump spectral densities centered at the origin and that is all they can do.

• RBF: $S(\omega) \propto e^{-\|\omega\|^2}$ (Gaussian).
• Matern: $S(\omega) \propto (1+\|\omega\|^2)^{-\nu-d/2}$ (Student-t / Cauchy-like). These densities fall off monotonically and they cannot hold negative covariances and they cannot hold periodic structure.

The set of Gaussian mixtures is dense in the space of distribution functions under weak convergence and so you can approximate any spectral density this way.

Taking the inverse Fourier transform hands you the Spectral Mixture kernel.

The parameters are easy to read off the expression.

• $w_q$: Contribution of the $q$-th mode.
• $\Sigma_q$: Inverse lengthscale (frequency variance).
• $\mu_q$: Frequency of the periodicity.

By learning the weights and means and variances of the spectral density a GP ends up doing non-parametric Fourier analysis and it handles noise and missing data at the same time. Learning the spectral density swaps out manual kernel selection for kernel learning straight from the data and this is a big deal.

5. Scaling: Inducing Points (Sparse GPs)

The $O(N^3)$ cost of exact inference pushes you toward low-rank approximations of the kernel matrix. Nystrom approximations that just subsample the matrix are numerically flaky and the modern fix is the Variational Sparse GP (Titsias, 2009) which treats the inducing points as variational parameters and tunes them.

The augmented model starts by dropping in $M \ll N$ inducing inputs $Z = \{z_1, \dots, z_M\}$ and inducing variables $u = f(Z)$ and the joint prior factors cleanly.

with $p(u) = \mathcal{N}(0, K_{uu})$ and $p(f|u) = \mathcal{N}(K_{fu} K_{uu}^{-1} u, K_{ff} - Q_{ff})$ where $Q_{ff} = K_{fu} K_{uu}^{-1} K_{uf}$.

For variational inference you approximate the true posterior $p(f, u | y)$ with a variational distribution $q(f, u)$ and the key move is to set $q(f, u) = p(f|u) q(u)$. You hold the exact conditional $p(f|u)$ fixed and only approximate the inducing distribution $q(u) = \mathcal{N}(m, S)$ and then you minimize the KL divergence $\text{KL}(q \| p)$ which is the same thing as maximizing the Evidence Lower Bound.

The $p(f|u)$ terms cancel top and bottom and you are left with this.

The bound splits into two pieces and both are clean.

1. Expected log likelihood. $$\sum_{i=1}^N \mathbb{E}_{q(f_i)} [\log \mathcal{N}(y_i | f_i, \sigma^2)]$$ Since $q(f_i)$ is Gaussian as the marginal of $q(f)$ this one is analytically tractable and it only touches the diagonal variances $k(x_i, x_i)$ and nothing else.
2. KL regularizer. $$\text{KL}(\mathcal{N}(m, S) \| \mathcal{N}(0, K_{uu}))$$ This pulls the inducing distribution back toward the prior whenever it drifts away.

Computing the ELBO forces you to invert $K_{uu}$ which costs $O(M^3)$ and multiply by $K_{uf}$ which costs $O(NM^2)$ and the total comes out to $O(NM^2)$. This stretches GPs out to millions of points once you mini-batch the expected log-likelihood sum.

6. Deep Kernel Learning and Implementation

Standard GPs run with fixed kernel families. Deep Kernel Learning glues a base kernel onto a learned feature map so $k(x, y) = k_{base}(\phi_\theta(x), \phi_\theta(y))$ where $\phi_\theta$ is a neural network and the parameters $\theta$ get tuned by pushing up the marginal likelihood.

Computing the gradient forces you to differentiate through the Cholesky decomposition and use $\partial K^{-1} = -K^{-1} (\partial K) K^{-1}$ and modern frameworks like JAX and GPyTorch handle this for you.

import jax
import jax.numpy as jnp
from jax import random, jit, grad, value_and_grad
import matplotlib.pyplot as plt
 
# -------------------------------------------------------------
# 1. Spectral Mixture Kernel Implementation
# -------------------------------------------------------------
def spectral_mixture_kernel(x1, x2, weights, means, variances):
    """
    Computes SM Kernel: k(tau) = sum w_q * exp(-2pi^2 tau^T Sigma tau) * cos(2pi mu^T tau)
    x1, x2: (N, D)
    weights: (Q,)
    means: (Q, D)
    variances: (Q, D)
    """
    # Tau: distance vectors
    tau = x1[:, None, :] - x2[None, :, :] # (N, N, D)
    tau_sq = tau ** 2
    
    # We sum over Q components
    # Vectorized computation for speed
    def component_func(w, m, v):
        # exp term: exp(-2 pi^2 * tau^2 * v)
        exp_term = jnp.exp(-2.0 * jnp.pi**2 * jnp.sum(tau_sq * v, axis=-1))
        # cos term: cos(2 pi * mu * tau)
        cos_term = jnp.cos(2.0 * jnp.pi * jnp.sum(tau * m, axis=-1))
        return w * exp_term * cos_term
        
    # vmap over the Q mixture components
    components = jax.vmap(component_func)(weights, means, variances)
    return jnp.sum(components, axis=0)
 
# -------------------------------------------------------------
# 2. Exact GP Implementation
# -------------------------------------------------------------
class ExactGP:
    def __init__(self, kernel_fn):
        self.kernel_fn = kernel_fn
        
    def loss_fn(self, params, X, y):
        # Unpack params
        noise = jnp.exp(params['log_noise'])
        # Kernel
        K = self.kernel_fn(X, X, params) + noise * jnp.eye(X.shape[0])
        # Cholesky
        L = jnp.linalg.cholesky(K)
        # Solve
        alpha = jax.scipy.linalg.cho_solve((L, True), y)
        
        # MLL
        term1 = -0.5 * jnp.dot(y.flatten(), alpha.flatten())
        term2 = -jnp.sum(jnp.log(jnp.diag(L)))
        term3 = -0.5 * len(y) * jnp.log(2 * jnp.pi)
        return -(term1 + term2 + term3) # Return Negative Log Likelihood
 
    def predict(self, params, X_train, y_train, X_test):
        noise = jnp.exp(params['log_noise'])
        K_xx = self.kernel_fn(X_train, X_train, params) + noise * jnp.eye(X_train.shape[0])
        K_xs = self.kernel_fn(X_train, X_test, params)
        K_ss = self.kernel_fn(X_test, X_test, params)
        
        L = jnp.linalg.cholesky(K_xx)
        alpha = jax.scipy.linalg.cho_solve((L, True), y_train)
        
        mu = K_xs.T @ alpha
        v = jax.scipy.linalg.solve_triangular(L, K_xs, lower=True)
        var = K_ss - v.T @ v
        return mu, jnp.diag(var)
 
# -------------------------------------------------------------
# 3. Training Loop
# -------------------------------------------------------------
def train_sm_gp(X, y, num_components=4, steps=1000):
    key = random.PRNGKey(42)
    D = X.shape[1]
    
    # Initialize parameters randomly (crucial for SM kernel)
    # Mixture weights, frequencies, lengthscales
    params = {
        'log_weights': jnp.log(jnp.ones(num_components) / num_components),
        'log_means': random.normal(key, (num_components, D)),
        'log_vars': random.normal(key, (num_components, D)) - 2.0, # Start with small bandwidth
        'log_noise': jnp.log(0.1)
    }
    
    # Wrapper to unpack dict to args for kernel
    def wrapped_kernel(x1, x2, p):
        return spectral_mixture_kernel(
            x1, x2, 
            jnp.exp(p['log_weights']), 
            jnp.exp(p['log_means']), 
            jnp.exp(p['log_vars'])
        )
        
    gp = ExactGP(wrapped_kernel)
    optimizer = jax.experimental.optimizers.adam(1e-2)
    
    @jit
    def update(i, opt_state):
        p = optimizer.params_fn(opt_state)
        grads = grad(gp.loss_fn)(p, X, y)
        return optimizer.update_fn(i, grads, opt_state)
        
    opt_state = optimizer.init_fn(params)
    for i in range(steps):
        opt_state = update(i, opt_state)
        
    final_params = optimizer.params_fn(opt_state)
    return gp, final_params
 
# The SM Kernel allows the GP to extrapolate patterns (like cardiac rhythm or airline passengers)
# far beyond the training window, unlike the RBF kernel which reverts to the mean.

7. Approximate Inference for Classification

Everything above assumed Gaussian observation noise $y = f(x) + \epsilon$. For binary classification with $y \in \{-1, 1\}$ the likelihood stops being Gaussian.

The posterior $p(f|y) \propto p(y|f) p(f)$ stops being Gaussian and the marginal likelihood $\int p(y|f) p(f) df$ turns intractable and you have to fall back on approximate inference.

7.1 The Laplace Approximation. You approximate the posterior with a Gaussian $q(f) = \mathcal{N}(\hat{f}, A^{-1})$ centered at the MAP estimate and that is the whole idea.

1. Find the mode $\hat{f} = \text{argmax} (\log p(y|f) + \log p(f))$.
2. Compute the Hessian of the negative log posterior at the mode. $$A = - \nabla \nabla \log p(f|y) |_{\hat{f}} = K^{-1} + W$$ where $W$ is the diagonal matrix of second derivatives of the log-likelihood and for the logistic likelihood you get $W_{ii} = \pi_i(1-\pi_i)$ and nothing messier than that.
3. Use $q(f)$ for predictions. $$\mathbb{E}[f_*] = k_*^T K^{-1} \hat{f}$$ The predictive variance soaks up the curvature through $W$ and that is where the uncertainty lives.

7.2 Expectation Propagation. The Laplace approximation is local and only grabs the curvature at the mode while EP gives you a global approximation by matching moments. Each likelihood factor $p(y_i|f_i)$ gets replaced by an unnormalized Gaussian site function $t_i(f_i)$.

You iterate on the site parameters by matching the first and second moments of the cavity distributions. EP usually gives better calibrated predictive probabilities than the Laplace approximation for GP classification and you pay for it with extra iterations and more bookkeeping.

8. The Bayesian Ideal for Uncertainty

Gaussian Processes give you a principled Bayesian framework for regression and function estimation and they handle uncertainty and smoothness and finite-sample corrections in a mathematically exact way. The $O(N^3)$ cost used to be a blocker and sparse variational GPs and spectral kernels and deep kernel learning have stretched them out to datasets with millions of observations.

The deeper point is the link to neural networks. Training a wide enough neural network in the infinite-width limit is the same thing as kernel regression in an architecture-dependent RKHS and so understanding GPs is a prerequisite for understanding the inductive biases and the generalization behavior of modern deep learning.

Timeline

Deriving the Matern Spectral Density

The Matern kernel pops out as the Green’s function of a stochastic partial differential equation. Take the SPDE driven by white noise $W(x)$ and watch what happens.

This is a fractional Helmholtz equation and taking the Fourier transform of both sides untangles it right away.

Since $W$ is white noise its spectral density is constant and $S_W(\omega) = \text{const}$ and the power spectral density $S_f(\omega)$ of the solution drops out as the squared modulus of the transfer function.

This is exactly the spectral density of the Matern class and the spatial kernel $k(r)$ comes back by running the inverse Fourier transform.

Using Schwinger parametrization and Bessel function identities this integral collapses down to a clean form.

The physics reads simply. At $\nu=1/2$ the SPDE is first-order and you get the Ornstein-Uhlenbeck process which is Brownian motion with friction. As $\nu \to \infty$ the operator slides toward the heat operator and the kernel drifts toward the Gaussian RBF which is the diffusion kernel.

The Neural Tangent Kernel (NTK)

Neal’s Theorem (1994). A single-layer Bayesian neural network at infinite width turns into a Gaussian Process. Let $f(x) = \frac{1}{\sqrt{H}} \sum_{j=1}^H v_j \sigma(w_j^T x + b_j)$ and by the Central Limit Theorem at fixed $x$ the output $f(x)$ is a sum of i.i.d. terms and drops to Gaussian in the limit.

The Modern Extension (Jacot et al., 2018). Look at the training dynamics of an infinite-width network under gradient descent. The network moves as a linear model against a fixed kernel $\Theta(x, x')$ which is the Neural Tangent Kernel.

In the infinite-width limit $\Theta$ holds constant all the way through training and the training dynamics drop out exactly as the kernel regression solution.

So training an infinite-width network is really kernel regression with an architecture-dependent kernel and you can work these kernels out analytically for specific activations and architectures like ReLU networks and CNNs. For a ReLU network the arc-cosine kernel falls out.

with $\phi$ the angle between inputs. The kernel nails down the inductive bias of the architecture by picking out the RKHS norm that the network implicitly regularizes.

To derive the arc-cosine kernel look at a single neuron $h(x) = \sigma(w^T x)$ with $w \sim \mathcal{N}(0, I)$. The kernel is $k(x, y) = \mathbb{E}_w [\sigma(w^T x) \sigma(w^T y)]$ and you take $\sigma(z) = \max(0, z)$ for ReLU. By rotational invariance you can spin $x$ onto the first coordinate axis and drop $y$ into the plane spanned by the first two axes so $x = \|x\| e_1$ and $y = \|y\| (\cos \phi \, e_1 + \sin \phi \, e_2)$. Writing $w = (w_1, w_2, \dots)$ the product of activations turns into a concrete expression. $\max(0, \|x\| w_1) \max(0, \|y\| (w_1 \cos \phi + w_2 \sin \phi))$ You then have to evaluate the integral that comes next.

The integration domain is the sector where both arguments come out positive and passing to polar coordinates $(u, v) = r(\cos \theta, \sin \theta)$ splits it wide open.

The radial integral gives $\int_0^\infty r^3 e^{-r^2/2} dr = 2$ and the angular integral hands back the factor $\frac{1}{2\pi} (\sin \phi + (\pi - \phi) \cos \phi)$. This result (Cho & Saul, 2009) says ReLU networks measure similarity by the angle between data points and not by Euclidean distance the way the RBF kernel does. This angular dependence is part of why they do so well on high-dimensional data like images where Euclidean distances lose information because of concentration of measure.

Gaussian Process Latent Variable Models (GPLVM)

Standard PCA assumes data $Y$ comes from a linear mapping of latent variables $X$ and $Y = WX + \epsilon$. Dual Probabilistic PCA (Lawrence, 2004) integrates $W$ out and optimizes $X$ instead and this gives you a Gaussian Process mapping from latent space to data space and $Y \sim \mathcal{N}(0, K(X, X) + \sigma^2 I)$. The GPLVM tunes the latent locations $X$ to push up the likelihood of the data $Y$.

The objective looks like this and it is short.

Since $X$ enters through the non-linear kernel $K$ this optimization problem is not convex. The kernel you pick sets the class of manifold structures you can recover and a linear kernel recovers PCA exactly and an RBF kernel gives you a non-linear embedding that holds onto local distances and it behaves a lot like t-SNE but sitting inside a probabilistic framework.

Bayesian GPLVM (Titsias & Lawrence, 2010) goes further. The standard GPLVM treats $X$ as point estimates at the MAP and the Bayesian version drops a prior $X \sim \mathcal{N}(0, I)$ on it and approximates the posterior $q(X)$ variationally. This forces you to build a variational lower bound that accounts for the kernel matrix $K_{ff}$ being a random variable because it depends on the stochastic inputs $X$. Titsias worked out closed-form expressions for the required $\Psi$-statistics.

For the RBF kernel these expectations are analytically tractable and they collapse to Gaussian convolutions. This framework opens the door to Automatic Relevance Determination on the latent dimensions and it picks out the intrinsic dimensionality of the data on its own.

Terminology

References

1. Rasmussen, C. E., & Williams, C. K. I. (2006). “Gaussian Processes for Machine Learning”. The standard text and it covers regression and classification and covariance functions and theoretical properties.

2. Wilson, A. G., & Adams, R. P. (2013). “Gaussian Process Kernels for Pattern Discovery and Extrapolation”. Introduced Spectral Mixture Kernels and showed that standard kernels are just local averages while spectral kernels can hold long-range structure and negative covariance.

3. Titsias, M. (2009). “Variational Learning of Inducing Variables in Sparse Gaussian Processes”. The breakthrough that made GPs scalable. Instead of Nystrom which ignores data Titsias treated inducing points as variational parameters to be tuned.

4. Jacot, A. et al. (2018). “Neural Tangent Kernel: Convergence and Generalization in Neural Networks”. The landmark paper linking infinite width Neural Networks to Kernel methods and giving exact solvability for Deep Learning dynamics.

5. Lawrence, N. (2004). “Gaussian Process Latent Variable Models for Visualisation of High Dimensional Data”. Introduced the GPLVM framework and flipped the GP regression problem to tune inputs rather than weights.

6. MacKay, D. J. C. (1998). “Introduction to Gaussian Processes”. A classic tutorial by the late Sir David MacKay and it connects GPs to Neural Networks early on and gives a Bayesian perspective on regularization.

7. Quinonero-Candela, J., & Rasmussen, C. E. (2005). “A Unifying View of Sparse Approximate Gaussian Process Regression”. The paper that organized the confusing zoo of sparse approximations like FITC and DTC and PITC into a single effective theory.