How to calculate the inverse of the sum of an identity and a Kronecker product efficiently?

Question

I have a matrix $K$ which is the sum of a identity and a Kronecker product of two symmetric matrices as following and I want to calculate the inverse of it $K^{-1}$. \begin{eqnarray} K=\mathbf{I}_{mn}+XX^\top\otimes YY^\top \end{eqnarray} where $X$ and $Y$ are $m\times p$ and $n\times q$ matrices respectively and $m>p,n>q$.

Since the size of Kronecker product of two matrices is the product of their sizes, the matrix to be inverted is very large. Can this calculation be simplified?

Using Woodbury matrix identity, the matrix to be inverted can be reduced from $mn\times mn$ to $pq\times pq$ as \begin{eqnarray} \left(\mathbf{I}+XX^\top\otimes YY^\top\right)^{-1} &=&\left(\mathbf{I}+(X\otimes Y)(X\otimes Y)^\top\right)^{-1}\\\ &=&\mathbf{I}-(X\otimes Y)\left(\mathbf{I}+(X\otimes Y)^\top(X\otimes Y)\right)^{-1}(X\otimes Y)^\top \end{eqnarray} However, in my application, $p$ and $q$ are not much less than $m$ and $n$. Thus inverting a $pq\times pq$ matrix is still a consuming work and I want to find a more simple calculation of $K^{-1}$.

Consider the following problem. If there is not the identity matrix $\mathbf{I}_{m\times n}$, we can invert $XX^\top$ and $YY^\top$ respectively and then calculate the Kronecker product as $$ (XX^\top\otimes YY^\top)^{-1}=(XX^\top)^{-1}\otimes(YY^\top)^{-1} $$ using the property of Kronecker product. However, if the identity item exists, this property cannot be utilized.

Then my question is that if there is any way to decompose the calculation of $K^{-1}$ into inverses of some smaller matrices whose size is linear with $m$ and $n$ but not their products? Or is it possible to extract the Kronecker product $\otimes$ out of the inverse?

If you have any suggestion or idea, please let me know. Thank you very much for your help!

Answer 1

I assume spectral decomposition of a symmetric matrix does not cost too much. Let $PXX^TP=\Lambda_1$, $QYY^TQ=\Lambda_2$ be the spectral decomposition of $XX^T, YY^T$, respectively. Here $\Lambda_1, \Lambda_2$ are diagonal and $P, Q$ are orthogonal. Then $$K^{-1}=(P\otimes Q)^T(I+\Lambda_1\otimes \Lambda_2)^{-1}(P\otimes Q).$$

$(I+\Lambda_1\otimes \Lambda_2)^{-1}$ can be read directly.

Answer 2

Excellent answer by Betrand. But it's easy to verify that the final form should be $K^{-1} = (P\otimes Q)(I + \Lambda_1 \otimes \Lambda_2)^{-1} (P \otimes Q)^T$.

Here is some R code:

A = matrix(c(3, 2, 2, 3), ncol = 2)
B = matrix(c(2, 2, 3, 4), ncol = 2)
AA = A %*% t(A)
BB = B %*% t(B)
AA_decomp = eigen(AA)
P = AA_decomp$vectors
LAMB_A = diag(AA_decomp$values)
P %*% LAMB_A %*% t(P)
BB_decomp = eigen(BB)
Q = BB_decomp$vectors
LAMB_B = diag(BB_decomp$values)
Q %*% LAMB_B %*% t(Q)
K = diag(dim(A)[1] * dim(B)[1]) + kronecker(AA, BB)
solve(K)
kronecker(P, Q) %*% solve(diag(dim(A)[1] * dim(B)[1]) + kronecker(LAMB_A, LAMB_B)) %*% t(kronecker(P, Q))