numpy 矩阵乘法的奇怪性能结果

Posted 2023-02-16

【中文标题】numpy 矩阵乘法的奇怪性能结果

【英文标题】：Strange performance results for numpy matrix multiplication

【发布时间】：2014-07-09 19:54:49

【问题描述】：

最近我发现了一个案例，其中 numpy 的矩阵乘法显示出非常奇怪的性能（至少对我而言）。为了说明这一点，我创建了一个此类矩阵的示例和一个简单的脚本来演示时序。两者都可以从the repo下载，这里不包含脚本，因为没有数据就没什么用了。

脚本使用dot 函数和einsum 以不同的方式将两对矩阵相乘（每对矩阵在shape 和dtype 方面相同，只是数据不同）。实际上，我注意到了一些异常情况：

第一对 (A * B) 的倍增速度比第二对 (C * D) 快得多。 当我将所有矩阵转换为 float64 时，两对的时间变得相同：比乘以 A * B 的时间长，但比 C * D 的时间短。 einsum（据我了解，numpy 实现）和dot（在我的机器上使用 BLAS）都存在这些效果。 为了完整起见，此脚本在我的笔记本电脑上的输出：

With np.dot:
A * B: 0.142910003662 s
C * D: 4.9057161808 s
A * D: 0.20524597168 s
C * B: 4.20220398903 s
A * B (to float32): 0.156805992126 s
C * D (to float32): 5.11792707443 s
A * B (to float64): 0.52608704567 s
C * D (to float64): 0.484733819962 s
A * B (to float64 to float32): 0.255760908127 s
C * D (to float64 to float32): 4.7677090168 s
With einsum:
A * B: 0.489732980728 s
C * D: 7.34477996826 s
A * D: 0.449800014496 s
C * B: 4.05954909325 s
A * B (to float32): 0.411967992783 s
C * D (to float32): 7.32073783875 s
A * B (to float64): 0.80580997467 s
C * D (to float64): 0.808521032333 s
A * B (to float64 to float32): 0.414498090744 s
C * D (to float64 to float32): 7.32472801208 s

如何解释这样的结果，以及如何更快地将C * D 相乘，例如A * B？

【问题讨论】：

***.com/q/13964606/270986 或多或少是同一个问题，@EricPostpischil 给出了很好的回答。从本质上讲，C*D 的计算达到了float32 类型的次正常范围，并且根据您的处理器，次正常数的运算可能比正常数慢很多倍。在乘法之前转换为float64 时，不再涉及次正规。

另见***.com/a/9314926/270986

@MarkDickinson：感谢您的指点！您能否将其发布为答案，因为它实际上回答了问题的第一部分？在这种情况下，可以做些什么来将非规范化刷新为零？我尝试编译一个 Cython 函数，它使用 gcc 参数 -ffast-math -msse -mfpmath=sse 调用 np.dot，但没有区别。我只关心 BLAS 的点积。

会的。对不起;没有时间更早地写一个深思熟虑的回复。

【参考方案1】：

您看到的减速是由于涉及subnormal numbers 的计算。许多处理器在执行具有次正规输入或输出的算术运算时要慢得多。有几个现有的 *** 问题是相关的：请参阅 this related C# question（尤其是 Eric Postpischil 的 answer）和 this answer 到 C++ 问题以了解更多信息。

在您的特定情况下，矩阵C（dtype 为float32）包含几个次正规数。对于单精度浮点数，次正规/正规边界是2^-126，或在1.18e-38 附近。这是我看到的C：

>>> ((0 < abs(C)) & (abs(C) < 2.0**-126)).sum()  # number of subnormal entries
44694
>>> C.size
682450

因此，大约 6.5% 的 C 条目是次正规的，这足以减慢 C*B 和 C*D 乘法的速度。相比之下，A 和 B 不会靠近次正规边界：

>>> abs(A[A != 0]).min()
4.6801152e-12
>>> abs(B[B != 0]).min()
4.0640174e-07

因此，A*B 矩阵乘法中涉及的中间值都不是次正规的，也不会应用速度损失。

至于你问题的第二部分，我不知道该建议什么。如果您足够努力，并且您使用的是 x64/SSE2（而不是 x87 FPU），则可以从 Python 设置刷新为零和非正规为零标志。请参阅this answer 以获取基于 ctypes 的粗略且不可移植的 hack；如果您真的想遵循这条路线，那么编写一个自定义 C 扩展来执行此操作可能是一个更好的选择。

我很想尝试缩放 C 以使其完全进入正常范围（并将 C*D 中的单个产品也带入正常范围），但如果@ 987654340@ 的值也位于浮点范围的上限。或者，简单地将 C 中的微小值替换为零可能会起作用，但由此产生的精度损失是否显着和/或可接受将取决于您的应用程序。

【讨论】：

我会说，如果可能的话，很难猜测这种性能异常的原因，而没有听说过非正规在硬件上要慢得多（实际上，为什么不默认禁用它们......） .

感谢您提出将 C 和 D 缩放几个数量级的想法 - 在我的情况下这很容易实现，因为那里的值不能大于 1。我自己通过手动将 1e-20 以下的数字刷新为零来解决了这个问题——它也比我需要的精度高得多。但是，缩放似乎更正确方法，所以我会坚持下去。

【参考方案2】：

Mark Dickinson 已经回答了你的问题，但只是为了好玩，试试这个：

Cp = np.array(list(C[:,0]))
Ap = np.array(list(A[:,0]))

这消除了拼接延迟并确保数组在内存中相似。

%timeit Cp * Cp   % 34.9 us per loop
%timeit Ap * Ap   % 3.59 us per loop

哎呀。