在SymPy中加速符号行列式的计算



我有一个4x4矩阵A,它的每个条目中都有相当长但简单的符号表达式。大约涉及30种不同的符号。我所说的"简单"是指这些符号只使用加法/减法、乘法/除法和整数幂进行组合。我所说的"长"是指,如果我打印出矩阵,它覆盖了三到四个屏幕。

我需要这个矩阵的行列式。或者,更具体地说,我知道行列式是一个特定符号中的四阶多项式,我需要这个多项式的系数。A.det()在运行数小时后不会终止,所以我需要一种不同的方法。有什么想法吗?到目前为止,我已经尝试在A的每个元素中抛出各种simplify函数,但没有成功。

有没有什么策略可以让SymPy知道我的表达式的简单结构,或者我知道结果是其中一个符号中的多项式?

也许可以创建4x4行列式的通用表达式

In [30]: A = Matrix(4, 4, symbols('A:4:4'))
In [31]: A
Out[31]:
⎡A₀₀  A₀₁  A₀₂  A₀₃⎤
⎢                  ⎥
⎢A₁₀  A₁₁  A₁₂  A₁₃⎥
⎢                  ⎥
⎢A₂₀  A₂₁  A₂₂  A₂₃⎥
⎢                  ⎥
⎣A₃₀  A₃₁  A₃₂  A₃₃⎦
In [32]: A.det()
Out[32]:
A₀₀⋅A₁₁⋅A₂₂⋅A₃₃ - A₀₀⋅A₁₁⋅A₂₃⋅A₃₂ - A₀₀⋅A₁₂⋅A₂₁⋅A₃₃ + A₀₀⋅A₁₂⋅A₂₃⋅A₃₁ + A₀₀⋅A₁₃⋅A₂₁⋅A₃₂ - A₀₀⋅A₁₃⋅A₂₂⋅A₃₁ - A₀₁⋅A₁₀⋅A₂₂⋅A₃₃ + A₀₁⋅A₁₀⋅A₂₃⋅A₃₂ + A₀₁⋅A₁₂⋅A₂₀⋅
A₃₃ - A₀₁⋅A₁₂⋅A₂₃⋅A₃₀ - A₀₁⋅A₁₃⋅A₂₀⋅A₃₂ + A₀₁⋅A₁₃⋅A₂₂⋅A₃₀ + A₀₂⋅A₁₀⋅A₂₁⋅A₃₃ - A₀₂⋅A₁₀⋅A₂₃⋅A₃₁ - A₀₂⋅A₁₁⋅A₂₀⋅A₃₃ + A₀₂⋅A₁₁⋅A₂₃⋅A₃₀ + A₀₂⋅A₁₃⋅A₂₀⋅A₃₁ - A₀₂⋅A₁
₃⋅A₂₁⋅A₃₀ - A₀₃⋅A₁₀⋅A₂₁⋅A₃₂ + A₀₃⋅A₁₀⋅A₂₂⋅A₃₁ + A₀₃⋅A₁₁⋅A₂₀⋅A₃₂ - A₀₃⋅A₁₁⋅A₂₂⋅A₃₀ - A₀₃⋅A₁₂⋅A₂₀⋅A₃₁ + A₀₃⋅A₁₂⋅A₂₁⋅A₃₀

然后用类似的东西替换条目

A.det().subs(zip(list(A), list(your_matrix)))

不过,SymPy生成4x4行列式的速度较慢是一个缺陷。你应该在https://github.com/sympy/sympy/issues/new.

编辑(这不适合评论)

看起来Matrix.det正在调用一个简化函数。对于3x3及以下的矩阵,行列式公式是明确写出的,但对于较大的矩阵,它是使用Bareis算法计算的。您可以看到简化函数(cancel)在这里的调用位置,这是计算的一部分,但最终会做很多工作,因为它试图简化非常大的表达式。可能更明智的做法是只做必要的简化,以抵消行列式本身的项。我为此开了一期。

另一种加快速度的可能性是选择一种不同的行列式算法,我不确定它是否有效。选项为Matrix.det(method=alg),其中alg"bareis"(默认值)、"berkowitz""det_LU"之一。

