How to efficiently invert $K \otimes M+I_T\otimes \Sigma$?

Question

I'm looking for a way to efficiently invert $$K \otimes M+I_T\otimes \Sigma$$ where the inverses for $M,K$ exist. $I_T$ is the identity matrix of dimension $T$, and $\Sigma$ is a diagonal matrix, with positive only elements along the main diagonal.

If I had only $K \otimes M$, I would just do $K^{-1}\otimes M^{-1}$. However, with the extra term, I don't see how to efficiently invert it.

Answer 1

Generally there isn't a way to compute the inverse of a sum of Kronecker products. However, suppose there is a factor in common, let's say $I_T$ here and your sum is

$$ A = K \otimes I_T + I_T \otimes \Sigma $$

Then solving linear systems with $A$ becomes equivalent to solving a Sylvester equation. Thus perhaps one way to approach your problem would be to try and introduce $I_T$ as a common factor by say

$$ (I \otimes M)^{-1} (K \otimes M + I_T \otimes \Sigma) = K \otimes I_T + I_T\otimes \Sigma M^{-1} $$

You can find a list of helpful formulas for Kronecker products on this wikipedia page

Answer 2

Expanding on my previous comment. There are efficient algorithms to solve linear systems of the form $$ (K \otimes M+I_T\otimes \Sigma)\operatorname{vec}(X) = \operatorname{vec}(B), $$ so in practice you can just use those instead of an explicit inverse. Explicit inverses are overrated, in some fields. :)

The trick is that these linear systems are equivalent to the matrix equation $$ MXK^\top + \Sigma X I_T^\top = B, $$ which is known as generalized Sylvester equation.

There are around efficient algorithms to solve these equations in $O(\text{larger_dimension_of_B}^3)$ For instance, see https://dl.acm.org/citation.cfm?id=146929 or https://people.cs.umu.se/isak/recsy/ . Basically the idea is reducing all coefficients to triangular form simultaneously using two QZ decompositions, and then the resulting system can be solved by back-substitution.

In your case you can get a slight simplification, by exploiting the fact that one of the coefficients is an identity. Set $\Sigma X = Y$, to get $$ (M\Sigma^{-1})YK^\top + Y = B, $$ which is another classical matrix equation (a discrete-time Sylvester equation), and can be solved with a slightly faster algorithm, which is also implemented in a Matlab one-liner as Y = dlyap(-M/Sigma, K', -B). Unfortunately, a solver for this variant is not available in Numpy/Scipy, up to my knowledge; it has only a solver for the classical Sylvester equation, which you can use with the approach described in @Reid.Atcheson's answer.

In the case in which $M\Sigma^{-1}$ and $K^\top$ have all eigenvalues inside the unit circle, you can also write a closed-form solution for this equation as the infinite series $$ Y = \sum_{k=0}^\infty (-M\Sigma^{-1})^kB(K^\top)^k. $$ This follows from applying the Neumann series $(I-M)^{-1}=\sum_{k=0}^\infty M^k$ to the Kroneckerized matrix $(I_{T^2} - K\otimes (-M\Sigma^{-1}))$.