CODEWITHSUNDEEP
pymc3
swmfg
swmfg

Reputation: 1359

Solving Gaussian Process posterior analytically in PyMC3

I'm used to doing a Gaussian Process regression in GPFlow which lets you do this to solve for the posterior analytically:

import gpflow as gp
from gpflow.kernels import RBF, White, Periodic, Linear
k = RBF(x.shape[1]) + White(x.shape[1])
m = gp.models.GPR(x, y, k)
self.model = m
m.compile()
opt = gp.train.ScipyOptimizer()
opt.minimize(m)

I've recently moved to PyMC3 and am trying to accomplish the same thing as above. I've found in the documentation this bit of code (https://docs.pymc.io/notebooks/GP-slice-sampling.html#Examine-actual-posterior-distribution):

# Analytically compute posterior mean
L = np.linalg.cholesky(K_noise.eval())
alpha = np.linalg.solve(L.T, np.linalg.solve(L, f))
post_mean = np.dot(K_s.T.eval(), alpha)

Ultimately I want to be using the GP for a regression on unseen data. Is using np.linalg the right way to analytically solve a GP posterior?

Upvotes: 1

Views: 279

Answers (1)

merv
merv

Reputation: 76780

Sure. As stated in the tutorial, they implemented Algorithm 2.1 from Rasmussen's GPML, and he explicitly uses left matrix divide notation (\), which indicates to use a linear solve. For example, in theory (i.e., real number system),

A\b === A^(-1) * b === x

where x solves A*x = b. But in a practical computational domain (e.g., IEEE floating point), this equivalence breaks down because solve(A, b) is faster and more numerically stable than inv(A) * b.

The left matrix divide (\) notation is commonplace in Numerical Linear Algebra, and I'd venture that the most salient explanation for its preference is that it tacitly reminds students to never compute a matrix inverse, unnecessarily.

Upvotes: 3

Related Questions