如何在 SymPy 中加速慢速矩阵乘法？答案

【问题标题】：How to accelerate slow matrix multiplication in SymPy?如何在 SymPy 中加速慢速矩阵乘法？
【发布时间】：2014-10-07 23:02:33
【问题描述】：

我正在编写一个工具来使用 SymPy 求解特定的递推方程，并发现其中一个涉及矩阵乘法的步骤花费了非常长的时间。例如，如果我在 iPython 控制台中尝试以下操作，

In [1]: from sympy import *

In [2]: A = Matrix(500, 500, lambda i,j: 2 + abs(i-j) if i-j in [-1, 0, 1] else 0)

In [3]: A[0:7, 0:7]
Out[3]: 
Matrix([
[2, 3, 0, 0, 0, 0, 0],
[3, 2, 3, 0, 0, 0, 0],
[0, 3, 2, 3, 0, 0, 0],
[0, 0, 3, 2, 3, 0, 0],
[0, 0, 0, 3, 2, 3, 0],
[0, 0, 0, 0, 3, 2, 3],
[0, 0, 0, 0, 0, 3, 2]])

In [4]: A * A
...

我从来没有等待足够长的时间来完成它。

我很欣赏符号操作比数值计算慢得多，但这似乎很荒谬。使用 Octave 或其他线性代数包可以在几分之一秒内完成此计算。

有没有人使用 SymPy 的矩阵类来处理大约 1000 行和列的 Rational 条目？

【问题讨论】：

您确定在构建矩阵后需要保持“完全符号化”吗？将A-values 移动到numpy.array 怎么样？提到 Octave 后，无法直接将其与 sympy symbolic-solver 引擎进行比较。在Octavetic(); ... ;toc()-segment 中对您的代码进行基准测试，然后尝试numpy。在对 A[500,500] 个浮点数进行测试时，A*A 在 python (numpy) 中花费了大约 19 [msec]
n.b.当然，可能是我错过了使用“完全符号”符号的要点。你写了关于递归方程的求解器。如果它用于生成三对角矩阵（如演示代码中所示），那肯定是对充分使用这种强大的“全符号”求解器工具的极大误解。
我对 SymPy 很感兴趣，因为它可以精确地处理有理算术，根据需要扩展分母以不丢失精度。在我的研究领域中，必须不通过浮点近似截断解决方案。他们在 Matrix 类中的示例封装了这种 Rational 算术。有没有办法在 SymPy 的 Rational 类中使用 numpy？
说实话，当求解器达到矩阵乘法阶段时，我已经计算出矩阵中的所有符号的数值，所以我只是对任意线性代数感兴趣-精确有理算术。
感谢麦克的笔记。有用于此目的的 pythonic 类。参考。详情如下。

标签： performance matrix multiplication sympy

【解决方案1】：

任意精度是必须的，速度是最好的

有一些用于此目的的 Pythonic 类，可能会对您的任意精度计算感兴趣

decimal.Decimal()
fractions.Fraction( numerator = 0, denominator = 1 )

这些对于精确计算是必不可少的，无论是对于大规模天文学、DEP 模拟还是其他领域，在这些领域，精度不得随着计算策略在长时间/事件/递归以及对标准数的类似威胁中的进展而降低- 表示。

我个人使用decimal.Decimal()（在加密序列随机性分析中的精度高达 5.000.000 位），但并没有专注于将它们“置于”numpy 矩阵中。

Numpy 可以在它的矩阵数据结构（dtype = decimal.Decimal 语法等）中承载这些类型，但是需要进行测试，以验证它的矢量化函数操作和处理这些可爱的 Python 类实例后的整体速度。

性能

作为对标准上 numpy 速度的初步观察（“dense”对于这种规模来说听起来很有趣）2x2 decimal.Decimal-s：

>>> aClk.start();m*m;aClk.stop()
array([[Decimal('0.01524157875323883675019051999'),
        Decimal('5.502209507697009702374335655')],
       [Decimal('1.524157875323883675019051999'),
        Decimal('11.94939027587381419881628113')]], dtype=object)
5732L # 5.7 msec_______________________________________________________________

虽然dtype = numpy.float96 上也是这样

>>> aClk.start();f*f;aClk.stop()
array([[ 0.042788046,  0.74206772],
       [ 0.10081096,  0.46544855]], dtype=float96)
