通常没有办法计算克罗内克积之和的倒数。但是，假设有一个共同因素，假设这里的和你的总和是$I_T$

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

求解线性系统就等同于求解Sylvester 方程。因此，解决您的问题的一种方法可能是尝试将作为一个共同因素引入$A$$I_T$

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

您可以在此维基百科页面上找到 Kronecker 产品的有用公式列表

扩展我之前的评论。有有效的算法来解决以下形式的线性系统

$$(K \otimes M+I_T\otimes \Sigma)\operatorname{vec}(X) = \operatorname{vec}(B),$$ 所以在实践中你可以只使用那些而不是显式的逆。在某些领域，显式逆被高估了。:)

诀窍是这些线性系统等价于矩阵方程

$$MXK^\top + \Sigma X I_T^\top = B,$$ 这被称为广义西尔维斯特方程。

有一些有效的算法可以解决这些方程$O(\text{larger_dimension_of_B}^3)$ 例如，请参阅https://dl.acm.org/citation.cfm?id=146929或https://people.cs.umu.se/isak/recsy/。基本上，这个想法是使用两个 QZ 分解同时将所有系数减少为三角形形式，然后可以通过反向替换来求解得到的系统。

在您的情况下，您可以通过利用其中一个系数是一个单位这一事实来稍微简化一下。放$\Sigma X = Y$， 要得到

$$(M\Sigma^{-1})YK^\top + Y = B,$$ 这是另一个经典矩阵方程（离散时间 Sylvester 方程），可以用稍微快一点的算法求解，该算法也在 Matlab 单线中实现为Y = dlyap(-M/Sigma, K', -B)。不幸的是，据我所知，在 Numpy/Scipy 中没有此变体的求解器；它只有一个经典西尔维斯特方程的求解器，您可以将其与@Reid.Atcheson 的答案中描述的方法一起使用。

在这种情况下$M\Sigma^{-1}$和$K^\top$具有单位圆内的所有特征值，您也可以将此方程的封闭形式的解写为无限级数

$$Y = \sum_{k=0}^\infty (-M\Sigma^{-1})^kB(K^\top)^k.$$ 这是从应用诺依曼级数得出的$(I-M)^{-1}=\sum_{k=0}^\infty M^k$到 Kroneckerized 矩阵$(I_{T^2} - K\otimes (-M\Sigma^{-1}))$.