Suppose that both $A$ and $B$ are $p \times p$ positive definite matrices.
I want to compute $X=\left[ (A\otimes I_p) + (I_p \otimes A) \right] \left[ (B \otimes I_p) + (I_p \otimes B) \right]^{-1} (B^2 \otimes B^2)$.
When $p$ is small, it is straightforward to obtain $X$. However, when $p$ is large, since $X$ includes a huge matrix of dimension $p^2 \times p^2$, computation of $X$ is very time consuming or almost impossible.
I'm wondering if there are efficient ways to compute $X$ even for large $p$.
By the spectral theorem we can write $$ B = U D U^* $$ for some $D>0$ and unitary $U$. Then $$ B \otimes I + I \otimes B = (U \otimes U) (D \otimes I + I \otimes D) (U^* \otimes U^*). $$ Thus To compute the inverse of $B \otimes I + I \otimes B$ we need only compute $$ (U \otimes U) (D \otimes I + I \otimes D)^{-1} (U^* \otimes U^*) $$ and now the inverse is on a diagonal matrix (which we compute from $B$ only) and hence is much faster. We can also absorb the multiplication by $B^2 \otimes B^2$ into this computation as $$ (B \otimes I + I \otimes B)^{-1}(B^2 \otimes B^2) = (U \otimes U) (D \otimes I + I \otimes D)^{-1}(D^2 \otimes D^2) (U^* \otimes U^*), $$ which then implies the final multiplication can be reduced to multplication by a diagonal matrix.