2979L # 2.9 msec_______________________________________________________________

对于 500 x 500 完全填充的 dtype = fractions.Fraction

>>> aClk.start();M*M;aClk.stop()
array([[Fraction(9, 64), Fraction(1, 4), Fraction(64, 25), ...,
        Fraction(64, 81), Fraction(16, 81), Fraction(36, 1)],
        ..,
       [Fraction(1, 1), Fraction(9, 4), Fraction(4, 1), ...,
        Fraction(1, 4), Fraction(25, 36), Fraction(1, 1)]], dtype=object)
2692088L # 2.7 sec_<<<_Fraction_______________________________vs. 19 msec float96

对于 500 x 500 完全填充的密集 dtype = decimal.Decimal

>>> aClk.start();D*D;aClk.stop()
array([[Decimal('0.140625'), Decimal('0.25'), Decimal('2.56'), ...,
        Decimal('0.7901234567901234567901234568'),
        Decimal('0.1975308641975308641975308642'), Decimal('36')],
       [Decimal('3.24'), Decimal('0.25'), Decimal('0.25'), ...,
        Decimal('0.02040816326530612244897959185'), Decimal('0.04'),
        Decimal('0.1111111111111111111111111111')],
       [Decimal('0.1111111111111111111111111111'), Decimal('0.25'),
        Decimal('2.25'), ..., Decimal('0.5102040816326530612244897959'),
        Decimal('0.25'), Decimal('0.0625')],
       ...,
       [Decimal('0'), Decimal('5.444444444444444444444444443'),
        Decimal('16'), ..., Decimal('25'), Decimal('0.81'), Decimal('0.04')],
       [Decimal('1'), Decimal('7.111111111111111111111111113'),
        Decimal('1'), ..., Decimal('0'), Decimal('81'), Decimal('2.25')],
       [Decimal('1'), Decimal('2.25'), Decimal('4'), ..., Decimal('0.25'),
        Decimal('0.6944444444444444444444444444'), Decimal('1')]], dtype=object)
4789338L # 4.8 sec_<<<_Decimal_______________________________vs. 19 msec float96
2692088L # 2.7 sec_<<<_Fraction______________________________vs. 19 msec float96

对 1000x1000 的三对角矩阵的预期可能低于 50 毫秒吗？

由于有 3,000 个非零（稀疏表示）元素，与上面测试的完全填充的 500x500 矩阵中的 250,000 个单元相比，使用这些 Pythonic 类进行任意精度计算会有巨大的性能提升。一旦引擎可以使用numerator/denominator 构造的优势进行 MUL / DIV 操作而不是 Decimal，Fraction 就会有更多空间，但是应该在体内测试确切的性能包络，在您的计算方法使用的实际情况下。

`SparseMatrix` 的 Sympy 语法和 1000x1000 三对角线的测试

@tmyklebu 提出的 1000x1000 SparseMatrix 上的真实测试，由于安装问题，需要更长的时间才能详细说明，但是可能会让您对现实世界的实施项目有一些进一步的了解：

>>> F = sympy.SparseMatrix( 1000, 1000, { (0,0): 1} )        # .Fraction()
>>> D = sympy.SparseMatrix( 1000, 1000, { (0,0): 1} )        # .Decimal()

>>> for i in range( 1000 ):                                  # GEN to have F & D hold
...     for j in range( 1000 ):                              #     SAME values,
...         if i-j in [-1,0,1]:                              # but DIFF representations
...            num = int( 100 * numpy.random.random() )      #     
...            den = int( 100 * numpy.random.random() ) + 1  # + 1 to avoid DIV!0
...            F[i,j] = fractions.Fraction( numerator = num, denominator = den )
...            D[i,j] = decimal.Decimal( str( num ) ) / decimal.Decimal( str( den ) )

# called in Zig-Zag(F*F/D*D/F*F/D*D/...) order to avoid memory-access cache artifacts

>>> aClk.start();VOID=F*F;aClk.stop()
770353L                                      # notice the 1st eval took  TRIPLE LONGER
205585L                                      # notice the 2nd+
205364L # 0.205 sec_<<<_Fraction()____________________________vs. 0.331 sec Decimal()


