How to calculate the inverse of sum of a Kronecker product and a diagonal matrix

Question

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.

Answer 1

Here's a method that might work well. Let $C$ have diagonal elements $d_1,\dots,d_{mn}$. Write $$ C = d_1 e_1e_1^T + \cdots + d_{mn} e_{mn} e_{mn}^T $$ Now, we can find the inverse by applying the Sherman-Morrison formula (iteratively) $mn$ times. That is, we've broken the computation into several rank-1 updates.

What's nice is that this allows you to leverage your formula for $(A \otimes B)^{-1}$. I'm not sure how efficient that will be in practice, but it's a start.

Answer 2

Here is an answer that only works for computing $S^{-1}a$ and only works when S itself is positive definite (which is the case if $A$ and $B$ is positive definite).

The conjugate gradient algorithm computes $x=S^{-1}a$ by doing at most $mn$ multiplications of the form $S x_k$ where $x_0$ is an initial guess. You can get the exact algorithm from http://en.wikipedia.org/wiki/Conjugate_gradient_method. You have
$$ S x_k = S \mathrm{vec}(X_k) = (A\otimes B + C) \mathrm{vec}(X_k) = \mathrm{vec}(BX_kA^\top + K\odot X_k) $$ where $\odot$ is the elementwise matrix product, $x_k=\mathrm{vec}(X_k)$ and $C=\mathrm{diag}(\mathrm{vec}(K))$. Therefore $Sx_k$ can be computed in $m^2n+n^2m+nm$ operations. Since we at most need to do $mn$ of these vector matrix multiplications the total computation is $\mathcal{O}(m^3n^2+n^3m^2)$, which is at least a little better than $m^3n^3$. Actually it's a lot better when you consider how much memory the algorithms require -- $\mathcal{O}(m^2+n^2)$ instead of $m^2n^2$.