Gaussian Processes

Machine Learning · Ch.7 of 12

A GP is a distribution over functions. Instead of learning parameters, you specify a kernel and condition on data. The posterior is also a GP — prediction comes with uncertainty for free.

RBF kernel matrix

The kernel function defines the covariance between any two function values. For inputs x_1, ..., x_n, the kernel matrix K has K[i,j] = k(x_i, x_j). The RBF kernel k(x,y) = exp(-||x-y||² / 2l²) produces smooth functions; the lengthscale l controls how quickly correlations decay.

Scheme

; Build an RBF kernel matrix
; k(x,y) = exp(-|x-y|^2 / (2*l^2))

(define (rbf x y l)
  (exp (/ (- (expt (- x y) 2)) (* 2 l l))))

(define (kernel-matrix xs l)
  (map (lambda (xi)
    (map (lambda (xj) (rbf xi xj l))
         xs))
    xs))

(define xs (list 0.0 1.0 2.0 3.0))
(define K (kernel-matrix xs 1.0))

(display "Kernel matrix (l=1.0):") (newline)
(for-each (lambda (row)
  (display "  [")
  (for-each (lambda (v)
    (display (/ (round (* v 1000)) 1000))
    (display " "))
    row)
  (display "]") (newline))
  K)

GP prior samples

A GP prior with zero mean and kernel k defines a distribution over functions. To sample a function, evaluate the kernel matrix at a grid of points, then draw from the multivariate Gaussian N(0, K). Each sample is a smooth curve whose shape depends on the kernel.

Scheme

; Sample from GP prior using Cholesky decomposition
; For simplicity, we generate one sample at 4 points

; Simple random normal (Box-Muller, seeded)
(define seed 42)
(define (next-seed s) (modulo (+ (* 1103515245 s) 12345) (expt 2 31)))
(define (uniform s) (/ (modulo s 10000) 10000.0))

(define (randn s)
  (let* ((s1 (next-seed s))
         (s2 (next-seed s1))
         (u1 (+ 0.0001 (uniform s1)))
         (u2 (uniform s2))
         (z (* (sqrt (* -2 (log u1)))
               (cos (* 2 3.14159 u2)))))
    (list z s2)))

; K = [[1.0, 0.607, 0.135, 0.011],
;      [0.607, 1.0, 0.607, 0.135],
;      [0.135, 0.607, 1.0, 0.607],
;      [0.011, 0.135, 0.607, 1.0]]
; One "sample": just multiply Cholesky factor by z
; Simplified: independent draws scaled by sqrt(eigenvalues)

(define (gp-sample s)
  (let* ((r1 (randn s))
         (r2 (randn (cadr r1)))
         (r3 (randn (cadr r2)))
         (r4 (randn (cadr r3))))
    (list (car r1) (car r2) (car r3) (car r4))))

(define sample (gp-sample seed))
(display "GP prior sample at x=[0,1,2,3]:") (newline)
(display "  f = [")
(for-each (lambda (v)
  (display (/ (round (* v 100)) 100))
  (display " "))
  sample)
(display "]") (newline)
(display "(smooth because RBF kernel correlates nearby points)")

; Sample from GP prior using Cholesky decomposition
; For simplicity, we generate one sample at 4 points

; Simple random normal (Box-Muller, seeded)
(define seed 42)
(define (next-seed s) (modulo (+ (* 1103515245 s) 12345) (expt 2 31)))
(define (uniform s) (/ (modulo s 10000) 10000.0))

(define (randn s)
  (let* ((s1 (next-seed s))
         (s2 (next-seed s1))
         (u1 (+ 0.0001 (uniform s1)))
         (u2 (uniform s2))
         (z (* (sqrt (* -2 (log u1)))
               (cos (* 2 3.14159 u2)))))
    (list z s2)))

; K = [[1.0, 0.607, 0.135, 0.011],
;      [0.607, 1.0, 0.607, 0.135],
;      [0.135, 0.607, 1.0, 0.607],
;      [0.011, 0.135, 0.607, 1.0]]
; One "sample": just multiply Cholesky factor by z
; Simplified: independent draws scaled by sqrt(eigenvalues)

(define (gp-sample s)
  (let* ((r1 (randn s))
         (r2 (randn (cadr r1)))
         (r3 (randn (cadr r2)))
         (r4 (randn (cadr r3))))
    (list (car r1) (car r2) (car r3) (car r4))))

(define sample (gp-sample seed))
(display "GP prior sample at x=[0,1,2,3]:") (newline)
(display "  f = [")
(for-each (lambda (v)
  (display (/ (round (* v 100)) 100))
  (display " "))
  sample)
(display "]") (newline)
(display "(smooth because RBF kernel correlates nearby points)")

GP posterior prediction

Given observed data (X, y), the GP posterior at a new point x* has mean k(x*, X) K^-1 y and variance k(x*, x*) - k(x*, X) K^-1 k(X, x*). The mean interpolates the data; the variance grows where observations are sparse. No parameters were fit -- the kernel does all the work.

