1
$\begingroup$

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)

enter image description here

(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.

$\endgroup$

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.