为什么numpy的求和比"+"运算符慢10倍？

Question

我发现很奇怪的是，np.sum比手写的求和慢了10倍。 应用axis的np.sum：

p1 = np.random.rand(10000, 2)
def test(p1):
    return p1.sum(axis=1)
%timeit test(p1)

每次循环平均耗时186微秒±4.21微妙，共7次循环，每次循环执行1000次

未指定轴的np.sum：

p1 = np.random.rand(10000, 2)
def test(p1):
    return p1.sum()
%timeit test(p1)

每次循环平均耗时为17.9微秒，标准差为236纳秒(基于7次运行，每次运行10000次循环)

+:

p1 = np.random.rand(10000, 2)
def test(p1):
    return p1[:,0] + p1[:,1]
%timeit test(p1)

每次循环耗时15.8微秒±328纳秒（7次运行的平均值和标准差，每次循环运行100000次）

乘法：

p1 = np.random.rand(10000, 2)
def test(p1):
    return p1[:,0]*p1[:,1]
%timeit test(p1)

每个循环的平均时间为15.7微秒，标准偏差为701纳秒（7次运行，每次10000次循环）

我不知道这是为什么。有什么想法吗？我的numpy版本是1.15.3。

编辑： 使用10000000：

np.sum (with axis): 202 ms (5 x)
np.sum (without axis): 12 ms
+ : 46 ms (1 x)
* : 44.3 ms 

所以我猜测玩耍会有一些开销，到某种程度上...

- beesleep

---

只有一小部分的开销可能与成对求和实现有关。不过只是很小的一部分 - 就我所知，prod 不会进行成对操作，而且在我的测试中，prod 的运行时间约为 sum 的 5/6。我认为 NumPy 也在使用 SIMD 进行 + 而不是 sum，但我还不确定。 - user2357112

你的“乘法”做了一些不同的事情...其他人只是使用p1，并基本忽略了p2。 - Sam Mason

这是一对求和例程的源代码，可以在此处找到：https://github.com/numpy/numpy/blob/v1.15.3/numpy/core/src/umath/loops.c.src#L1662。被求和的元素数量足够小，以至于该例程应立即进入直接非成对循环情况，但代码路径仍似乎比用于像“prod”之类的东西的代码路径具有一些开销。正如先前所述，相对于“+”，这只是开销的一小部分。 - user2357112

@beesleep 浮点数中乘法和加法并没有太大的区别；如果有的话，乘法可能会更容易一些。当然，对于整数来说情况是不同的。 - Cubic

这很奇怪，主要是因为轴应该返回非连续元素的内存视图（这几乎不是优化代码，因为你在这里非常不友好地使用了缓存）。实际上，我可以通过将p1 = np.random.rand(10000, 2)更改为p1 = np.random.rand(2, 10000)和p1.sum(axis=1)更改为p1.sum(axis=0)来改变代码性能约10倍。 - Alex Huszagh

Answer 1

主要区别在于计算a.sum(axis=1)时的开销更大。计算缩减（在此情况下为sum）不是一件琐碎的事情：

• 必须考虑舍入误差，因此使用成对求和来减少误差。
• 对于更大的数组，平铺很重要，因为它可以充分利用可用的缓存。
• 为了能够使用现代CPU的SIMD指令/乱序执行能力，应同时计算多行。

我在此处和此处更详细地讨论了上述主题。

但是，如果只有两个元素需要相加，则所有这些都是不必要的，也不如朴素求和好 - 您将获得相同的结果，但开销更小且速度更快。

对于仅有1000个元素的情况，调用NumPy功能的开销可能比实际执行这1000次加法（或乘法）还要高 - 如您所见，对于10^4，运行时间仅大约增加了2倍，这是开销对于10^3起更大作用的明确标志!在此回答中更详细地研究了开销和缓存未命中的影响。

让我们查看分析器结果，看看上述理论是否成立（我使用perf）：

对于a.sum(axis=1)：

  17,39%  python   umath.cpython-36m-x86_64-linux-gnu.so       [.] reduce_loop
  11,41%  python   umath.cpython-36m-x86_64-linux-gnu.so       [.] pairwise_sum_DOUBLE
   9,78%  python   multiarray.cpython-36m-x86_64-linux-gnu.so  [.] npyiter_buffered_reduce_iternext_ite
   9,24%  python   umath.cpython-36m-x86_64-linux-gnu.so       [.] DOUBLE_add
   4,35%  python   python3.6                                   [.] _PyEval_EvalFrameDefault
   2,17%  python   multiarray.cpython-36m-x86_64-linux-gnu.so  [.] _aligned_strided_to_contig_size8_src
   2,17%  python   python3.6                                   [.] lookdict_unicode_nodummy
   ...

