Numpy：除以 0.5 有什么特别之处？

Question

@Dunes 的 answer 指出，由于流水线化，浮点乘法和除法之间（几乎）没有区别。但是，根据我对其他语言的经验，我预计该部门会变慢。

我的小测试如下：

A=np.random.rand(size)
command(A)

对于不同的命令和size=1e8，我在我的机器上得到以下时间：

Command:    Time[in sec]:
A/=0.5      2.88435101509
A/=0.51     5.22591209412
A*=2.0      1.1831600666
A*2.0       3.44263911247  //not in-place, more cache misses?
A+=A        1.2827270031

最有趣的部分：除以0.5 的速度几乎是除以0.51 的两倍。可以假设，这是由于一些智能优化，例如用A+A 替换除法。但是A*2 和A+A 的时间太远了，无法支持这种说法。

一般来说，值(1/2)^n 的浮点数除法更快：

Size: 1e8
    Command:    Time[in sec]:
    A/=0.5      2.85750007629
    A/=0.25     2.91607499123
    A/=0.125    2.89376401901
    A/=2.0      2.84901714325
    A/=4.0      2.84493684769
    A/=3.0      5.00480890274
    A/=0.75     5.0354950428
    A/=0.51     5.05687212944

如果我们看一下size=1e4，它会变得更有趣：

Command:    1e4*Time[in sec]:
A/=0.5      3.37723994255
A/=0.51     3.42854404449
A*=2.0      1.1587908268
A*2.0       1.19793796539
A+=A        1.11329007149

现在，除以.5 和除以.51 之间没有区别！

我针对不同的 numpy 版本和不同的机器进行了尝试。在某些机器上（例如 Intel Xeon E5-2620）可以看到这种效果，但在其他一些机器上看不到 - 这与 numpy 版本无关。

使用@Ralph Versteegen 的脚本（请参阅他的精彩回答！）我得到以下结果：

• i5-2620 的时序（Haswell，2x6 内核，但不使用 SIMD 的非常旧的 numpy 版本）：

• i7-5500U（Broadwell，2 核，numpy 1.11.2）的时序：

问题是：如果数组大小很大（>10^6 )。

@nneonneo 的回答指出，对于某些英特尔处理器，除以 2 的幂时存在优化，但这并不能解释为什么我们只能看到它对大型阵列的好处。

---

最初的问题是“如何解释这些不同的行为（0.5 与0.51 的划分）？”

这里也是我的原始测试脚本，它产生了时间：

import numpy as np
import timeit

def timeit_command( command, rep):
    print "\t"+command+"\t\t", min(timeit.repeat("for i in xrange(%d):"
        %rep+command, "from __main__ import A", number=7))    

sizes=[1e8,  1e4]
reps=[1,  1e4]
commands=["A/=0.5", "A/=0.51", "A*=2.2", "A*=2.0", "A*2.2", "A*2.0",
          "A+=A", "A+A"]

for size, rep in zip(sizes, reps):
    A=np.random.rand(size)
    print "Size:",size
    for command in commands:
        timeit_command(command, rep)

---

Answer 1

起初我怀疑 numpy 正在调用 BLAS，但至少在我的机器上（python 2.7.13、numpy 1.11.2、OpenBLAS），它没有，正如对 gdb 的快速检查所揭示的那样：

> gdb --args python timing.py
...
Size: 100000000.0
^C
Thread 1 "python" received signal SIGINT, Interrupt.
sse2_binary_scalar2_divide_DOUBLE (op=0x7fffb3aee010, ip1=0x7fffb3aee010, ip2=0x6fe2c0, n=100000000)
    at numpy/core/src/umath/simd.inc.src:491
491 numpy/core/src/umath/simd.inc.src: No such file or directory.
(gdb) disass
   ...
   0x00007fffe6ea6228 <+392>:   movapd (%rsi,%rax,8),%xmm0
   0x00007fffe6ea622d <+397>:   divpd  %xmm1,%xmm0
=> 0x00007fffe6ea6231 <+401>:   movapd %xmm0,(%rdi,%rax,8)
   ...
(gdb) p $xmm1
$1 = {..., v2_double = {0.5, 0.5}, ...}

事实上，无论使用什么常量，numpy 都在运行完全相同的通用循环。所以所有的时间差异都纯粹是由于 CPU 造成的。

实际上，除法是一条执行时间可变的指令。要完成的工作量取决于操作数的位模式，也可以检测和加速特殊情况。根据these tables（我不知道其准确性），在您的 E5-2620（Sandy Bridge）上，DIVPD 具有 10-22 个周期的延迟和反向吞吐量，而 MULPS 具有 10 个周期和 5 个周期的反向吞吐量.

现在，A*2.0 比 A*=2.0 慢。 gdb 显示完全相同的函数被用于乘法，除了现在输出op 与第一个输入ip1 不同。所以它必须纯粹是额外的内存被吸入缓存的人为，减慢了大输入的非就地操作（即使 MULPS 每个周期只产生 2*8/5 = 3.2 字节的输出！）。使用 1e4 大小的缓冲区时，所有内容都适合缓存，因此不会产生显着影响，并且其他开销主要掩盖了 A/=0.5 和 A/=0.51 之间的差异。

