SymPy中的符号行列式计算缓慢答案

【问题标题】：Symbolic Determinant Calculation Slow in SymPySymPy中的符号行列式计算缓慢
【发布时间】：2017-09-16 12:49:38
【问题描述】：

我正在进行的一些研究需要象征性地采用大矩阵的行列式；矩阵范围从 18x18 到 318x318。矩阵条目是同一变量omega 中的数值或二次多项式。

目前，我正在尝试在 SymPy 中使用 .det() 方法，但速度很慢；一个 18x18 矩阵现在已经运行了 45 多分钟，并且在我写这篇文章时仍在计算。我意识到行列式计算非常密集，但是我能做些什么来加快速度吗？我已经阅读了Speeding up computation of symbolic determinant in SymPy 的帖子，但没有从帖子中删除任何关于实际可以做什么来加快进程的内容。我能做些什么？

【问题讨论】：

是什么让您认为这是可能的？密集符号 NxN 矩阵的行列式具有math.factorial(N) 项。因此，例如，对于 18x18 矩阵，这是 6,402,373,705,728,000 个术语。我想说我们会在那之前死掉，除非你在那之前很久就会用完 RAM ;-)
@TimPeters：我认为这应该是一个答案。
@TimPeters，我不确定你从哪里得到 N 阶乘，但是有比标准的辅因子展开更有效的算法来计算行列式。例如，使用 LU 分解可以大大减少计算大型矩阵的行列式所需的计算开销。对于 15x15 矩阵，使用辅因子展开计算行列式将需要约 2.3 万亿次乘法。使用 LU 分解将只需要大约 1100 次乘法来计算行列式。必须有一些东西我可以改变以使它运行得更快。
看来你的矩阵有特殊的结构：它不像带有符号 a,b,c,d 的 [[a,b],[c,d]] 而是更像 [[x+1 , x], [x-1, x-3]]，矩阵中有一个符号。这很重要，因为答案确实要小得多，并且可以希望优化计算。请编辑您的问题以描述符号矩阵的结构。
好的，但如果这些完全是 omega 的通用函数，则不能指望行列式有任何简化。它们是欧米茄的线性函数，还是其他次数的多项式？系数是整数还是浮点数？

标签： python matrix sympy determinants

【解决方案1】：

SymPy 对行列式并不天真（请参阅MatrixDeterminant class），但在整个计算过程中处理符号表达式似乎是一个缓慢的过程。当行列式已知为一定程度的多项式时（因为矩阵条目是），结果证明为变量的几个值计算其数值并进行插值会更快。

我的测试用例是一个密集的 15 x 15 矩阵，其中包含变量 omega 的二次多项式，具有整数系数。对于数值行列式，我仍然使用 SymPy 的 .det 方法，因此无论哪种方式，系数最终都是完全相同的长整数。

import numpy as np
from sympy import *
import time
n = 15
omega = Symbol('omega')
A = Matrix(np.random.randint(low=0, high=20, size=(n, n)) + omega*np.random.randint(low=0, high=20, size=(n, n)) + omega**2 * np.random.randint(low=0, high=20, size=(n, n)))
start = time.time()
p1 = A.det()       # direct computation 
print('Time: ' + str(time.time() - start))

start = time.time()
xarr = range(-n, n+1)    # 2*n+1 points to get a polynomial of degree 2*n
yarr = [A.subs(omega, x).det() for x in xarr]  # numeric values
p2 = expand(interpolating_poly(len(xarr), omega, xarr, yarr))  # interpolation
print('Time: ' + str(time.time() - start))

p1 和 p2 都是同一个多项式。运行时间（在相当慢的机器上，来自亚马逊的 t2.nano）：

直接计算需要 74.6 秒，
插值需要 5.4 秒。

如果您的系数是浮点数，并且您在处理它们时不期望精确的算术结果，则可以通过将矩阵评估为 NumPy 数组并使用 NumPy 方法作为行列式来进一步加快速度：

Anum = lambdify(omega, A)
yarr = [np.linalg.det(Anum(x)) for x in xarr]

【讨论】：

【解决方案2】：

作为其他关注此主题的人的后续行动：自从几年前尝试解决这个问题以来，我学到了更多关于数值方法和一般计算的知识，并意识到采用符号行列式是多么不可行矩阵那么大。我最终通过将其转换为特征值问题以数值方式解决了这个问题。故事的寓意...解决问题的方法通常有多种，有些方法可能比其他方法更可行。

【讨论】：