Decomposing a locally stationary covariance matrix

Question

Say I have a non-stationary Gaussian Process with a square exponential covariance whose shape varies throughout space. The covariance entries are:

$$ K_{ij} = N(|x_i-x_j|,\sigma_i^2+\sigma_j^2) $$

Where $\sigma_i^2$ is the local variance at point $i$. (This covariance is constructed from a non-uniform convolution--for more info, see Paciorek & Schervish, Environmetrics 2006).

Evaluating the Gaussian Process likelihood requires inverting (or at least computing the Cholesky/Eigen decomposition of) K:

$$ p(D)\sim N(0,K+\sigma_n^2I) $$

Where $\sigma_n$ is the independent noise term. That is potentially a very difficult task when K is large, and I want to be able to do this for ~$10^6$ data points. Are there any tricks to speeding this decomposition? Here are two simplifying assumptions:

• the $\sigma_i$ terms change slowly, so the covariance is locally stationary (not sure how to rigorously define that - Paciorek has the $\sigma$ terms come from another GP).
• the data are on a regular 3D grid (if it were stationary I could use super fast Kronecker techniques, see Gilboa 2015).

Current approach: "brute force" the decomposition with sparse inversion techniques, but this really only works in a very small number of cases.

Answer 1

There are a few options:

• A very fast method that kind of breaks GPs but seems to work well in many practical cases are local expert models such as robust Bayesian Committee Machines (paper here, very easy to implement!)

• A more mathematically beautiful approach is to use hierarchical matrices eg HODLR matrix approximations (paper here, fast code available here)

• Low rank approximations eg FITC, Nystrom, etc. Work well when the model requires only a few degrees of freedom (many implementations out there like GPy so just google)

• Stochastic Variational Inference scales based on the number of iterations applied in an iterative method (paper here, code here )

I just named a few approaches above but there are also spectral methods, KISSgps, etc. In fact it is one of the most heavily researched areas in Gaussian Process literature!

Good luck and let us know which works best for you :)