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? Banded covariance

  • $\begingroup$ 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? $\endgroup$
    – g g
    Commented Dec 11, 2015 at 21:20
  • $\begingroup$ 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. $\endgroup$
    – cgreen
    Commented Dec 11, 2015 at 21:32
  • $\begingroup$ 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! $\endgroup$
    – j__
    Commented Dec 12, 2015 at 10:27
  • $\begingroup$ @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? $\endgroup$
    – cgreen
    Commented Dec 12, 2015 at 17:07

2 Answers 2


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.


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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.