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

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.

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}$$

• There are also generalized Sylvester equations $AXB+CXD=E$ which correspond exactly to OP's problem. One can solve them with generalizations of the Bartels-Stewart algorithm, based on the QZ decomposition rather than the Schur decomposition. – Federico Poloni May 24 '19 at 8:13
• In OP's case, after some scaling to replace $\Sigma$ with an identity, one has a so-called discrete-time Sylvester equation $AXB+C=E$, which can be solved, for instance, with Matlab's dlyap. – Federico Poloni May 24 '19 at 8:15
• @FedericoPoloni, Thanks for your comments. Could you elaborate a bit more on them? What are A,B,X,C,D and E? Could you provide a link so that I can study what you purpose I should use? Thanks – An old man in the sea. May 24 '19 at 15:18
• Reid, thanks for your answer. ;) – An old man in the sea. May 24 '19 at 15:26
• I forgot to +1 Reid ;) – An old man in the sea. May 25 '19 at 7:27

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}))$$.