Inverse of a diagonal matrix plus a Kronecker product? Given two matrices $X$ and $Y$, it's easy to take the inverse of their Kronecker product:
$(X\otimes Y)^{-1} = X^{-1}\otimes Y^{-1}$
Now, suppose we have some diagonal matrix $\Lambda$ (or more generally an easily inverted matrix, or one for which we already know the inverse).  Is there a closed-form expression or efficient algorithm for computing $(\Lambda + (X\otimes Y))^{-1}$?
 A: yes there is.
See equation 5 in http://books.nips.cc/papers/files/nips24/NIPS2011_0443.pdf
Stegle et al. Efficient inference in matrix-variate Gaussian models
with iid observation noise
A: Let $C := \left\{ c_{i,j} \right\}_{i,j=1}^N$ and $A := \left\{ a_{i,j} \right\}_{i,j=1}^n$ be symmetric matrices. The spectral decompositions of the two matrices read $A = O^T D_A O$ and $C = U^T  D_C U$ where $D_A := Diag(\lambda_k)_{k=1}^n$ and $D_C := Diag(\mu_k)_{k=1}^N$ and $O \cdot O^T = O^T \cdot O = 1_n$ and $U \cdot U^T = U^T \cdot U = 1_N$. We use equation 5 from the cited paper in order to compute the resolvent $G_{A \otimes C}(z) := \left(z 1_{ n N} - A \otimes C\right)^{-1}$. We have:
\begin{eqnarray}
G_{A \otimes C}(z) &=& \left( O^T \otimes U^T\right) \cdot \left( z 1_{ n N} - D_A \otimes D_C \right)^{-1} \cdot \left(O \otimes U \right) \\
&=& \left\{ \sum\limits_{k=1}^n O_{k,i} O_{k,j}  U^T \cdot \left( z 1_{N} - \lambda_k D_C\right)^{-1} U \right\}_{i,j=1}^n \\
&=& \sum\limits_{p=0}^\infty \frac{1}{z^{1+p}} A^p \otimes C^p \\
&=& \frac{\sum\limits_{p=0}^{d-1} z^{d-1-p} \sum\limits_{l=d-p}^d (-1)^{d-l} {\mathfrak a}_{d-l} \left(A \otimes C\right)^{p-d+l}}{\sum\limits_{p=0}^d z^{d-p} (-1)^p {\mathfrak a}_p}
\end{eqnarray}
where 
$\det( z 1_{n N} - A \otimes C) := \sum\limits_{l=0}^{n N} (-1)^l {\mathfrak a}_l z^{n N-l}$.
The first two lines from the top are straightforward. In the third line we expanded the inverse matrix in a series and finally in the fourth line we summed up the infinite series using Cayley-Hamilton's theorem.
A: Use the following trick:
$$(\Lambda+A)^{-1}=A^{-1}(A^{-1}\Lambda+I)^{-1}$$
where $A=X\bigotimes Y$.
Then we  can decompose $\Lambda=\Lambda_1\bigotimes \Lambda_2$ by any feasible $\Lambda_1$ and $\Lambda_2$ because $\Lambda$ is diagonal. Now we have
$$A^{-1}\Lambda=X^{-1}\Lambda_1 \bigotimes Y^{-1}\Lambda_2$$
and use the method mentioned by David Rohde
