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.
