GP Tutorial Master List

A reconciled, exhaustive curriculum spanning what currently exists in gaussx, pyrox, and research_notebook, plus gaps surfaced from the gaussx + pyrox public APIs, open GitHub issues, and pyrox design_docs/. Goal: the most complete GP tutorial sequence we could ship.

Bayesian NN / NeRF / basis-function-regression tutorials live in ../bayesian_nns/TUTORIAL_MASTER_LIST.md. Cross-listed items (RFF, deep kernels, BLR, last-layer-Bayes) are flagged 🔁.

Legend — Source columns:

• G = exists in gaussx (docs/notebooks/<name>)
• P = exists in pyrox (docs/notebooks/<name>)
• R = exists in research_notebook (projects/gaussian_processes/notebooks/<path>)
• — = does not exist yet (gap)

Scope tag: 🧱 fundamental · 🔬 research · 🌉 bridge · 🔁 cross-listed

Refs column: gh#N = open GitHub issue · dd:path = pyrox design_docs/pyrox/<path> · api:foo = gaussx exported symbol.

Curriculum at a glance¶

A bird’s-eye view of the parts and their subparts. Skim this first to orient; the detailed per-tutorial tables live below.

• Part 0 — Linear Algebra & Gaussian Foundations
  • 0.A — The Multivariate Gaussian
  • 0.B — Parameterizations
  • 0.C — Bayesian Updates & Conditioning
  • 0.D — Numerical Mechanics
• Part 1 — Structured Linear Operators
  • 1.A — Operator Zoo (catalog)
  • 1.B — Matrix Identities & Decompositions
  • 1.C — Matrix-Free / Implicit
  • 1.D — Solvers
  • 1.E — Trace, Log-Det, Roots
• Part 2 — Kernels
  • 2.A — Standard kernels
  • 2.B — Spectral & deep kernels
  • 2.C — Multi-output kernels
  • 2.D — Spherical / localized kernels
  • 2.E — Kernel-based statistics & utilities
  • 2.F — Non-Euclidean & operator-valued kernels
• Part 3 — Exact GP Regression
  • 3.A — Foundations
  • 3.B — Diagnostics
  • 3.C — Heteroscedastic noise
  • 3.D — High-level API
  • 3.E — Constrained & Physics-informed GPs
• Part 4 — Structured GPs
  • 4.A — Kronecker GPs
  • 4.B — Grid / Toeplitz GPs
  • 4.C — Sparse-precision (mesh / GMRF)
• Part 5 — Approximations & Scalability
  • 5.A — Random features
  • 5.B — Inducing-point fundamentals
  • 5.C — Inter-domain features
  • 5.D — Iterative-solver scaling
  • 5.E — Deep GPs
• Part 6 — Non-Conjugate Likelihoods & Inference
  • 6.A — Likelihood & integrator zoos
  • 6.B — Classification
  • 6.C — Newton / Gauss-Newton family
  • 6.D — Variational inference
  • 6.E — Expectation Propagation
  • 6.F — Bayesian linear regression & non-standard outputs
  • 6.G — Aggregate Bayesian methods
• Part 7 — Spectral GPs
  • 7.A — Spectral foundations
  • 7.B — Spectral kernel models
  • 7.C — Random Fourier features
  • 7.D — Hilbert-space methods
  • 7.E — Variational spectral methods
  • 7.F — Spherical / periodic spectral
• Part 8 — Markov / State-Space GPs
  • 8.A — Foundations
  • 8.B — SDE kernel zoo
  • 8.C — Markov GP workflows
  • 8.D — Parallel & scalable filtering
  • 8.E — Nonlinear filtering
  • 8.F — Ensemble methods
  • 8.G — Steady-state & structured-Gaussian surfaces
  • 8.H — Non-conjugate temporal case studies
• Part 9 — Sampling, Pathwise, Conditioning
  • 9.A — Pathwise sampling
  • 9.B — Matheron’s-rule conditioning
• Part 10 — Uncertainty Propagation & UQ
  • 10.A — Foundations
  • 10.B — Uncertain inputs
  • 10.C — Analytic moments
  • 10.D — BGPLVM
  • 10.E — Special integrators & quantiles
• Part 11 — Probabilistic Programming Integration
  • 11.A — gaussx + NumPyro
  • 11.B — pyrox patterns
  • 11.C — Hierarchical & sampling
