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

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? • Using the tensor decomposition seems to be a good idea. But what about your predictions? Given the sparsity and peculiarity of your design how can you be sure the extrapolation involved in predicting will make sense? Your problem seems to be highly structured why do you think a radially symmetric kernel makes sense? And what about your prior mean - is it zero? – g g Dec 11 '15 at 21:20
• Tensor decomposition is the same as the Kronecker one, right? For your questions: 1) I'm not sure it's so peculiar, I would imagine this is quite common. Though now I'm realizing that I need to add an independent noise component to the covariance, which maybe means it can't be separated. 2) The radially symmetric kernel is just the way the data is produced, it's basically a convolution of a white noise process plus deterministic signal. 3) For simplicity I've set the mean to 0 since you can always subtract it from the data. – cgreen Dec 11 '15 at 21:32
• Great idea with the tensor product! May I sugest at the 1D level then using Toeplitz matrix inversion. Together they will be very fast indeed! – j__ Dec 12 '15 at 10:27
• @j__thanks! And I'll try that! Problem is I don't think the decomposition applies for added independent noise: $(K + \sigma_n I)$. Any ideas? – cgreen Dec 12 '15 at 17:07

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.