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

I recently found that using numpyroc sampling in pymc will use much more memory than pure numpyro. So I had to use numpyro instead. They also have an example of HSGP. Itā€™s not clean and easy to use, but itā€™s better than the limited memory I have (8GB).

Oh, thatā€™s a bummer. You could also try nutpie, blackjax or just the built-in PyMC sampler too.

2 Likes

@bwengals has now posted an HGSP ā€œstarterā€ notebook. The notebook uses the well-known cherry blossom data and walks through a variety of approaches, culminating in the Hibert space Gaussian processes he discussed in his live talk (soon to be posted to the PyMC YouTube channel).

The notebook can be found here. If you try it out, either with the cherry blossom data or with your own data, feel free to post comments or questions about how it worked out for you.

4 Likes

@bwengals I just saw the video of your talk and your HSGP approach looks very promising to handle (my) common performance problems with GPs! :slightly_smiling_face:

However, there were just two questions coming into my mind:

  1. Is HSGP also working with WrapedInput kernels? Iā€™m frequently using WrapedInput to model GPs on non-euclidean domains, e.g. data on a 2-D sphere.

  2. Iā€™m not an expert on Hilbert space methods, but if I understood your approach correctly, HSGP is basically approximating the GP via a Fourier series with random coefficients. Then, it should be quite natural to extend this method to use any kind of orthonormal basis functions for approximation, e.g. orthogonal polynomials, wavelets, etc. - shouldnā€™t it? Intuitively, I would have maybe chosen Hermite polynomials to approximate GPs on an euclidean domain, to avoid the artificial periodicity that seems to be currently introduced if you predict values outside [-L, L].

I hope it helps!

  1. Is HSGP also working with WrapedInput kernels? Iā€™m frequently using WrapedInput to model GPs on non-euclidean domains, e.g. data on a 2-D sphere.

Yes and no. You can absolutely transform your X inputs, but youā€™ll have to do it outside of the WarpedInput covariance kernel unfortunately. HSGPs can only be applied to stationary kernels with a power spectral density. So for now, ExpQuad, Matern52 and Matern32.

  1. Iā€™m not an expert on Hilbert space methods, but if I understood your approach correctly, HSGP is basically approximating the GP via a Fourier series with random coefficients. Then, it should be quite natural to extend this method to use any kind of orthonormal basis functions for approximation, e.g. orthogonal polynomials, wavelets, etc. - shouldnā€™t it? Intuitively, I would have maybe chosen Hermite polynomials to approximate GPs on an euclidean domain, to avoid the artificial periodicity that seems to be currently introduced if you predict values outside [āˆ’L,L].

Itā€™s certainly not my method! I merely ported it into PyMC. It was first proposed in this paper. I donā€™t think you can trivially extend the method to use any type of orthonormal basis functions because it depends on the expression of a kernel as a power spectral density.