• Part 12 — Ensembles
• Part 13 — Data Pipelines
• Part 14 — Applied Case Studies (research_notebook)
  • 14.A — Spatial extremes
  • 14.B — SVGP applied
  • 14.C — Geophysics & emulation
  • 14.D — Optimization & decision
  • 14.E — Causal & event data
  • 14.F — Practical
• Part 15 — Metrics & Calibration

Part 0 — Linear Algebra & Gaussian Foundations¶

0.A — The Multivariate Gaussian¶

Key equations / models:

• Density: $\mathcal{N}(x;\mu,\Sigma) = (2\pi)^{-d/2}|\Sigma|^{-1/2}\exp\!\big({-}\tfrac{1}{2}(x-\mu)^\top\Sigma^{-1}(x-\mu)\big)$
• Reparameterized sample: $x = \mu + L\epsilon$, $LL^\top = \Sigma$, $\epsilon\sim\mathcal{N}(0,I)$
• Entropy: $H = \tfrac{1}{2}\log|2\pi e\,\Sigma|$
• KL: $\mathrm{KL}(p\Vert q) = \tfrac{1}{2}\big[\mathrm{tr}(\Sigma_q^{-1}\Sigma_p) + \Delta\mu^\top\Sigma_q^{-1}\Delta\mu - d + \log\tfrac{|\Sigma_q|}{|\Sigma_p|}\big]$

0.B — Parameterizations¶

Key equations / models:

• Natural parameters: $\eta_1 = \Sigma^{-1}\mu$, $\eta_2 = -\tfrac{1}{2}\Sigma^{-1}$
• Expectation parameters: $m_1 = \mu$, $m_2 = \Sigma + \mu\mu^\top$

0.C — Bayesian Updates & Conditioning¶

Key equations / models:

• Sequential conjugate update: $p(\theta\mid y_{1:n}) \propto p(y_n\mid\theta)\,p(\theta\mid y_{1:n-1})$
• Schur conditional: $\mu_{a\mid b} = \mu_a + \Sigma_{ab}\Sigma_{bb}^{-1}(x_b-\mu_b)$, $\Sigma_{a\mid b} = \Sigma_{aa} - \Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba}$
• Structured-MVN sample: $x = \mu + \mathrm{root}(\Sigma)\,\epsilon$ with dispatched root per operator type

0.D — Numerical Mechanics¶

Key equations / models:

• Joseph-form covariance update: $P^+ = (I-KH)\,P\,(I-KH)^\top + KRK^\top$ (PSD-preserving)
• Cholesky: $A = LL^\top$, $\log|A| = 2\sum_i \log L_{ii}$
• Implicit diff through solve: $\partial_\theta x = A^{-1}(\partial_\theta b - (\partial_\theta A)\,x)$
• Jacobi’s formula: $\partial_\theta \log|A| = \mathrm{tr}(A^{-1}\partial_\theta A)$
• Jitter / safe Cholesky: $A + \epsilon I$, doubling ε until SPD
• Stable squared distances: mixed-precision $\|x-z\|^2$ to avoid catastrophic cancellation

Part 1 — Structured Linear Operators¶

1.A — Operator Zoo (catalog)¶

Key equations / models:

• Kronecker: $A\otimes B$, with Roth’s lemma $(A\otimes B)\,\mathrm{vec}(X) = \mathrm{vec}(BXA^\top)$
• BlockDiag: $\mathrm{diag}(A_1,\dots,A_K)$
• LowRankUpdate: $L + UDV^\top$
• Toeplitz: $T_{ij} = t_{i-j}$, matvec via FFT in $O(n\log n)$
• BlockTriDiag: tridiagonal precision Λ for Markov chains
• KroneckerSum: $A\oplus B = A\otimes I + I\otimes B$ (separable Laplacian)

1.B — Matrix Identities & Decompositions¶

Key equations / models:

• Kron eigendecomp: $A\otimes B = (Q_A\otimes Q_B)(\Lambda_A\otimes\Lambda_B)(Q_A\otimes Q_B)^\top$
• Sherman–Morrison–Woodbury: $(A+UCV)^{-1} = A^{-1} - A^{-1}U(C^{-1}+VA^{-1}U)^{-1}VA^{-1}$
• Det lemma: $|A+UCV| = |C^{-1}+VA^{-1}U|\,|C|\,|A|$
• Operator sandwich: $A\,P\,A^\top$ assembled lazily
• UDL of block-tridiagonal: $\Lambda = U^\top D^{-1} U$
• Discrete Lyapunov: $P_\infty = AP_\infty A^\top + Q$

