加速 SymPy 中符号行列式的计算
Speeding up computation of symbolic determinant in SymPy
我有一个 4x4 矩阵 A，其每个条目中都有相当长但简单的符号表达式。涉及大约 30 个不同的符号。 "simple" 我的意思是这些符号仅使用 addition/subtraction、multiplication/division 和 整数 幂组合。 "long" 我的意思是，如果我打印出矩阵，它会覆盖三到四个屏幕。
我需要这个矩阵的行列式。或者，更具体地说，我知道行列式是一个特定符号的四阶多项式，我需要这个多项式的系数。 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) 的位置 here, which is necesssary as part of the computation, but end up doing a lot of work because it tries to simplify your very large expressions. It would probably be smarter to only do the simplifications that are needed to cancel terms of the determinant itself. I opened an issue。
另一种加快速度的可能性是 select 一种不同的行列式算法，我不确定是否可行。选项是 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 + '
在某些时候，应该可以将其集成到主矩阵中 class 或者以其他方式使 DomainMatrix class 更易于公开访问。
我有一个 4x4 矩阵 A，其每个条目中都有相当长但简单的符号表达式。涉及大约 30 个不同的符号。 "simple" 我的意思是这些符号仅使用 addition/subtraction、multiplication/division 和 整数 幂组合。 "long" 我的意思是，如果我打印出矩阵，它会覆盖三到四个屏幕。
我需要这个矩阵的行列式。或者，更具体地说，我知道行列式是一个特定符号的四阶多项式，我需要这个多项式的系数。 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) 的位置 here, which is necesssary as part of the computation, but end up doing a lot of work because it tries to simplify your very large expressions. It would probably be smarter to only do the simplifications that are needed to cancel terms of the determinant itself. I opened an issue。
另一种加快速度的可能性是 select 一种不同的行列式算法，我不确定是否可行。选项是 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 + '
在某些时候，应该可以将其集成到主矩阵中 class 或者以其他方式使 DomainMatrix class 更易于公开访问。