我必须对相对较小的整数进行大量操作(加法),我开始考虑哪种数据类型在 64 位机器上提供最佳性能。
我确信,将4个<code>uint16</code>相加所需的时间与一个<code<uint64</code〕相加所需时间相同,因为ALU只需使用1个<code>uint64>加法器即可进行4个<code>uint16>/code>加法运算。(进位传播意味着这对于单个64位加法器来说不那么容易,但整数SIMD指令就是这样工作的。)
显然情况并非如此:
In [3]: data = np.random.rand(10000)
In [4]: int16 = data.astype(np.uint16)
In [5]: int64 = data.astype(np.uint64)
In [6]: int32 = data.astype(np.uint32)
In [7]: float32 = data.astype(np.float32)
In [8]: float64 = data.astype(np.float64)
In [9]: %timeit int16.sum()
13.4 µs ± 43.3 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [10]: %timeit int32.sum()
13.9 µs ± 347 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [11]: %timeit int64.sum()
9.33 µs ± 47.8 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [12]: %timeit float32.sum()
5.79 µs ± 6.51 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [13]: %timeit float64.sum()
6 µs ± 3.54 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
所有 int 操作都花费相同的时间(比浮点操作更多),没有加速。这是由于 numpy 的 C 实现没有完全优化还是存在一些我不知道的硬件限制?
TL;DR:我在Numpy 1.21.1上做了一个实验分析。实验结果表明< code>np.sum并没有(真正)利用SIMD指令:没有SIMD指令用于整数,标量SIMD指令用于浮点数!此外,默认情况下,Numpy将整数转换为较小整数类型的64位值,以避免溢出!
请注意,这可能不会反映所有的Numpy版本,因为目前正在为常用函数提供SIMD支持(尚未发布的版本Numpy 1.22.0rc1继续这项长期工作)。此外,使用的编译器或处理器可能会显著影响结果。下面的实验是使用从带有i5-9600KF处理器的Debian Linux上的pip中检索的Numpy完成的。
对于浮点数,Numpy使用一种成对算法,这种算法在速度相对较快的同时,在数值上非常稳定。这可以在代码中看到,也可以简单地使用一个分析器:< code>TYPE_pairwise_sum是在运行时被调用来计算总和的C函数(其中< code>TYPE是< code>DOUBLE或< code>FLOAT)。
对于整数,Numpy 使用经典的朴素约简。调用的 C 函数在 AVX2 兼容计算机上ULONG_add_avx2
。如果类型不是 np.int64
,它还会令人惊讶地将项目转换为 64 位项目。
下面是由< code>DOUBLE_pairwise_sum函数执行的汇编代码的热部分
3,65 │ a0:┌─→add $0x8,%rcx ; Loop iterator
3,60 │ │ prefetcht0 (%r8,%rax,1) ; Prefetch data
9,46 │ │ vaddsd (%rax,%rbp,1),%xmm1,%xmm1 ; Accumulate an item in xmm1[0]
4,65 │ │ vaddsd (%rax),%xmm0,%xmm0 ; Same for xmm0[0]
6,91 │ │ vaddsd (%rax,%rdx,1),%xmm4,%xmm4 ; Etc.
7,77 │ │ vaddsd (%rax,%rdx,2),%xmm7,%xmm7
7,41 │ │ vaddsd (%rax,%r10,1),%xmm2,%xmm2
7,27 │ │ vaddsd (%rax,%rdx,4),%xmm6,%xmm6
6,80 │ │ vaddsd (%rax,%r11,1),%xmm3,%xmm3
7,46 │ │ vaddsd (%rax,%rbx,1),%xmm5,%xmm5
3,46 │ │ add %r12,%rax ; Switch to the next items (x8)
0,13 │ ├──cmp %rcx,%r9 ; Should the loop continue?
3,27 │ └──jg a0 ; Jump to the beginning if so
Numpy编译的代码显然使用标量SIMD指令vaddsd
(仅计算单个双精度项),尽管它成功地展开了8次循环。为FLOAT_pairwise_sum
生成相同的代码:vaddss
调用8次。
对于 np.uint32
,下面是生成的程序集代码的热门部分:
2,37 │160:┌─→add $0x1,%rax ; Loop iterator
95,95 │ │ add (%rdi),%rdx ; Accumulate the values in %rdx
0,06 │ │ add %r10,%rdi ; Switch to the next item
│ ├──cmp %rsi,%rax ; Should the loop continue?
1,08 │ └──jne 160 ; Jump to the beginning if so
Numpy显然不会对np.uint32
类型使用SIMD指令。它甚至不会展开循环。由于对累加器的数据依赖性,执行累加的添加 (%rdi),%rdx
指令是这里的瓶颈。在 np.uint64 中可以看到相同的循环(尽管函数的名称是
ULONG_add_avx2')。
然而,np.uint32
版本调用C函数_aligned_contig_cast_uint_to_long
,以便将整数项转换为更宽的类型。Numpy这样做是为了避免整数溢出。对于np.uint8
和np.uint16
类型,可以看到相同的情况(尽管函数的名称不同)。希望该函数使用了SIMD指令(SSE),但仍占用了很大一部分执行时间(约占np.sum
时间的30%)。
EDIT:正如@user2357112supportsMonica所指出的,可以显式指定np.sum
的dtype
参数。当它与输入数组的dtype
匹配时,则不执行转换。这会导致较小的执行时间,但可能会导致溢出。
以下是我的机器上的结果:
uint16: 7.17 µs ± 80 ns per loop (mean ± std. dev. of 7 runs, 20000 loops each)
uint32: 7.11 µs ± 12.3 ns per loop (mean ± std. dev. of 7 runs, 20000 loops each)
uint64: 5.05 µs ± 8.57 ns per loop (mean ± std. dev. of 7 runs, 20000 loops each)
float32: 2.88 µs ± 9.27 ns per loop (mean ± std. dev. of 7 runs, 20000 loops each)
float64: 3.06 µs ± 10.6 ns per loop (mean ± std. dev. of 7 runs, 20000 loops each)
首先,并不是说结果与问题中提供的结果非常相似,这意味着在我的机器上看到的行为可以在另一台机器上成功复制。因此,解释也应该是一致的。
如您所见,64位版本比其他基于整数的版本更快。这是由于转换的开销。前两个同样快,因为标量循环和< code>add指令对于8位、16位和32位整数同样快(对于大多数64位主流平台应该是这样)。整数实现比浮点实现慢,因为缺少(适当的)循环展开。
由于采用了标量 SIMD 指令,浮点实现速度同样快。事实上,指令 vaddss(对于 np.float32)和 vaddsd
(对于
np.float64
)具有相同的延迟和吞吐量(至少在所有现代英特尔处理器上)。因此,两个实现的吞吐量是相同的,因为两个实现的循环是相似的(相同的展开)。
遗憾的是,np.sum
没有充分利用 SIMD 指令,因为这会大大加快使用它的计算速度(尤其是小整数)。
[更新]查看 Numpy 代码,我发现代码没有矢量化,因为数组步幅是运行时值,编译器不会生成步幅为 1 的专用版本。实际上,这可以在前面的汇编代码中部分看到:编译器使用指令添加 %r10,%rdi,因为 %r10
(目标数组的步幅)在编译时是未知的。目前还没有针对 Numpy 代码中这种减少的特定情况进行优化(函数相对通用)。这种情况在不久的将来可能会改变。
除了跨距问题之外,还有一个很大的问题使得编译器很难自动向量化代码:浮点加法不是可交换的,也不是关联的(除非使用了像< code>-ffast-math这样的标志)。