1.C — Matrix-Free / Implicit¶

Key equations / models:

• Matrix-free matvec: $y = K(X,X')\,v$ via nested vmap over rows, never forming $K$
• Cross-kernel matvec: $y = K(X,Z)\,v$ for prediction without $|X|\times|Z|$ allocation

1.D — Solvers¶

Key equations / models:

• Cholesky solve: $LL^\top x = b$
• Conjugate Gradient: minimize $\tfrac{1}{2}x^\top A x - b^\top x$ in Krylov subspace
• BBMM: batched matvec drives solve + logdet + grad simultaneously (Gardner et al. 2018)
• MINRES: indefinite symmetric · LSMR: rectangular least squares
• Preconditioned CG: solve $M^{-1}Ax = M^{-1}b$ with $M\approx A$

1.E — Trace, Log-Det, Roots¶

Key equations / models:

• Hutchinson trace: $\mathrm{tr}(A) \approx \tfrac{1}{m}\sum_{i=1}^m z_i^\top A z_i$, $z_i\sim\mathrm{Rademacher}$
• Stochastic Lanczos Quadrature: $\log|A| \approx \tfrac{n}{m}\sum z_i^\top \log(A) z_i$ via Lanczos tridiagonalization
• Contour-integral root: $A^{1/2} = \tfrac{1}{2\pi i}\oint \sqrt{z}(zI-A)^{-1}dz$
• Joint inv-quad-logdet: shared CG/Lanczos passes return $y^\top A^{-1}y$ and $\log|A|$ together

Part 2 — Kernels¶

2.A — Standard kernels¶

Key equations / models:

• RBF: $k(x,x') = \sigma^2\exp(-\tfrac{1}{2\ell^2}\|x-x'\|^2)$
• Matérn-ν: $k_\nu(r) = \sigma^2 \tfrac{2^{1-\nu}}{\Gamma(\nu)}(\sqrt{2\nu}\,r/\ell)^\nu K_\nu(\sqrt{2\nu}\,r/\ell)$, $\nu\in\{1/2,3/2,5/2\}$
• Periodic: $\sigma^2\exp(-2\sin^2(\pi|x-x'|/p)/\ell^2)$
• Linear: $\sigma^2 x^\top x'$ · Polynomial: $(\sigma^2 x^\top x' + c)^d$
• ARD lengthscale: $\|x-x'\|^2_{\Lambda} = \sum_d (x_d-x'_d)^2/\ell_d^2$

2.B — Deep kernels¶

Spectral kernels (Bochner / Spectral Mixture) live in Part 7 — Spectral GPs. This subsection covers neural-network-warped kernels only.

Key equations / models:

• Deep kernel: $k_\theta(x,x') = k_{\mathrm{base}}(\phi_\theta(x), \phi_\theta(x'))$
• ArcCosine-$n$ (Cho & Saul 2009): $k_n(x,x') = \tfrac{1}{\pi}\|x\|^n\|x'\|^n J_n(\theta)$, NN-correspondence

2.C — Multi-output kernels¶

Key equations / models:

• LMC (Linear Model of Coregionalization): $f_p(x) = \sum_q a_{pq}\,g_q(x)$, $k_{pq}(x,x') = \sum_q a_{pq}a_{qq'} k_q(x,x')$
• ICM (rank-1 LMC): $K = B \otimes K_X$ with output covariance $B$
• OILMM: orthogonal projection $H$ such that $H^\top H = I$ decouples outputs

2.D — Spherical / localized kernels (moved)¶

Spherical / harmonic / Slepian kernels are now grouped under Part 7.F — Spherical / periodic spectral, where they sit alongside the full spectral toolkit (Bochner, RFF, HSGP). Riemannian / graph kernels remain in 2.F.

2.E — Kernel-based statistics & utilities¶

Key equations / models:

• Centered kernel: $\tilde K = HKH$, $H = I - \tfrac{1}{n}\mathbf{1}\mathbf{1}^\top$
• HSIC: $\tfrac{1}{n^2}\,\mathrm{tr}(K_x H K_y H)$
• MMD²: $\mathbb{E}_{x,x'}[k(x,x')] + \mathbb{E}_{y,y'}[k(y,y')] - 2\,\mathbb{E}_{x,y}[k(x,y)]$

2.F — Non-Euclidean & operator-valued kernels¶

Key equations / models:

