Writing the SVDs of symmetric positive-definite matrices $A=U\_AD\_AU\_A^T$ and $B=U\_BD\_BU\_B^T$, the following inversion can conveniently be simplified:
\begin{align\*}
(I + A \otimes B)^{-1}
&= ((U\_A \otimes U\_B)(U\_A \otimes U\_B)^T + (U\_A \otimes U\_B)(D\_A \otimes D\_B)(U\_A \otimes U\_B)^T)^{-1} \\
&=(U\_A \otimes U\_B)(I + D\_A \otimes D\_B)^{-1}(U\_A \otimes U\_B)^T.
\end{align\*}
That is, computing the inverse of $I + A \otimes B$ can be reduced to computing the SVDs of $A$ and $B$ and a cheap inversion of a diagonal matrix.


It is possible to similarly simplify $I + A \otimes B + C \otimes D$, where $A$, $B$, $C$, and $D$ are symmetric positive-definite matrices and the Kronecker products are compatible? That is, how can one favourably exploit the structure in $I + A \otimes B + C \otimes D$ in computing its inverse?


Edit: To clarify, for my particular use case, I am interested in computing $(I + A \otimes B + C \otimes D)^{-1} E$ for some matrix $E$; I don't necessarily need the inverse itself.



Inversion of the sum of an identity matrix and two Kronecker products
Solving a linear system with that matrix is equivalent to solving the linear matrix equation $X + B^TXA + D^TXC = E$. As far as I know, it is an open problem how to exploit the structure of that matrix product to solve it (with a direct algorithm) in fewer than $O(n^5)$ operations (which is the complexity of $n^2$ steps of GMRES using the Kronecker structure to compute the products). It would be a big surprise to find a cheaper algorithm.


If one of the terms is small in norm or in rank, you can obtain a preconditioner for an iterative method by dropping it. That is currently your best bet to solve this problem. Pointers to some research in this direction: <https://arxiv.org/abs/1704.02167>, <https://doi.org/10.1016/j.laa.2017.06.027>.


