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

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

  • $\begingroup$ 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. $\endgroup$ May 24 '19 at 8:13
  • 1
    $\begingroup$ 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. $\endgroup$ May 24 '19 at 8:15
  • $\begingroup$ @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 $\endgroup$ May 24 '19 at 15:18
  • $\begingroup$ Reid, thanks for your answer. ;) $\endgroup$ May 24 '19 at 15:26
  • $\begingroup$ I forgot to +1 Reid ;) $\endgroup$ 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}))$.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.