• Operator-valued kernel: $K(x,x') \in \mathcal{L}(\mathcal{Y},\mathcal{Y})$, predicts function-valued outputs (velocity fields, spectral curves)
• Graph heat kernel: $k(u,v) = \exp(-tL)_{uv}$ for graph Laplacian $L = D - A$
• Geodesic RBF on manifold: $k(x,x') = \sigma^2\exp(-d_g(x,x')^2/2\ell^2)$, $d_g$ = geodesic distance

Part 3 — Exact GP Regression¶

3.A — Foundations¶

Key equations / models:

• GP prior: $f \sim \mathcal{GP}(0, k)$ with $y = f(x) + \epsilon$, $\epsilon \sim \mathcal{N}(0,\sigma^2)$
• Posterior mean: $m_*(x_*) = k_*^\top (K + \sigma^2 I)^{-1} y$
• Posterior variance: $v_*(x_*) = k(x_*,x_*) - k_*^\top (K + \sigma^2 I)^{-1} k_*$
• Log marginal likelihood: $\log p(y) = -\tfrac{1}{2}y^\top(K+\sigma^2 I)^{-1}y - \tfrac{1}{2}\log|K+\sigma^2 I| - \tfrac{n}{2}\log 2\pi$

3.B — Diagnostics¶

Key equations / models:

• LOO-CV (LOVE): $\mu_{-i} = y_i - [(K+\sigma^2 I)^{-1}y]_i / [(K+\sigma^2 I)^{-1}]_{ii}$
• Probability integral transform (PIT) for calibration
• Coverage at $1-\alpha$: empirical fraction of $y_i \in [\mu_i \pm z_{\alpha/2}\sigma_i]$

3.C — Heteroscedastic noise¶

Key equations / models:

• Two-GP heteroscedastic: $y = f(x) + g(x)\,\epsilon$, with $f \sim \mathcal{GP}(0,k_f)$ and $\log g \sim \mathcal{GP}(0,k_g)$
• Joint ELBO over $(f, \log g)$ via posterior linearization or coupled cubature

3.D — High-level API¶

Key equations / models:

• GPEstimator(kernel, ...).fit(X, y).predict(X*, quantiles=...) — sklearn-style facade

3.E — Constrained & Physics-informed GPs¶

Key equations / models:

• Monotone GP: derivative observations $f'(x)\ge 0$ via linear operator on prior; or projection to monotone function space
• Convex GP: Hessian positivity $\nabla^2 f(x) \succeq 0$ as inequality constraint
• Boundary-condition GP: zero mean at domain boundary $\partial\Omega$ via Dirichlet eigenfunction basis $\phi_j$ with $\phi_j\vert_{\partial\Omega}=0$
• PDE-constrained GP (Raissi et al.): encode $\mathcal{L}f = 0$ as derivative / linear-operator observations; cross-covariance $k_{\mathcal{L}}(x,x') = \mathcal{L}_x k(x,x')$
• Student-t process: $f \sim \mathcal{T}\nu(\mu, K)$; posterior analytically tractable with heavier tails than GP

Part 4 — Structured GPs¶

GPs whose covariance has direct algebraic structure (Kronecker, Toeplitz, grid, sparse-precision) — exploited by Part 1 operators.

4.A — Kronecker GPs¶

Key equations / models:

• 2D-grid GP: $K = K_x \otimes K_t$, solve in $O(n_x^3 + n_t^3)$ via Roth’s lemma
• Kronecker + low-rank: $K = K_0 + UU^\top$, posterior via Woodbury
• Sum-of-Kronecker: $K = \sum_i K_i^x \otimes K_i^t$ (additive separable)
• Additive decomposition: $f(x,t) = f_\text{trend}(t) + f_\text{seasonal}(t) + f_\text{spatial}(x,t)$

4.B — Grid / Toeplitz GPs¶

Key equations / models:

• KISS-GP / SKI: $K \approx W^\top K_U W$ with cubic local interpolation $W$ to grid points $U$
• Toeplitz matvec via FFT: $T v = \mathrm{IFFT}(\hat t \odot \mathrm{FFT}(v))$, $O(n\log n)$

4.C — Sparse-precision (mesh / GMRF)¶

Key equations / models:

• SPDE (Lindgren et al. 2011): $(\kappa^2 - \Delta)^{\alpha/2} f(x) = \mathcal{W}(x)$ → Matérn kernel
• FEM precision: $Q = \kappa^4 C + 2\kappa^2 G + GC^{-1}G$ on triangulated mesh, sparse banded

