Decomposing large symmetric banded sparse matrices

Question

I'm investigating 3D image deblurring and one of the approaches I'm interested in is applying spatial regularisation. To do this I have generated a matrix $A$ which encodes the 6-connectivity of each image voxel (i.e. the adjacent voxels in $x$,$y$, and $z$ planes). For the specific images I am working with this results in a large (5242880x5242880), sparse, symmetric, banded, matrix with $1$s on the first, 512-nd and 262144-th super- and sub- diagonals. All other entries are $0$.

To compute the regularised solution I need to obtain the inverse of $A$. It seems intuitive that it should be possible to factorise or otherwise easily decompose this highly structured binary matrix into smaller matrices to make this problem more tractable?

Thanks

Answer 1

One efficient method (using Hessenberg matrices) is described here: The inverse of banded matrices (2013).

Let $B_{r,n}$ ($1 \leq r \leq n$) be an $n \times n$ matrix of entries $\{a_{ij}\}$, $−r \leq i \leq r$, $1 \leq j \leq r$, with the remaining un-indexed entries all zeros. In this paper we give the LU factorization and the inverse of the matrix $B_{r,n}$, using a method based on Hessenberg submatrices associated to $B_{r,n}$.