SymPy 中的符号行列式计算缓慢
Symbolic Determinant Calculation Slow in SymPy
我正在进行的一些研究需要象征性地取大矩阵的行列式;矩阵范围从 18x18 到 318x318。矩阵条目是同一变量中的数值或二次多项式,omega
.
目前,我正在尝试使用 SymPy 中的 .det()
方法,但速度很慢;一个 18x18 矩阵现在已经 运行 超过 45 分钟,并且在我写这篇文章时仍在计算。我意识到行列式计算非常密集,但我能做些什么来加快速度吗?
我已经阅读了 中的 post,但并没有从 post 中摘取任何关于实际可以做些什么来加快进程的内容。我能做什么?
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]
作为对查看此线程的其他人的跟进:自从几年前尝试解决这个问题以来,我学到了更多关于数值方法和一般计算的知识,并意识到采用符号行列式是多么不可行一个大的矩阵。我最终通过将其转换为特征值问题以数值方式解决了这个问题。这个故事的寓意...通常有多种解决问题的方法,有些方法可能比其他方法更可行。
我正在进行的一些研究需要象征性地取大矩阵的行列式;矩阵范围从 18x18 到 318x318。矩阵条目是同一变量中的数值或二次多项式,omega
.
目前,我正在尝试使用 SymPy 中的 .det()
方法,但速度很慢;一个 18x18 矩阵现在已经 运行 超过 45 分钟,并且在我写这篇文章时仍在计算。我意识到行列式计算非常密集,但我能做些什么来加快速度吗?
我已经阅读了
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]
作为对查看此线程的其他人的跟进:自从几年前尝试解决这个问题以来,我学到了更多关于数值方法和一般计算的知识,并意识到采用符号行列式是多么不可行一个大的矩阵。我最终通过将其转换为特征值问题以数值方式解决了这个问题。这个故事的寓意...通常有多种解决问题的方法,有些方法可能比其他方法更可行。