How to invert a sparse covariance matrix with only nearest-neighbor entries

Question

I have a set of random variables that have associated 3D coordinates (actually this is describing a 3D "image", where each R.V. is a voxel). Basically we can model it as a Markov random field. Under the Markov assumption, voxels are independent of non-neighbors given their nearest neighbors. So I could lay out a covariance matrix for this (huge) set of variables, but it would be quite sparse:

$$ \Sigma_{ij} = \begin{cases} \alpha,& \text{if } d_1(x_i,x_j)\leq 2\\ 0, & \text{otherwise} \end{cases} $$ where $d_1$ is the Manhattan distance. Are there any closed form ways to invert a matrix like this?

Use case: I want to evaluate a likelihood for a model of these R.V.'s using a multivariate normal: $p(D|M) = (2\pi)^{-3/2}|\Sigma|^{-1/2}\exp[-\frac{1}{2}(D-M)^T\Sigma^{-1}(D-M)]$ where $D$ are the measured values of the R.V.'s mentioned above and $M$ is my model of them. So I'm open to other ideas on how to speed up this computation when $D$ has potentially millions of voxels. E.g. perhaps it's not necessary to invert it for this computation?

Answer 1

No there are not closed form methods taking anything like $O(N^3)$ operations, for general sparse structured matrices. Finding such an inversion method would be the breakthrough of the half-century for fields like lattice gauge theory.

The best techniques involve some sort of Krylov-space sequence (for example, bi-conjugate gradient) usually with preconditioning (at least red-black preconditioning in the problem you pose). This works with just the length $N^3$ vectors and does not take steps involving the $N^6$ total matrix entries. If the condition number of the matrix (the ratio of largest to smallest eigenvalues) is not too big, this sort of scheme converges in roughly $O(N)$ iterations.

However, for the super-structured matrix $\Sigma_{ij}$ you pose, there may conceivably be a closed-form answer in terms of the single possible non-zero value $\alpha$. You can start by dividing by $\alpha$ to get a matrix of ones along all links and zero's elsewhere. Even for that matrix, which can be inverted in two dimensions, I don't think there is a closed form inverse in three dimensions.