[PyMCon Web Series 03a] Introduction to Hilbert Space GPs (HSGPs) in PyMC (Mar 15, 2023) (Bill Engels)

[bwengals]

@jiajie_wu

Hello! May I ask when HSGP will be officially launched? Are there any learning examples available now?

Not sure, it’ll be part of the next release. As of now, there’s not other PyMC HSGP material out there – working on some for the next event and for the docs. If you want to give it a try sooner, the examples in the docstrings of HSGP.prior, HSGP.conditional and HSGP.prior_linearized

@ScroobiusTrip

Thank you for the great talk! I was wondering, is it possible to tack on an arbitrary mean function onto a HSGP in the same way that you would for a vanilla GP? Thinking of your example here , where you add a piecewise linear mean onto the existing GP.

@ScroobiusTrip I wish I’d mentioned this in the talk. You don’t really need to worry about whether something is specifically a mean function if you’re using HSGPs. Just add the (mean zero) HSGP to the piecewise linear part of your model – this is equivalent to using the piecewise linear model as the mean function within a GP. For the rest of the GP implementations in PyMC, you need to at least wrap things into a Mean function object to keep prediction via gp.conditional works correctly.

For your case though, just add them! If you use HSGP.prior_linearized you can use pm.set_data and it’ll work for both parts of the model.

@murphyk

Question for Bill.
How fast is it to generate the basis set given X and kernel hparams?
Can this conversion be released as a standalone funnction outside of pymc?
If you want to optimize the hparams, can you just optimize the marginal likelihood
of the corresponding induced linear model, ignoring the GP?

@murphyk

• It’s \mathcal{O}(m^2 n) to generate the basis (m basis functions, n data points). But, you don’t need to know anything about the kernel hyperparams or even the kernel to do it.
• No plans now, but the implementation of the tricky parts is here.
• I think so! f = phi @ (beta * sqrt_psd), with each \beta_i \sim N(0, 1). So you should be able to calculate the marginal likelihood integrating out beta just like one could with a vanilla Bayesian linear regression model. Sorry, I didn’t quite get this live.

@michaelosthege

Can I sample >one function< from the posterior and then deterministically evaluate that function on new X ?

…to run an optimizer on that function

…it’s not possible with GPs, but my prior is that this approximation allows it (like the RFF approximation which is very similar)

@michaelosthege
Does a function like this do the trick?

import pymc as pm
import pytensor.tensor as pt

def one_gp_function(X_new, beta_posterior_sample, cov_func, L, m, X_mean):
    Xs = X_new - X_mean
    L = pt.as_tensor_variable(L)
    eigvals = pm.gp.hsgp_approx.calc_eigenvalues(L, m, tl=pt)
    phi = pm.gp.hsgp_approx.calc_eigenvectors(Xs, L, eigvals, m, tl=pt)

    omega = pt.sqrt(eigvals)
    psd = cov_func.power_spectral_density(omega).eval()
    f = phi @ (beta_posterior_sample * pt.sqrt(psd))
    return f

I’m sure this could be cleaner, but I think one “function” from the GP posterior will correspond to one draw from the posterior of beta. Then you can regenerate the phi basis given arbitrary X and that’s your function.