您所描述的多项式比率就是所谓的有理函数:https://en.wikipedia.org/wiki/Rational_function

SymPy的polys模块确实有表示有理函数的方法,尽管它们可能很慢,尤其是在有很多变量的情况下。

sympy 1.7中有一个新的矩阵实现,它仍然是实验性的,但基于polys模块,可以处理有理函数。我们可以通过快速创建一个随机矩阵来测试它:

In [35]: import random                                                                                                            
In [36]: from sympy import random_poly, symbols, Matrix                                                                           
In [37]: randpoly = lambda : random_poly(random.choice(symbols('x:z')), 2, 0, 2)                                                  
In [38]: randfunc = lambda : randpoly() / randpoly()                                                                              
In [39]: M = Matrix([randfunc() for _ in range(16)]).reshape(4, 4)                                                                
In [40]: M                                                                                                                        
Out[40]: 
⎡     2              2            2              2       ⎤
⎢  2⋅z  + 1       2⋅z  + z     2⋅z  + z + 2     x  + 2   ⎥
⎢  ────────     ────────────   ────────────   ────────── ⎥
⎢   2            2                 2             2       ⎥
⎢  y  + 2⋅y     y  + 2⋅y + 1      x  + 1      2⋅z  + 2⋅z ⎥
⎢                                                        ⎥
⎢  2              2                  2        2          ⎥
⎢ y  + y + 1   2⋅x  + 2⋅x + 1       z        z  + 2⋅z + 1⎥
⎢ ──────────   ──────────────     ──────     ────────────⎥
⎢     2            2               2           2         ⎥
⎢  2⋅y  + 2       y  + 2⋅y        y  + 1      x  + x + 2 ⎥
⎢                                                        ⎥
⎢     2           2                2            2        ⎥
⎢  2⋅z  + 2    2⋅z  + 2⋅z + 2     y  + 1     2⋅y  + y + 2⎥
⎢────────────  ──────────────   ──────────   ────────────⎥
⎢   2             2                2             2       ⎥
⎢2⋅z  + z + 1  2⋅x  + 2⋅x + 2   2⋅y  + 2⋅y      x  + 2   ⎥
⎢                                                        ⎥
⎢    2               2            2             2        ⎥
⎢ 2⋅y  + 2⋅y      2⋅y  + y     2⋅x  + x + 1  2⋅x  + x + 1⎥
⎢ ──────────      ────────     ────────────  ────────────⎥
⎢    2              2                 2          2       ⎥
⎣   z  + 2         x  + 2          2⋅y          x  + 2   ⎦

如果我们将其转换为新的矩阵实现,那么我们可以使用charpoly方法计算行列式:

In [41]: from sympy.polys.domainmatrix import DomainMatrix                                                                        
In [42]: dM = DomainMatrix.from_list_sympy(*M.shape, M.tolist())                                                                  
In [43]: dM.domain                                                                                                                
Out[43]: ZZ(x,y,z)
In [44]: dM.domain.field                                                                                                          
Out[44]: Rational function field in x, y, z over ZZ with lex order
In [45]: %time det = dM.charpoly()[-1] * (-1)**M.shape[0]                                                                         
CPU times: user 22 s, sys: 231 ms, total: 22.3 s
Wall time: 23 s

这比@asmeurer上面建议的方法慢,但它以正则形式产生作为展开多项式比率的输出。特别是,这意味着你可以立即判断行列式是否为零(对于所有x,y,z)。该时间也由等效的cancel占用,但实现比Matrix.det更有效。

这需要多长时间在很大程度上取决于最终输出的复杂程度,你可以从它的字符串表示的长度中了解到这一点(我不会展示全部内容!):

In [46]: len(str(det))                                                                                                            
Out[46]: 54458
In [47]: str(det)[:80]                                                                                                            
Out[47]: '(16*x**16*y**7*z**4 + 48*x**16*y**7*z**2 + 32*x**16*y**7 + 80*x**16*y**6*z**4 + '

在某种程度上,应该可以将其集成到主Matrix类中,或者以其他方式使DomainMatrix类更易于公开访问。

最新更新