Numpy:划分0.5的特别之处是什么？

Question

这个答案 @Dunes的状态,即由于管道-ING有(几乎)没有浮点乘法和除法之间的差别.但是,从我对其他语言的看法来看,我认为这种划分会更慢.

我的小测试看起来如下:

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核心,但是一个非常古老的numpy版本,不使用SIMD)的时间安排:

• 与i7-5500U(Broadwell,2核,numpy 1.11.2)的时间安排:

问题是:如果阵列大小很大(> 10 ^ 6),那么与某些处理器 的除法0.51相比,除法成本更高的原因是什么？0.5

@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.

实际上,除法是一种具有高度可变执行时间的指令.要完成的工作量取决于操作数的位模式,也可以检测和加速特殊情况.根据这些表(我不知道其准确性),在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缓存大小.根据这些表,1/4.5周期(看起来可疑),蓝线是DIVPD的假设最大吞吐量.正如你所看到的,即使A+=2.0执行numpy操作的"开销"接近于零,也不会接近这一点.因此,大约有24个周期的开销只是循环,并从L2缓存读取和写入16个字节!非常令人震惊,也许内存访问没有对齐.

有很多有趣的效果需要注意:

• 在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),因此这种优化既易于测试又便宜实现.