>>> aClk.start();VOID=D*D;aClk.stop()
383137L # 0.383 sec_<<<_Decimal()____________________________vs. 0.770 sec 1st Fraction()
390164L # 0.390 sec_<<<_Decimal()____________________________vs. 0.205 sec 2nd Fraction()
331291L # 0.331 sec_<<<_Decimal()____________________________vs. 0.205 sec 3rd Fraction()

>>> F[0:4,0:4]
Matrix([
[ 1/52,  6/23,     0,     0],
[42/29, 29/12,     1,     0],
[    0, 57/88, 39/62, 13/57],
[    0,     0, 34/83, 26/95]])
>>> D[0:4,0:4]
Matrix([
[0.0192307692307692, 0.260869565217391,                 0,                 0],
[  1.44827586206897,  2.41666666666667,               1.0,                 0],
[                 0, 0.647727272727273, 0.629032258064516, 0.228070175438596],
[                 0,                 0, 0.409638554216867, 0.273684210526316]])

【讨论】：

感谢您的帮助，我将研究 numpy 和 Fraction。不过，这些性能示例有点毫无意义，因为即使 SymPy 版本在如此小的规模上也能快速执行。我将保留这个问题，以防有人知道如何让 SymPy 的 Matrix/Rational 类执行得更快。
当然，麦克。与 float96 相比，如果 numpy 在 Decimal/Fraction 计算操作上保持恒定的“双倍”开销，它仍然非常快。经过几次测试后，我将向您发布对 numpy 中的 Decimal 和 Fraction 类的一些观察

【解决方案2】：

为什么不使用sympy's sparse matrices 而不是密集矩阵？求解（线性）递归时出现的矩阵通常是稀疏的；一种技术为您提供了一个矩阵，其中第一个上对角线为 1，除最底行（即递归系数所在的位置）外，其他所有位置均为零。

【讨论】：

我用 SymPy 的稀疏实现尝试了一个类似的例子，但仍然需要几秒钟的时间。虽然这是朝着正确方向迈出的一步，但我仍然对是否有一种有效的方法可以让密集矩阵类快速处理数字条目感兴趣。
@mack：我有点怀疑，至少对于您正在查看的尺寸。一个 500x500 * 500x500 的密集矩阵乘法将需要 125M真正复杂的运算——每次乘法的开销将大大超过所进行的算术运算。
我的评论源于我对 Sage 的有限经验。它是 SymPy 的替代品，它也有一个执行任意精度有理算术的矩阵类。与 SymPy 不同，相比之下，Sage 几乎可以立即执行这个 500x500 的小乘法（约 0.5 秒），但我无法忍受 Sage 奇怪的 API。

【解决方案3】：

我没有使用sympy 进行矩阵运算，但我可以重现您在使用此代码时遇到的缓慢情况。看来sympy中的矩阵运算不是那么好。

我会推荐你使用numpy，它有很好的矩阵运算并且速度非常快。这是您在 numpy 中的代码副本，它可以在我的笔记本电脑上在 1 秒内完成乘法运算：

In [1]: import numpy as np

In [2]: A = np.matrix([[2 + abs(i-j) if i-j in [-1, 0, 1] else 0 for i in range(0, 500)] for j in range(0, 500)])

In [3]: A[0:7,0:7]
Out[3]:
matrix([[2, 3, 0, 0, 0, 0, 0],
        [3, 2, 3, 0, 0, 0, 0],
        [0, 3, 2, 3, 0, 0, 0],
        [0, 0, 3, 2, 3, 0, 0],
        [0, 0, 0, 3, 2, 3, 0],
        [0, 0, 0, 0, 3, 2, 3],
        [0, 0, 0, 0, 0, 3, 2]])

In [4]: A * A
Out[4]:
matrix([[13, 12,  9, ...,  0,  0,  0],
        [12, 22, 12, ...,  0,  0,  0],
        [ 9, 12, 22, ...,  0,  0,  0],
        ...,
        [ 0,  0,  0, ..., 22, 12,  9],
        [ 0,  0,  0, ..., 12, 22, 12],
        [ 0,  0,  0, ...,  9, 12, 13]])

【讨论】：

任意精度是必须的，速度是最好的

性能

对 1000x1000 的三对角矩阵的预期可能低于 50 毫秒吗？

SparseMatrix 的 Sympy 语法和 1000x1000 三对角线的测试

`SparseMatrix` 的 Sympy 语法和 1000x1000 三对角线的测试