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?