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!