Inverse of a diagonal matrix plus a Kronecker product?

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

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.