How to calculate the inverse of sum of a Kronecker product and a diagonal matrix I want to calculate the inverse of a matrix of the form  $S = (A\otimes B+C)$, where $A$ and $B$ are symetric and invertible, $C$ is a diagonal matrix with positive elements. Basically if the dimension is high, direct calculation can be expensive. So I wonder if it is possible to speed up the calculation taking advantage of some algebraic structures? 
For example, let $A$ be a $m\times m$ matrix and $B$ be a $n\times n$ matrix, then we have $(A\otimes B)^{-1}=A^{-1}\otimes B^{-1}$ so that the approximate cost in floating point operations (flops) can be reduced from $7m^3n^3/3$ to $7(m^3+n^3)/3+m^2n^2$.
Since the matrix $S$ also has some special structure, any idea on the efficient calculation of $S^{-1}$? 
Many thanks!


*

*In addition, what about $S^{-1}a$, where $a$ is a $mn\times 1$ vector.

 A: Well, we know that:
$$(\mathbf{A} \otimes \mathbf{B})^{-1} = \mathbf{A}^{-1} \otimes \mathbf{B}^{-1}$$ 
(and $(\mathbf{A} \otimes \mathbf{B})$ is invertible $\iff$ $\mathbf{A}$ 
 and $\mathbf{B}$ are invertible as well). So that takes care of the first term. 
$\mathbf{C}$ is invertible with inverse equal to the diagonal matrix with diagonal elements formed of the element-wise inverses of the diagonal entries of $\mathbf{C}$. 
We also know that 
if $\mathbf{C}$ is diagonal:
$$\left(\mathbf{D}+\mathbf{C}\right)^{-1} = \mathbf{D}^{-1} - \mathbf{D}^{-1} \left(\mathbf{C}^{-1}+\mathbf{D}^{-1} \right)^{-1}\mathbf{D}^{-1},$$
where $\mathbf{D}=\mathbf{A} \otimes \mathbf{B}$. Next, you use the fact that:
$$\left(\mathbf{C}^{-1}+\mathbf{D}^{-1} \right)^{-1}=\mathbf{C}\left(\mathbf{C}+\mathbf{D}\right)^{-1}\mathbf{D}.$$
Plugging this in the previous formula yields:
$$\left(\mathbf{D}+\mathbf{C}\right)^{-1} = \mathbf{D}^{-1} - \mathbf{D}^{-1} \mathbf{C}\left(\mathbf{C}+\mathbf{D}\right)^{-1},$$
moving things around:
$$\left(\mathbf{D}+\mathbf{C}\right)^{-1} = \left(\mathbf{I}+\mathbf{D}^{-1}\mathbf{C}\right)^{-1}\mathbf{D}^{-1},$$
where $\mathbf{I}$ is the identity matrix of same size as $\mathbf{D}$.
Now, given the spectral decomposition of $\mathbf{D}^{-1}\mathbf{C}$, that of $\left(\mathbf{I}+\mathbf{D}^{-1}\mathbf{C}\right)^{-1}$ is easy to compute because the spectral decomposition of $\mathbf{I}+\mathbf{D}^{-1}\mathbf{C}$ is related to that of $\mathbf{D}^{-1}\mathbf{C}$ (the eigenvector of the sum of the two terms are the same as the eigenvector of the second term and the eigen-value of the sum equals the sum of the eigenvalues of the two terms).
Because of the fact that $\mathbf{C}$ is diagonal, the spectral decomposition of $\mathbf{D}^{-1}\mathbf{C}$ itself is easy to obtain from that of $\mathbf{D}^{-1}$.
Now, you still have to compute the spectral decomposition of $\mathbf{D}$ but I think that the fact that $\mathbf{D}$ is the product of two Kronecker matrices helps considerably here.