使用 reduce_loop 和 pairwise_sum_DOUBLE 的开销很大。

对于 a[:,0]+a[:,1]：

   7,24%  python   python3.6                                   [.] _PyEval_EvalF
   5,26%  python   python3.6                                   [.] PyObject_Mall
   3,95%  python   python3.6                                   [.] visit_decref
   3,95%  python   umath.cpython-36m-x86_64-linux-gnu.so       [.] DOUBLE_add
   2,63%  python   python3.6                                   [.] PyDict_SetDef
   2,63%  python   python3.6                                   [.] _PyTuple_Mayb
   2,63%  python   python3.6                                   [.] collect
   2,63%  python   python3.6                                   [.] fast_function
   2,63%  python   python3.6                                   [.] visit_reachab
   1,97%  python   python3.6                                   [.] _PyObject_Gen

预料之中：Python开销扮演了重要角色，简单的DOUBLE_add被使用。

---

调用a.sum()时，开销较小。

• 一次只调用reduce_loop一次，而不是为每一行都调用一次，这意味着显著减少开销。
• 没有创建新的结果数组，不再需要将1000个双精度浮点数写入内存。

因此可以预期，a.sum()更快（尽管必须执行2000个加法而不是1000个 - 但正如我们所看到的，大部分时间都花在开销和实际工作上，而不是加法产生的负担）。

---

通过运行获得的数据:

perf record python run.py
perf report

以及

#run.py
import numpy as np
a=np.random.rand(1000,2)

for _ in range(10000):
  a.sum(axis=1)
  #a[:,0]+a[:,1]

- ead

---

有趣的是，这也让我有点惊讶。我本来期望会有一些开销，但没有想到会这么大。如果轴大小低于某个阈值，几乎让人觉得值得实现一个特殊情况的微小缩减。 - Eelco Hoogendoorn

我同意这感觉很令人惊讶。我注意到使用linspace比arange产生了10倍的开销。 - beesleep

Answer 2

关于.sum()函数的axis参数和无参数的区别，使用axis参数需要生成一个与输入长度相同的浮点数数组，每行都有一个元素。这意味着它必须沿axis=1调用reduce() 10,000次。而不使用axis参数，则将所有元素的总和计算为单个浮点数，只需通过数组的平面表示进行一次reduce调用。

我不确定为什么手动添加函数更快，也不想深入挖掘源代码，但我认为我有一个很好的猜测。我认为开销来自于它必须对每行执行axis=1的reduce操作，因此需要进行10,000次单独的reduce调用。在手动添加函数中，当定义“+”函数的参数时，仅执行一次轴拆分，然后可以并行地将拆分列的每个元素相加。

- mevers303

---

虽然猜测不错，但 NumPy 实际上不需要复制浮点数数组。NumPy 有一个广泛的内存视图系统，它允许它使用指向原始缓冲区、起始位置、长度和步幅的引用计数指针来存储数组。为了创建轴或进行高级索引，永远不需要数据复制... 您可以通过检查 a = np.arange(100); b = a[::-1]; a.ctypes.data == b.ctypes.data 来验证这一点，这将检查数组 a 和 b 的起始指针是否相同，并打印 true。 - Alex Huszagh

您也可以通过相同的代码进行验证，只需设置一个值：a = np.arange(100); b = a[::-1]; a[4] = 5; b[2]，这将打印出5，尽管该值通常为4。 - Alex Huszagh

所以开销并非来自于拆分数组，但加法函数确实能够并行地将数组相加，是吗？这比大量调用 reduce 要更快。 - mevers303

说实话，我得检查一下实现才能确定，但如果NumPy是那样做的话，那会让我很惊讶。这似乎极其低效，因为它可以仅仅加倍步幅，将1个指针增加8，调用高度优化的C加程序，而且在两种情况下都可以轻松完成（显然，它必须做更多的工作，因为索引（特别是花式索引）和处理潜在的多维数组是相当复杂的）。 - Alex Huszagh

等等...你说得对。它在底层调用了reduce函数。虽然不是10,000次，但是...我的想法被颠覆了... https://github.com/numpy/numpy/blob/9d0225b800df3c2b0bffee9960df87b15527509c/numpy/core/src/multiarray/calculation.c#L519 - Alex Huszagh