Decomposition of inverse of sum of block diagonal matrices

I have $$K$$ binary matrices $$Z_k$$ of dimension $$n \times q_k$$. These matrices are the kind you would get when performing one-hot encoding on a categorical variable, e.g. for $$n = 7$$ and a categorical variable vector $$[a,a,b,b,b,c,c]$$ (so $$q_1 = 3$$), in R:

Z1 <- model.matrix(~0 + rep(letters[1:3], times = c(2,3,2)))
Z1


$$\begin{bmatrix} 1 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 1 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 1 \end{bmatrix}$$

Obviously $$Z_kZ_k'$$ is a $$n \times n$$ block diagonal matrix:

Z1 %*% t(Z1)


$$\begin{bmatrix} 1 & 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 & 0 & 1 & 1 \end{bmatrix}$$

Further assume $$I_n$$ is the $$Z_0$$ $$n \times n$$ matrix.

I wish to calculate the inverse of $$V$$, which could be seen as a covariance matrix, and is a sum of block diagonal matrices:

$$V = \sum_{k = 0}^{K} \sigma^2_k \cdot Z_k Z_k',$$ where $$\sigma^2_k$$ are $$K + 1$$ positive variance constants.

Now if $$K = 1$$ (a single categorical variable) $$V$$ is in itself block-diagonal, and its inverse can be easily decomposed into $$q_1$$ inverses:

$$V^{-1} = \text{diag}(V_1^{-1}, \dots, V_{q_1}^{-1}),$$ where each $$V_j = \sigma^2_0 I_{n_j} + \sigma^2_j J_{n_j}$$. Here $$j = 1, \dots, q_1$$, and $$n_j$$ is the no. of rows (observations) in level $$j$$, so $$V_j$$ is a block of dimensions $$n_j \times n_j$$.

E.g. in R:

Z0 <- model.matrix(~0 + letters[1:7])
sig0 <- 1
ZZ0 <- sig0 * Z0 %*% t(Z0)

Z1 <- model.matrix(~0 + rep(letters[1:3], times = c(2,3,2)))
sig1 <- 1
ZZ1 <- sig1 * Z1 %*% t(Z1)

solve(ZZ0 + ZZ1)

          1          2     3     4     5          6          7
1  0.6666667 -0.3333333  0.00  0.00  0.00  0.0000000  0.0000000
2 -0.3333333  0.6666667  0.00  0.00  0.00  0.0000000  0.0000000
3  0.0000000  0.0000000  0.75 -0.25 -0.25  0.0000000  0.0000000
4  0.0000000  0.0000000 -0.25  0.75 -0.25  0.0000000  0.0000000
5  0.0000000  0.0000000 -0.25 -0.25  0.75  0.0000000  0.0000000
6  0.0000000  0.0000000  0.00  0.00  0.00  0.6666667 -0.3333333
7  0.0000000  0.0000000  0.00  0.00  0.00 -0.3333333  0.6666667



Can be computed with $$q_1 = 3$$ inverses:

jmat <- function(n_j) matrix(1, nrow = n_j, ncol = n_j)

vmat <- function(n_j) sig0 * diag(n_j) + sig1 * jmat(n_j)

Matrix::bdiag(solve(vmat(2)), solve(vmat(3)), solve(vmat(2)))

7 x 7 sparse Matrix of class "dgCMatrix"

[1,]  0.6666667 -0.3333333  .     .     .     .          .
[2,] -0.3333333  0.6666667  .     .     .     .          .
[3,]  .          .          0.75 -0.25 -0.25  .          .
[4,]  .          .         -0.25  0.75 -0.25  .          .
[5,]  .          .         -0.25 -0.25  0.75  .          .
[6,]  .          .          .     .     .     0.6666667 -0.3333333
[7,]  .          .          .     .     .    -0.3333333  0.6666667


When $$K > 1$$ this is not necessarily true anymore, because $$V$$ may not be block-diagonal anymore, e.g.:

Z2 <- model.matrix(~0 + rep(letters[1:2], times = c(3, 4)))
sig2 <- 2
ZZ2 <- sig2 * Z2 %*% t(Z2)

ZZ0 + ZZ1 + ZZ2

 1 2 3 4 5 6 7
1 4 3 2 0 0 0 0
2 3 4 2 0 0 0 0
3 2 2 4 1 1 0 0
4 0 0 1 4 3 2 2
5 0 0 1 3 4 2 2
6 0 0 0 2 2 4 3
7 0 0 0 2 2 3 4


We can see the inverse of $$V$$ is dense:

solve(ZZ0 + ZZ1 + ZZ2)

           1           2           3           4           5           6           7
1  0.60498221 -0.39501779 -0.11743772  0.02491103  0.02491103 -0.01423488 -0.01423488
2 -0.39501779  0.60498221 -0.11743772  0.02491103  0.02491103 -0.01423488 -0.01423488
3 -0.11743772 -0.11743772  0.41103203 -0.08718861 -0.08718861  0.04982206  0.04982206
4  0.02491103  0.02491103 -0.08718861  0.62455516 -0.37544484 -0.07117438 -0.07117438
5  0.02491103  0.02491103 -0.08718861 -0.37544484  0.62455516 -0.07117438 -0.07117438
6 -0.01423488 -0.01423488  0.04982206 -0.07117438 -0.07117438  0.61209964 -0.38790036
7 -0.01423488 -0.01423488  0.04982206 -0.07117438 -0.07117438 -0.38790036  0.61209964


However, I have hopes the inverse could also be decomposed somehow, or that someone would have anything useful or at least interesting to say about it! Because:

(a) $$V$$ is not block-diagonal but it is certainly "close" in some sense, as can be seen by its and its inverse actual values:

matplot <- function(A, ...) lattice::levelplot(A[, nrow(A):1], ...)

V <- matplot(ZZ0 + ZZ1 + ZZ2, main = "V")
V_inv <- matplot(solve(ZZ0 + ZZ1 + ZZ2), main = "V_inv")
gridExtra::grid.arrange(V, V_inv, ncol = 2)


(b) It is banded (here by 2, or in general by $$\lfloor \max(n_j)/2 \rfloor$$)

(c) It has $$\sum_{k = 0}^{K}\sigma^2_k$$, a constant value on its diagonal

(d) It is still the sum of such special matrices, and every element has a meaning (e.g. for $$\sigma^2_k = 1$$ for all $$k$$, the diagonal element is $$K + 1$$ and each $$[i,j]$$ element is the number of $$Z_kZ'_k$$ which have $$1$$ in their $$[i,j]$$ element.