Part 5 — Approximations & Scalability¶

GPs that scale to large N via inducing points, random features, or iterative solvers — distinct from Part 4 in that the structure is imposed rather than inherent to the data geometry.

5.A — Random features (moved)¶

Random Fourier Features, Nyström, and FastFood live under Part 7.C — Random Fourier features. This subsection is intentionally empty in Part 5 — the inducing-point / variational story (5.B–5.E) does not depend on RFF derivations, only on the spectral approximation result they deliver.

5.B — Inducing-point fundamentals¶

Key equations / models:

• Inducing variables: $u = f(Z)$ with $Z \in \mathbb{R}^{m\times d}$; FITC: $K \approx Q_{nn} + \mathrm{diag}(K_{nn} - Q_{nn})$, $Q_{nn} = K_{nm}K_{mm}^{-1}K_{mn}$
• SVGP ELBO (Hensman 2013): $\mathcal{L} = \sum_i \mathbb{E}_{q(f_i)}[\log p(y_i\mid f_i)] - \mathrm{KL}(q(u)\Vert p(u))$
• Whitened: $u = L_{mm}\tilde u$, $\tilde u \sim \mathcal{N}(0,I)$ → isotropic optimization
• Collapsed ELBO (Titsias 2009): $\mathcal{L} = \log\mathcal{N}(y;0, Q_{nn}+\sigma^2 I) - \tfrac{1}{2\sigma^2}\mathrm{tr}(K_{nn}-Q_{nn})$

5.C — Inter-domain features¶

Inter-domain features that are fundamentally spectral (VFF, VISH, Laplacian eigenfunctions) live under Part 7.E — Variational spectral methods and Part 7.F — Spherical / periodic spectral. This subsection now covers only the generic decoupled-basis pattern.

Key equations / models:

• Inter-domain inducing variables: $u_j = \langle f, g_j \rangle_\mathcal{H}$ for chosen basis $\{g_j\}$
• Decoupled: separate basis for posterior mean (large) and covariance (small)

5.D — Iterative-solver scaling¶

Key equations / models:

• CG-based GP: $\alpha = (K+\sigma^2 I)^{-1} y$ via CG; logdet via SLQ
• CGLB: $\log|A| \geq -\sum \log(1-\theta_k)$ from Lanczos eigenvalue bounds
• Preconditioned CG: $M^{-1}Ax = M^{-1}b$ with Nyström / pivoted-Cholesky preconditioner
• EigenPro: SGD with eigen-spectrum preconditioner $\Lambda^{-1}$
• Falkon (Rudi 2017): Newton iteration on Nyström-reduced system, $O(N\sqrt{N}\log(1/\epsilon))$
• LogFalkon: extends Falkon to GSC losses (logistic, exponential)

5.E — Deep GPs¶

Key equations / models:

• Deep GP (Salimbeni & Deisenroth 2017): $f^{(L)} \circ \cdots \circ f^{(1)}$, each layer $f^{(l)}\sim\mathcal{GP}$
• Doubly stochastic ELBO: $\mathcal{L} = \sum_n\mathbb{E}_{q(f^{1:L}_n)}[\log p(y_n\mid f_n^{(L)})] - \sum_l\mathrm{KL}(q(u^{(l)})\Vert p(u^{(l)}))$
• Convolutional GP (van der Wilk 2017): patch-level inducing features, $f(x) = \sum_p f_p(x_p)$, $f_p \sim \mathcal{GP}$ over patches
• Inter-domain inducing variables per patch; $K_{uu}$ exploits translation-equivariance

Part 6 — Non-Conjugate Likelihoods & Inference¶

6.A — Likelihood & integrator zoos¶

Key equations / models:

• Generic non-conjugate factorization: $p(y\mid f) = \prod_i p(y_i\mid f_i)$
• Bernoulli ($\sigma(f)$), Poisson ($\exp(f)$), Student-t, Beta, Gamma, Exponential, Softmax
• ELL: $\mathbb{E}_{q(f)}[\log p(y\mid f)]$ via integrator
• Gauss–Hermite: $\int g(z)\mathcal{N}(z;\mu,\sigma^2)\,dz \approx \sum_i w_i\, g(\mu+\sqrt 2\sigma\,\xi_i)$
• Sigma points (Unscented): $2d+1$ deterministic points
• 5th-order cubature: symmetric $2d^2+1$ points

