Efficient, stable inverse of a patterned covariance matrix for gridded data

Question

I have computed a stationary covariance matrix defined for data on a grid. The data y are regularly spaced in 3D, lexicographically ordered in the covariance matrix, and I'm using a using a square exponential kernel (EDIT: plus independent noise): $$ K_{ij} = \exp\left(-\frac{\lvert\mathbf{x_i}-\mathbf{x_j}\rvert^2}{2\sigma^2}\right) + \sigma_n\delta_{ij} $$ A picture of this covariance is below. As you can see it's highly banded. I would like to compute the likelihood of a Gaussian process with this covariance: $$ \begin{align} \mathcal{L(y;\sigma}) &= \mathcal{N}(\mathbf{y}|\mathbf{0},K) \\ &=(2\pi)^{-\frac{N}{2}}\lvert K\rvert^{-\frac{1}{2}}\exp\left(-\frac{1}{2}\mathbf{y}'K^{-1}\mathbf{y}\right) \end{align} $$ So in fact I don't need to store the inverse. For numerical stability I compute the Cholesky decomposition to get this likelihood, according to Rasmussen: $$ \begin{align} L&:=\text{cholesky}(K)\\ \alpha&:=L' \backslash (L\backslash \mathbf{y}) \\ \log \mathcal{L}(y)&:=-\frac{1}{2}\mathbf{y}'\mathbf{\alpha}-\sum_i L_{ii} - \frac{n}{2}\log(2\pi) \end{align} $$ Question: how can I take advantage of the sparseness, stationarity, and isotropy of this covariance to get a way more efficient likelihood calculation? Not just of the Cholesky decomposition (or whatever other decomposition you suggest) but also the operations that follow, leading to the final likelihood calculation.

What I've done so far

I've just learned about the Kronecker product, which has this helpful property: $$ (A\otimes B)^{-1} = A^{-1}\otimes B^{-1} $$ In my case the separability of the kernel into X,Y, and Z directions means (I think): $$ K = K^X \otimes K^Y \otimes K^Z $$ where: $$ K^X_{ij} = \exp\left(-\frac{(x_0(i-j))^2}{2\sigma^2}\right) $$ and similarly for Y and Z directions, where $x_0$ is the grid spacing in $x$, and $\lvert\cdot\rvert$ is the absolute value. Assuming my data is on a cube with side $D$, this reduces the problem from an inverse of a $D^3 \times D^3$ matrix to three $D \times D$ matrices. What else can I do to further simplify, while maintaining numerical stability?

EDIT can this decomposition still work with the independent noise along the diagonal? And if I change $\sigma_n$ during simulation, do I have to completely recalculate the cholesky decomposition?

Answer 1

I have found the best answer to be in "Scaling Multidimensional Inference for Structured Gaussian Processes" by Gilboa et al, IEEE Transactions on Pattern Analysis and Machine Intelligence, Feb 2015.

They show you can compute $y^TK^{-1}y$ very quickly without taking any inverses or decompositions, when the kernel for K is separable into a product of x, y, and z kernels (as it is when the kernel is square exponential) and the data is on a regular grid. They note that under these conditions, $K=K_x \otimes K_y \otimes K_z$, and thus the eigenvector/value matrices of K are also the Kronecker products of the eigenvector/value matrices along x, y, and z. So you can easily compute:

$$ (K+\sigma_n^2I)^{-1}y = Q(V+\sigma_n^2I)^{-1}Q^Ty $$

...where Q are the vectors and V are the values, if you have an algorithm for computing the product of a bunch of kronecker products of matrices and a vector: $(\otimes A_i)b$. And they suggest a very nice efficient algorithm for computing exactly that. Cool paper, highly recommended! They also show you can take advantage of these techniques even when you have missing data.

Answer 2

The paper http://arxiv.org/pdf/1403.6573.pdf deals with tensor product of factors with noise and provides computational complexity for the algorithm you are considering. The paper uses Tensor structure of the matrix to get better computational complexities, while it doesn't use Cholesky decomposition, but SVD decomposition.

The results from the paper seems to work if noise is not uniform, however, I am not sure.