我有一个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类更易于公开访问。