6.B — Classification¶

Key equations / models:

• Latent GP classification: $f\sim\mathcal{GP}$, $y\sim\mathrm{Bernoulli}(\sigma(f))$
• Multi-class: $y\sim\mathrm{Cat}(\mathrm{softmax}(f_1,\dots,f_K))$, $K$ latent GPs

6.C — Newton / Gauss-Newton family¶

Key equations / models:

• Laplace: $q(f) = \mathcal{N}(\hat f, -H^{-1})$, $\hat f = \arg\max p(f\mid y)$, $H = \nabla^2\log p(f\mid y)$
• Newton update: $f \leftarrow f - H^{-1}\nabla$ with damping
• Gauss-Newton / GGN: $H \approx J^\top R J$ (drops 2nd-order terms)
• Posterior linearization (SLR): $p(y\mid f) \approx \mathcal{N}(y; Af + b, \Omega)$ matched on Gaussian
• Hutchinson Hessian diag: $\mathrm{diag}(H) \approx \mathbb{E}[z\odot Hz]$

6.D — Variational inference¶

Key equations / models:

• Variational guides: delta · diagonal mean-field · low-rank ($S = VV^\top + \mathrm{diag}$) · full-rank Cholesky · normalizing flow · whitened
• ELBO: $\log p(y) \geq \mathbb{E}_q[\log p(y,f)] - \mathbb{E}_q[\log q(f)]$
• Natural gradient: $\tilde\nabla \mathcal{L} = F^{-1}\nabla \mathcal{L}$, $F$ = Fisher
• CVI sites (Khan & Lin 2017): $\eta^{(t+1)} = (1-\rho)\eta^{(t)} + \rho\,\hat\eta_{\mathrm{tilted}}$

6.E — Expectation Propagation¶

Key equations / models:

• Cavity: $q^{-i}(f_i) \propto q(f_i) / t_i(f_i)$
• Tilted: $\hat q(f_i) \propto q^{-i}(f_i)\, p(y_i\mid f_i)$
• Site update: choose $t_i$ such that $\mathrm{KL}(\hat q\,\Vert\, q^{-i} t_i)$ is minimized → moment match

6.F — Bayesian linear regression & non-standard outputs¶

Key equations / models:

• BLR posterior: $\Sigma = (\Phi^\top R^{-1}\Phi + S_0^{-1})^{-1}$, $\mu = \Sigma\,\Phi^\top R^{-1}y$
• Sequential update via Sherman–Morrison: rank-1 covariance update on each new observation
• Log-Gaussian Cox Process: $\lambda(x) = \exp(f(x))$, observations from Poisson process
• Warped GP (Snelson 2003): $g(y) = f(x)$ with monotone bijection $g$, transformed likelihood

6.G — Aggregate Bayesian methods¶

Key equations / models:

• INLA (Rue et al. 2009): $p(\theta\mid y) \approx p(y\mid x^*,\theta)p(x^*\mid\theta)p(\theta)/q(x^*\mid y,\theta)$ at Laplace mode $x^*$
• Numerical integration over θ on grid, marginal posteriors of latent field

Part 7 — Spectral GPs¶

Spectral methods unify several threads from earlier parts: stationary kernels via Bochner’s theorem (extending 2.A), spectral kernel families (former 2.B), random-feature approximations (former 5.A), and spectral inducing-feature methods (subset of former 5.C). This part collects them with a coherent through-line: every stationary kernel is the Fourier transform of a positive measure, and every approximation in this part is a clever way of sampling, parameterizing, or projecting that measure.

7.A — Spectral foundations¶

Key equations / models:

• Bochner’s theorem: stationary $k(\tau) = \int_{\mathbb{R}^d} e^{i\omega^\top\tau}\,S(\omega)\,d\omega$ with $S \geq 0$ a finite measure (the power spectral density).
• Wiener–Khinchin: covariance $\leftrightarrow$ spectral density via Fourier transform; Toeplitz–circulant duality on a regular grid (pairs with 1.4).
• Karhunen–Loève expansion: $f(x) = \sum_j \sqrt{\lambda_j}\,\xi_j\,\phi_j(x)$, $\xi_j \sim \mathcal{N}(0,1)$, $(\lambda_j, \phi_j)$ eigenpairs of the kernel integral operator.
• Mercer’s theorem: $k(x,x') = \sum_j \lambda_j \phi_j(x)\phi_j(x')$, basis for the spectral representations in 7.D and 7.E.
• Sampling theorem connection: bandlimited prior + grid spacing $\Delta x < \pi/\omega_{\max}$ → no aliasing in Toeplitz / FFT matvec.

