在 Sympy 中求解包含符号大小的 MatrixSymbols 的矩阵方程？

Question

作为介绍，我想指出，如果一个矩阵 A 由 2x2 模式中的 4 个子矩阵组成，其中对角矩阵是正方形，那么如果我们将其逆矩阵表示为 X，则子矩阵 @ 987654323@，很容易用手显示。

我试图对 4x4 子矩阵的矩阵做同样的事情，但手工操作相当乏味。所以我认为 Sympy 会有所帮助。但我不知道是怎么回事（我一开始只是试图重现 2x2 结果）。

我试过了：

import sympy as s

def blockmatrix(name, sizes, names=None):
    if names is None:
        names = sizes
    ll = []
    for i, (s1, n1) in enumerate(zip(sizes, names)):
        l = []
        for j, (s2, n2) in enumerate(zip(sizes, names)):
            l.append(s.MatrixSymbol(name+str(n1)+str(n2), s1, s2))
        ll.append(l)
    return ll

def eyes(*sizes):
    ll = []
    for i, s1 in enumerate(sizes):
        l = []
        for j, s2 in enumerate(sizes):
            if i==j:
                l.append(s.Identity(s1))
                continue
            l.append(s.ZeroMatrix(s1, s2))
        ll.append(l)
    return ll

n1, n2 = s.symbols("n1, n2", integer=True, positive=True, nonzero=True)
M = s.Matrix(blockmatrix("m", (n1, n2)))
X = s.Matrix(blockmatrix("x", (n1, n2)))
I = s.Matrix(eyes(n1, n2))
s.solve(M*X[:, 1:]-I[:, 1:], X[:, 1:])

但它只是返回一个空列表而不是结果。

我也试过了：

• 使用M*X==I，但只返回False（布尔值，而不是表达式）
• 输入方程式列表
• 使用带有commutative=False 而不是MatrixSymbols 的“普通”符号——这会导致GeneratorsError: non-commutative generators: (x12, x22) 出现异常

但都没有运气。

您能否展示如何使用 Sympy 得出类似于我为 X22 提供的示例的结果？

使用 MatrixSymbols 解决的最相似的其他问题似乎已经通过使用内部符号数组或类似的方法来解决。但由于我正在处理符号大小的 MatrixSymbol，这对我来说不是一个选择。

Answer 1

这就是你所说的 2x2 矩阵的矩阵吗？

>>> a = [MatrixSymbol(i,2,2) for i in symbols('a1:5')]
>>> A = Matrix(2,2,a)
>>> X = A.inv()
>>> print(X[1,1])  # [1,1] instead of [2,2] because indexing starts at 0
a1*(a1*a3 - a3*a1)**(-1)

[您表示不正确并指出上述不正确——这似乎是一个应该解决的问题。]

我不确定为什么没有实现，但我们可以手动解决如下：

>>> n = 2
>>> v = symbols('b:%s'%n**2,commutative=False)
>>> A = Matrix(n,n,symbols('a:%s'%n**2,commutative=False))
>>> B = Matrix(n,n,v)
>>> eqs = list(A*B - eye(n))
>>> for i in range(n**2):
...   s = solve(eqs[i],v[i])[0]
...   eqs[i+1:] = [e.subs(v[i],s) for e in eqs[i+1:]]
...
>>> s  # solution for v[3] which is B22
(-a2*a0**(-1)*a1 + a3)**(-1)

您可以将 n 更改为 3 并查看一个适度复杂的表达式。将其更改为 4 并手动检查结果，为“乏味”一词提供新定义；-)

待求解方程的特殊结构也可以加快求解速度：感兴趣的变量是包含它的每一项中的最后一个因子：

>>> for i in range(n**2):
...   c,d = eqs[i].expand().as_independent(v[i])
...   assert all(j.args[-1]==v[i] for j in Add.make_args(d))
...   s = 1/d.subs(v[i], 1)*-c
...   eqs[i+1:] = [e.subs(v[i], s) for e in eqs[i+1:]]