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.

+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.

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

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

yes

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.