7.B — Spectral kernel models¶

Key equations / models:

• Spectral mixture (Wilson & Adams 2013): $S(\omega) = \sum_q w_q\,\mathcal{N}(\omega; \mu_q, \Sigma_q)$ — closed-form $k(\tau)$ via inverse FT.
• Spectral mixture product (multi-D): outer product of 1-D SM along each axis.
• Sparse spectrum kernel (Lázaro-Gredilla 2010): finite mixture of point masses in $S(\omega)$.

7.C — Random Fourier features¶

Key equations / models:

• RFF (Rahimi & Recht 2007): $\phi(x) = \sqrt{2/D}\cos(\omega^\top x + b)$, $\omega\sim S(\omega)$, $k(x,x') \approx \phi(x)^\top\phi(x')$.
• SSGP / VSSGP: Bayesian linear regression in RFF space; hierarchical / variational priors over ω.
• Nyström: $K \approx K_{nm}K_{mm}^{-1}K_{mn}$, rank-$m$ approximation; data-dependent counterpart to RFF.
• FastFood: $W = SHG\Pi HB$ — structured frequency matrix, $O(D\log d)$ matvec.
• Orthogonal random features: $\omega_i$ drawn from a uniform measure on the sphere (variance reduction).

7.D — Hilbert-space methods¶

Key equations / models:

• Hilbert-space GP (Solin & Särkkä 2020): on a bounded domain Ω with Dirichlet boundary, Laplacian eigenpairs $(\lambda_j, \phi_j)$ give $k(x,x') \approx \sum_{j=1}^M S(\sqrt{\lambda_j})\,\phi_j(x)\phi_j(x')$ — diagonalized prior, $O(NM + M^3)$ inference.
• Periodic kernel via truncated Fourier series: orthonormal basis on $[0, 2\pi]$.
• Convergence of HSGP to exact GP as $M \to \infty$ for stationary kernels (Mercer rate).

7.E — Variational spectral methods¶

Key equations / models:

• Variational Fourier Features (Hensman 2018): $u_j = \int f(x)\cos(j\pi x/L)\,dx$ on $[-L,L]$ → diagonal $K_{uu}$, $O(NM)$.
• Inter-domain spectral inducing variables: choose basis $\{g_j\}$ to be eigenfunctions of the kernel operator → diagonalized $K_{uu}$.
• Variational Sparse Spectrum GP (VSSGP, Gal & Turner 2015): variational distribution over RFF frequencies.

7.F — Spherical / periodic spectral¶

Key equations / models:

• Spherical harmonics on $S^2$: $Y_{\ell m}$ orthonormal eigenfunctions of the Laplacian; zonal kernels $k(x,x') = \sum_\ell a_\ell P_\ell(x^\top x')$.
• Funk–Hecke: convolution of zonal function with spherical-harmonic basis is diagonal in $\ell$.
• Variational Inducing Spherical Harmonics (VISH, Dutordoir 2020): inducing variables = spherical-harmonic projections; diagonal $K_{uu}$.
• Spherical Slepian: solve $\int_R Y^*_{\ell m}Y_{\ell'm'}\,d\Omega \cdot c = \lambda c$ to maximize energy in region $R$.
• Fourier features on the torus / circle: periodic boundary conditions → exact spectral kernel.

Part 8 — Markov / State-Space GPs¶

8.A — Foundations¶

Key equations / models:

• Discrete SSM: $x_{k+1} = A_k x_k + w_k$, $y_k = H_k x_k + v_k$, $w_k\sim\mathcal{N}(0,Q_k)$, $v_k\sim\mathcal{N}(0,R_k)$
• Kalman predict: $\bar x = A x$, $\bar P = A P A^\top + Q$
• Kalman update: $S = H\bar P H^\top + R$, $K = \bar P H^\top S^{-1}$, $x^+ = \bar x + K(y - H\bar x)$
• RTS smoother: backward $G_k = P_k A_k^\top \bar P_{k+1}^{-1}$
• Joseph form: $P^+ = (I-KH)\bar P(I-KH)^\top + KRK^\top$ (numerically stable)

8.B — SDE kernel zoo¶

Key equations / models:

