如何在 SymPy 中加速缓慢的矩阵乘法？

2023-12-04

我正在编写一个工具来使用 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 条目的经验吗？

任意精度是必须的，速度是可有可无的

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

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

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

就我个人而言，我曾与decimal.Decimal()（在密码序列随机性分析中达到 5,000,000 位精度）但并不专注于将它们放在“内部”numpy matrix.

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

表现

作为对标准上 numpy 速度的初步观察（“dense“这个规模听起来很有趣）2x2decimal.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 took

>>> 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在 Decimal 上构建 MUL/DIV 运算，但是应根据您的计算方法使用的实际情况在体内测试确切的性能范围。

Sympy 语法为`SparseMatrix`& 在 1000x1000 三对角线上进行测试

1000x1000 的真实测试SparseMatrix然而，正如 @tmyklebu 所提议的那样，由于安装问题，需要更长的时间才能详细说明可能会让您进一步了解现实世界的实施项目：

>>> 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]])

本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)