Gaussian Processes — regression & classification

Gaussian processes (GPs) give us a nonparametric prior over functions and an analytically tractable posterior in the regression case — a small, well-conditioned linear-algebra problem at the heart of “Bayesian deep learning for the reasonable person.” This section walks through the two canonical GP inference modes — conjugate regression and non-conjugate latent-variable classification — using pyrox’s gp module. For the textbook treatment see Rasmussen & Williams (2006); for a broader probabilistic-modeling context see Murphy (2012).

Model¶

Let $x \in \mathbb{R}^{D}$ be an input and $f : \mathbb{R}^{D} \to \mathbb{R}$ the latent function of interest. A Gaussian process prior

is fully specified by a mean function $m : \mathbb{R}^D \to \mathbb{R}$ (typically zero) and a positive-definite kernel $k_\theta : \mathbb{R}^D \times \mathbb{R}^D \to \mathbb{R}$ parameterized by hyperparameters θ (length-scales, amplitude, kernel-specific extras). Evaluated at any finite collection $X = \{x_i\}_{i=1}^N$, the vector $\mathbf{f} = f(X)$ is jointly Gaussian with $\mathbf{f} \sim \mathcal{N}(\mathbf{m}_X,\, K_{XX})$ where $[K_{XX}]_{ij} = k_\theta(x_i, x_j)$.

Conjugate regression¶

For Gaussian noise $y_i = f(x_i) + \varepsilon_i$, $\varepsilon_i \overset{\text{iid}}{\sim} \mathcal{N}(0, \sigma_n^2)$, the posterior predictive at new inputs $X_\star$ is closed-form:

Hyperparameters are learned by maximizing the log marginal likelihood

Non-conjugate classification¶

For binary outputs $y_i \in \{0, 1\}$ with Bernoulli likelihood $y_i \mid f(x_i) \sim \text{Bernoulli}\bigl(\sigma(f(x_i))\bigr)$ (where σ is the logistic link), the posterior $p(\mathbf{f} \mid \mathbf{y})$ is no longer Gaussian. Classical options are a Laplace approximation MacKay (1992) — a Gaussian approximation centered at the MAP with precision equal to the Hessian of the negative log posterior — or stochastic variational inference Hensman et al. (2013). In pyrox the latent-GP model is the same equinox module as in the regression case; only the likelihood changes.

Numerical considerations¶

• Cholesky factorization is the workhorse. Factor $K_{XX} + \sigma_n^2 I = L L^\top$ once, then solve with $L$ via two triangular systems; this is $O(N^3)$ time / $O(N^2)$ memory and dominates both training and prediction.
• Jitter / nugget. $K_{XX}$ is positive definite in exact arithmetic but loses rank near the data limit. Adding $\epsilon I$ with $\epsilon \in [10^{-6}, 10^{-4}]$ is standard. Pyrox threads this through its kernel call.
• Parameterization of positives. Length-scales, amplitudes, and noise variances are constrained to $\mathbb{R}_{>0}$. Train in the unconstrained space — typically $\log$ or $\operatorname{softplus}^{-1}$ — so gradients stay finite near zero. The masterclass notebooks show three ways pyrox wires this up.
• Inducing points. For $N \gtrsim 10^4$, exact GPs are impractical. Sparse variational GPs Titsias (2009)Hensman et al. (2013) pick $M \ll N$ inducing locations and give a $O(N M^2)$ algorithm; this isn’t covered in the tutorials below but is worth knowing when they feel slow.
• Prediction variance conditioning. Computing $\Sigma_\star$ via $K_{\star X} (K + \sigma_n^2 I)^{-1} K_{X\star}$ is numerically safer when done as $K_{\star\star} - V^\top V$ with $V = L^{-1} K_{X\star}$, rather than inverting the kernel explicitly.

Notebooks¶

• exact_gp_regression — conjugate regression end-to-end: prior / marginal likelihood / posterior predictive on a 1D synthetic dataset, three patterns for wiring hyperparameters (see the masterclass sub-section).
• latent_gp_classification — latent-GP + Bernoulli likelihood via NumPyro’s NUTS sampler, with the same three patterns for parameter handling.

References¶