Comment by bhl | Hacker News Reader

Comment by bhl | Hacker News Reader

Implementing Gaussian Processes in Python is a fun way to practice your programming and mathematical abilities.

I was working on a course project that used GPs for regression, with the inputs taken from the output of a CNN. We were trying to predict crop yield from satellite imagery, and GPs helped to incorporate spatio-temporal locality: think of predictions for neighboring county, or for image from a previous year.

One issue we ran into while using the model from the original authors [1] was a memory constraint from building on top of Google Colab instead of a private AWS instance, and using a larger dataset to predict upon.

The Gaussian kernel matrix grows O(n^2) per data point, and we had up to 22,000 datapoints to use. To get around the memory constraint, we had to re-implement the GP from scratch - carefully avoiding memory spikes that came with computing the kernel and mean predictions.

It was a two-fold process. For understanding the exact equations, I would highly recommend the resource mentioned in the article [2], Gaussian Processes for Machine Learning Chapter 2. For a lecture video, I would recommend [3] which is from a ML course at Cornell.

As for implementing GPs efficiently, this is still a challenge: numpy has a bunch of matrix-vector functions that are helpful, but there isn't a way to exploit the symmetry of the kernel matrix and use half the require memory.

[1] https://aaai.org/ocs/index.php/AAAI/AAAI17/paper/view/14435/...

[2] http://www.gaussianprocess.org/gpml/chapters/RW.pdf

[3] https://www.youtube.com/watch?v=R-NUdqxKjos

Another easy solution, if your kernel is appropriate (goes to 0 at infinity such as a Gaussian, Matèrn or Exponential kernel), is to set a threshold under which your covariance is set to zero and to use a sparse representation for the covariance matrix (you can even be more extreme and set a fixed number of entry to be non-zero giving you a fixed memory use).

This does not solve the issue. You still need to invert this sparse matrix for inference. The inverse of a sparse covariance matrix is dense.

The raw mathematical formula uses an inversion but what you truly need is to solve a linear system (which is why most implementations store a cholesky decomposition of the covariance matrix) which can be done without transforming the matrix into a dense matrix.

What would you use to solve the system? Iterative methods or gradient descent even?

yes

Something like MINRES or GMRES or CGSTAB. Correct me if I'm wrong, but the covariance matrix in GPR should be diagonally dominant, so the diagonal of the matrix can be used as a somewhat effective preconditioner.

With most kernels the matrix is symmetric and positive-definite so I would probably use a conjugate gradient: https://en.wikipedia.org/wiki/Conjugate_gradient_method

+1 This is how I solved for the optimal parameter: instead of inverting the matrix, which is memory-heavy and computationally expensive, I used the solve function from scipy.linalg.

I looked into numpy and scipy methods for symmetry and sparsity, but my kernel matrix is still too dense.

Regarding the challenge of an efficient implementation, pytorch+pyro attempt to tackle it. It's still a work in progress, but you might want to have a look at: https://pyro.ai/examples/gp.html

The real issue is that inversion / solving is O(n^3).

This really hammers you for high dimensional settings (e.g. brain imaging n>100k).

Solutions:

* Pick a very special kernel with an known / approximable solution.

* Reduce the dimension of the problem (sparse gps).

The n parameter there was the number of datapoints, not the size of the inputs (that would've been way worse). For brain imaging, is the 100k parameter also referencing datapoints, or dimensionality of the input?

Chapter 8 of the Rasmussen reference does talk about approximation methods for large datasets, but I couldn't find any reference implementations to go off of.

"number of datapoints, not the size of the inputs"

For brain imaging, 100k was 'the dimensionality of the input' (but also the number of datapoints). This could easily be O(10^6) though with higher res imaging.

I'd check out, for discussion:

https://www.prowler.io/blog/sparse-gps-approximate-the-poste...

For code:

https://github.com/GPflow/GPflow

Based on the people involved knowing what they are talking about, rather than on experience using that particular work.

Seeing the complications implied by the O(n^2) problem, I'm curious why it was worth it? What did using GPs to solve the problem buy you that you couldn't get with other methods?

I think GPs are interesting but the kernel matrix is a real problem with them in practice. I have yet to understand why they are worth the cost compared to just using least squares or some other regression method. Confidence intervals are an obvious benefit, but since you can get them via monte-carlo, at what point does the GP method pay off?

Compared to training the CNN, running a GP on the outputs seemed like a worthwhile tradeoff of time if it resulted in a discernible difference in validation error; and the original authors of the paper we were following already had an open implementation.

I am open to suggestions if you know of any "post-processing" methods that would help to incorporate locality into the overall prediction.

What I mean is, why GP instead of, say polynomial regression, or training a NN on the CNN outputs?

I'm curious if the decision was due to some feature of GPs that you really needed, and could not easily get otherwise. I'm not doubting your choice, really, I'm just trying to get a feeling for what niche GPs fill (given their difficult computational requirements), as I have never found the opportunity to use them.

> that would help to incorporate locality into the overall prediction.

Not sure what you mean by this so I am sure I am missing something.

We used a Gaussian process on the outputs from a CNN, not separately.

What I meant by locality is the idea that two points which are close to each other should also have predictions close to each other, assuming that there is some underlying Gaussian variable we can't see immediately from the data or from the CNN.

One example could be housing prices: if your neighbor sold their house, controlled for other factors like square ft or the number of bedrooms, you should value your own house similarly. We assume that there is this underlying quality of the neighborhood, possibly given by distance between houses.

The niche that GPs fill is broad. Two applications I've seen were their use in trying to find the missing MH370 flight, and another is in hyperparameter tuning of ML models.

Reply on news.ycombinator.com