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