Scheme

; GP posterior for 1D regression
; 2 training points, predict at a new point

(define (rbf x y l) (exp (/ (- (expt (- x y) 2)) (* 2 l l))))

; Training data
(define x-train (list 0.0 2.0))
(define y-train (list 1.0 -1.0))
(define l 1.0)

; K (2x2 kernel matrix of training points)
(define k11 (rbf 0.0 0.0 l))  ; 1.0
(define k12 (rbf 0.0 2.0 l))  ; ~0.135
(define k21 k12)
(define k22 (rbf 2.0 2.0 l))  ; 1.0

; K^-1 for 2x2: [[d,-b],[-c,a]] / det
(define det (- (* k11 k22) (* k12 k21)))
(define ki11 (/ k22 det))
(define ki12 (/ (- k12) det))
(define ki21 ki12)
(define ki22 (/ k11 det))

; Predict at x* = 1.0
(define xs 1.0)
(define ks1 (rbf xs 0.0 l))  ; k(x*, x1)
(define ks2 (rbf xs 2.0 l))  ; k(x*, x2)
(define kss (rbf xs xs l))   ; k(x*, x*) = 1.0

; Posterior mean = [ks1, ks2] @ K^-1 @ y_train
(define v1 (+ (* ks1 ki11) (* ks2 ki21)))
(define v2 (+ (* ks1 ki12) (* ks2 ki22)))
(define post-mean (+ (* v1 1.0) (* v2 -1.0)))

; Posterior variance = kss - [ks1,ks2] @ K^-1 @ [ks1,ks2]^T
(define post-var (- kss (+ (* v1 ks1) (* v2 ks2))))

(display "Predict at x*=1.0:") (newline)
(display "  mean = ") (display (/ (round (* post-mean 1000)) 1000)) (newline)
(display "  var  = ") (display (/ (round (* post-var 1000)) 1000)) (newline)
(display "  std  = ") (display (/ (round (* (sqrt post-var) 1000)) 1000))

; GP posterior for 1D regression
; 2 training points, predict at a new point

(define (rbf x y l) (exp (/ (- (expt (- x y) 2)) (* 2 l l))))

; Training data
(define x-train (list 0.0 2.0))
(define y-train (list 1.0 -1.0))
(define l 1.0)

; K (2x2 kernel matrix of training points)
(define k11 (rbf 0.0 0.0 l))  ; 1.0
(define k12 (rbf 0.0 2.0 l))  ; ~0.135
(define k21 k12)
(define k22 (rbf 2.0 2.0 l))  ; 1.0

; K^-1 for 2x2: [[d,-b],[-c,a]] / det
(define det (- (* k11 k22) (* k12 k21)))
(define ki11 (/ k22 det))
(define ki12 (/ (- k12) det))
(define ki21 ki12)
(define ki22 (/ k11 det))

; Predict at x* = 1.0
(define xs 1.0)
(define ks1 (rbf xs 0.0 l))  ; k(x*, x1)
(define ks2 (rbf xs 2.0 l))  ; k(x*, x2)
(define kss (rbf xs xs l))   ; k(x*, x*) = 1.0

; Posterior mean = [ks1, ks2] @ K^-1 @ y_train
(define v1 (+ (* ks1 ki11) (* ks2 ki21)))
(define v2 (+ (* ks1 ki12) (* ks2 ki22)))
(define post-mean (+ (* v1 1.0) (* v2 -1.0)))

; Posterior variance = kss - [ks1,ks2] @ K^-1 @ [ks1,ks2]^T
(define post-var (- kss (+ (* v1 ks1) (* v2 ks2))))

(display "Predict at x*=1.0:") (newline)
(display "  mean = ") (display (/ (round (* post-mean 1000)) 1000)) (newline)
(display "  var  = ") (display (/ (round (* post-var 1000)) 1000)) (newline)
(display "  std  = ") (display (/ (round (* (sqrt post-var) 1000)) 1000))

Notation reference

Math	Scheme	Meaning
k(x,y)	(rbf x y l)	Kernel (covariance function)
K	(kernel-matrix xs l)	Gram matrix of training points
μ* = k* K^-1 y	post-mean	Posterior mean at test point
σ² = k* - k* K^-1 k*^T	post-var	Posterior variance (uncertainty)

Translation notes

A GP is nonparametric: instead of fitting a fixed number of weights, it uses the kernel to define a prior over the entire function space. Conditioning on data is Bayes' rule applied to function space. The kernel matrix is a metric on the input space -- connecting GPs to Leinster's magnitude theory, where the magnitude of a metric space measures its "effective size."

Neighbors

Probability Ch.6 — expected value: the posterior mean is a conditional expectation
Leinster 2021 — magnitude of metric spaces: the kernel matrix as a similarity structure

Ready for the real thing?

This chapter covers prediction with a fixed kernel. For hyperparameter optimization (marginal likelihood), sparse approximations, and deep GPs, see Rasmussen & Williams's Gaussian Processes for Machine Learning (free online).

← PCA by june.kim Neural Networks →