不过，这些时间还是有很多奇怪的效果，所以我绘制了一些图表（生成它的代码如下）

我已经根据每个 DIVPD/MULPD/ADDPD 指令的 CPU 周期数绘制了 A 数组的大小。我在 3.3GHz AMD FX-6100 上运行它。黄色和红色垂直线是 L2 和 L3 缓存大小。根据这些表，蓝线是假设的 DIVPD 最大吞吐量，1/4.5 个周期（这似乎是可疑的）。如您所见，即使执行 numpy 操作的“开销”接近于零，即使 A+=2.0 也无法接近这一点。所以大约有 24 个周期的开销只是循环和读取和写入 16 个字节到/从 L2 缓存！非常令人震惊，也许内存访问没有对齐。

许多有趣的效果需要注意：

• 低于 30KB 的数组大部分时间是 python/numpy 中的开销
• 乘法和加法的速度相同（在 Agner 的表格中给出）
• A/=0.5 和A/=0.51 之间的速度差异向图表右侧下降；这是因为当读取/写入内存的时间增加时，它会重叠并掩盖执行除法所需的一些时间。因此，A/=0.5、A*=2.0 和 A+=2.0 成为相同的速度。
• 比较A/=0.51、A/=0.5和A+=2.0之间的最大差异表明除法的吞吐量为4.5-44个周期，与Agner表中的4.5-11不匹配。
• 但是，当 numpy 开销变大时，A/=0.5 和 A/=0.51 之间的差异大部分会消失，尽管仍然存在一些周期差异。这很难解释，因为 numpy 开销无法掩盖进行除法的时间。
• 当远大于 L3 缓存大小时，非就地操作（虚线）会变得异常缓慢，但就地操作不会。它们需要双倍的内存带宽到 RAM，但我无法解释为什么它们会慢 20 倍！
• 虚线在左侧发散。这可能是因为除法和乘法由不同的 numpy 函数处理，具有不同的开销。

不幸的是，在另一台具有不同 FPU 速度、缓存大小、内存带宽、numpy 版本等的 CPU 的机器上，这些曲线可能看起来完全不同。

我的收获是：使用 numpy 将多个算术运算链接在一起会比在 Cython 中执行相同的操作要慢很多倍，同时只对输入进行一次迭代，因为没有“最佳位置”可以让算术运算的成本支配着其他成本。

import numpy as np
import timeit
import matplotlib.pyplot as plt

CPUHz = 3.3e9
divpd_cycles = 4.5
L2cachesize = 2*2**20
L3cachesize = 8*2**20

def timeit_command(command, pieces, size):
    return min(timeit.repeat("for i in xrange(%d): %s" % (pieces, command),
                             "import numpy; A = numpy.random.rand(%d)" % size, number = 6))

def run():
    totaliterations = 1e7

    commands=["A/=0.5", "A/=0.51", "A/0.5", "A*=2.0", "A*2.0", "A+=2.0"]
    styles=['-', '-', '--', '-', '--', '-']

    def draw_graph(command, style, compute_overhead = False):
        sizes = []
        y = []
        for pieces in np.logspace(0, 5, 11):
            size = int(totaliterations / pieces)
            sizes.append(size * 8)  # 8 bytes per double
            time = timeit_command(command, pieces, (4 if compute_overhead else size))
            # Divide by 2 because SSE instructions process two doubles each
            cycles = time * CPUHz / (size * pieces / 2)
            y.append(cycles)
        if compute_overhead:
            command = "numpy overhead"
        plt.semilogx(sizes, y, style, label = command, linewidth = 2, basex = 10)

    plt.figure()
    for command, style in zip(commands, styles):
        print command
        draw_graph(command, style)
    # Plot overhead
    draw_graph("A+=1.0", '-', compute_overhead=True)

    plt.legend(loc = 'best', prop = {'size':9}, handlelength = 3)
    plt.xlabel('Array size in bytes')
    plt.ylabel('CPU cycles per SSE instruction')

    # Draw vertical and horizontal lines
    ymin, ymax = plt.ylim()
    plt.vlines(L2cachesize, ymin, ymax, color = 'orange', linewidth = 2)
    plt.vlines(L3cachesize, ymin, ymax, color = 'red', linewidth = 2)
    xmin, xmax = plt.xlim()
    plt.hlines(divpd_cycles, xmin, xmax, color = 'blue', linewidth = 2)

Answer 2

英特尔 CPU 在除以 2 的幂时进行了特殊优化。例如，请参阅http://www.agner.org/optimize/instruction_tables.pdf，其中声明

FDIV 延迟取决于控制字中指定的精度：64 位精度 给出延迟 38、53 位精度给出延迟 32、24 位精度给出延迟 18。除以 2 的幂需要 9 个时钟。

尽管这适用于 FDIV 而不是 DIVPD（正如 @RalphVersteegen 的回答所指出的那样），但如果 DIVPD 没有实现此优化，那将是相当令人惊讶的。

---

除法通常是一件非常缓慢的事情。但是，除以 2 的幂只是指数移位，尾数通常不需要更改。这使得操作非常快。此外，在浮点表示中很容易检测到 2 的幂，因为尾数将全为零（隐含前导 1），因此这种优化既易于测试又成本低廉。