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.