• LTI SDE: $dx = Fx\,dt + L\,dW$, observation $f(t) = Hx(t)$
• Kernel ↔ SDE map (Hartikainen & Särkkä 2010): Matérn-3/2 → 2D LTI; Matérn-5/2 → 3D
• Discretization: $A_k = \exp(F\Delta t_k)$, $Q_k = P_\infty - A_k P_\infty A_k^\top$ via Lyapunov
• Periodic: truncated Fourier, block-diagonal $F$ · QuasiPeriodic: Matérn × Periodic via $\oplus$
• Drift-KL: $\mathrm{KL}(p\Vert q) = \int_0^T \tfrac{1}{2}\|f_p - f_q\|^2_{\Sigma^{-1}}\,dt$

8.C — Markov GP workflows¶

Key equations / models:

• Marginal log-likelihood: $\log p(y) = \sum_k \log\mathcal{N}(y_k; H\bar x_k, S_k)$ from filter pass
• Hyperparameter learning: gradient through Kalman filter
• Sparse variational Markov GP: ELBO via filter over inducing time points $u_{1:M}$
• KalmanGuide: pseudo-observations $\tilde y_k$, $\tilde R_k$ → standard Kalman + RTS for posterior

8.D — Parallel & scalable filtering¶

Key equations / models:

• Parallel scan (Särkkä & García-Fernández 2021): associative op $\otimes$ on $(A_k, b_k)$ pairs, $O(\log N)$ depth
• Square-root form: propagate $\sqrt P$ instead of $P$ for numerical stability
• SpInGP: parallel-in-time + sparse (banded) state representation
• Mean-field Kalman: block-diagonal $P$ across independent dimensions

8.E — Nonlinear filtering¶

Key equations / models:

• EKF: linearize $f, h$ around $\hat x$ via Jacobian
• UKF: propagate $2d+1$ sigma points through $f$, recompute moments
• CKF: $2d$ cubature points, no tuning parameter
• Innovation cov as LowRankUpdate when $R$ has structure: $S = H\bar P H^\top + R$

8.F — Ensemble methods¶

Key equations / models:

• EnKF analysis: $X^a = X^f + P_e H^\top (HP_eH^\top + R)^{-1}(Y - HX^f)$, $P_e = \tfrac{1}{J-1}(X^f - \bar X^f)(X^f - \bar X^f)^\top$
• Bessel-corrected covariance: divide by $J-1$ not $J$
• Ensemble Kalman gain: low-rank $K_e = X^f X^{f\top}H^\top S^{-1}$

8.G — Steady-state & structured-Gaussian surfaces¶

Key equations / models:

• DARE: $P = APA^\top + Q - APH^\top(HPH^\top+R)^{-1}HPA^\top$ → unique SPD solution
• Infinite-horizon Kalman: solve DARE once, reuse $K_\infty$
• MarkovGaussian surface: structured Gaussian with $(A_k, Q_k)$ as PyTree, exposes filter/smoother/sample

8.H — Non-conjugate temporal case studies¶

Key equations / models:

• Laplace + Kalman: iterate site moments $(\mu_k, \tau_k)$ via Newton on each $p(y_k\mid f_k)$, run filter+smoother per iteration
• EP + Kalman: same loop with EP cavity / tilted moments
• Changepoints (additive): $f(t) = f_\text{trend}(t) + f_\text{fast}(t)$, Matérn-5/2 + Matérn-1/2
• Streaming filter-only: discard $x_{k<N-W}$ for fixed window $W$
• Time-varying lengthscale: $\ell(t)$ as random walk in numpyro.scan

Part 9 — Sampling, Pathwise, Conditioning¶

9.A — Pathwise sampling¶

Key equations / models:

• Pathwise (Wilson 2020): $f^*(x) = \underbrace{\sum_i w_i\phi_i(x)}_\text{prior basis} + \underbrace{k(x,X)\beta}_\text{update}$ with $\beta = (K+\sigma^2 I)^{-1}(y - \Phi w - \epsilon)$
• Decoupled SVGP sampling: parametric prior + non-parametric update from inducing variables only

9.B — Matheron’s-rule conditioning¶

Key equations / models:

• Matheron’s rule: $f\mid y \stackrel{d}{=} f + K_*(K+\sigma^2 I)^{-1}(y - f(X) - \epsilon)$ where $f, \epsilon$ are independent prior samples
• Partitioned joint: factor $p(f_a, f_b\mid y)$ as $p(f_a\mid y)\,p(f_b\mid f_a, y)$ with shared Cholesky