使用 NumPy 对 uint16 与 uint64 数组求和时没有加速？
No speedup when summing uint16 vs uint64 arrays with NumPy?
我必须对相对较小的整数进行大量运算（加法），我开始考虑哪种数据类型在 64 位机器上性能最佳。
我确信将 4 uint16 相加所花费的时间与一个 uint64 相加所花费的时间相同，因为 ALU 仅使用 1 uint64 就可以进行 4 uint16 次相加加法器。 （进位传播意味着这对于单个 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 操作都花费相同的时间（比 float 操作多），没有加速。这是因为 numpy 的 C 实现没有完全优化，还是存在一些我不知道的硬件限制？
TL;DR: 我在Numpy 1.21.1上做了实验分析。实验结果表明 np.sum 没有（真正）使用 SIMD 指令：没有 SIMD 指令用于整数，标量 SIMD 指令用于浮点数！此外，对于较小的整数类型，Numpy 默认将整数转换为 64 位值，以避免溢出!
请注意，这可能无法反映所有 Numpy 版本，因为正在进行的工作是为常用功能提供 SIMD 支持（尚未发布的版本 Numpy 1.22.0rc1 继续这项长期工作）。此外，所使用的编译器或处理器可能会对结果产生重大影响。以下实验是使用从具有 i5-9600KF 处理器的 Debian Linux 上的 pip 检索的 Numpy 完成的。
在 np.sum
的幕后
对于浮点数，Numpy 使用 成对算法，该算法在相对较快的同时非常数值稳定。这可以在代码中看到，但也可以简单地使用分析器：TYPE_pairwise_sum 是在运行时调用以计算总和的 C 函数（其中 TYPE 是 DOUBLE 或 FLOAT ).
对于整数，Numpy 使用经典的朴素归约。在 AVX2 兼容机器上调用的 C 函数是 ULONG_add_avx2。如果类型不是 np.int64.
，它还会令人惊讶地将项目转换为 64 位项目
这里是DOUBLE_pairwise_sum函数
执行的汇编代码的热点部分
3,65 │ a0:┌─→add [=10=]x8,%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 编译代码清楚地使用了 scalar SIMD 指令 vaddsd（仅计算单个双精度项）尽管它成功地展开了循环 8 次。为 FLOAT_pairwise_sum 生成了相同的代码：vaddss 被调用了 8 次。
对于np.uint32，这里是生成的汇编代码的热点部分：
2,37 │160:┌─→add [=11=]x1,%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 指令。它甚至不展开循环。由于对累加器的数据依赖性，执行累加的 add (%rdi),%rdx 指令是这里的瓶颈。对于 np.uint64(despite the name of the function isULONG_add_avx2`).
可以看到相同的循环
但是，np.uint32 版本调用 C 函数 _aligned_contig_cast_uint_to_ulong 以便 将整数项转换为更宽的类型 。 Numpy 这样做是为了 避免整数溢出 。对于类型 np.uint8 和 np.uint16 可以看到同样的事情（尽管函数名称不同）。希望这个函数使用 SIMD 指令 (SSE)，但仍然占用执行时间的很大一部分（约 np.sum 时间的 30%）。
编辑：正如@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 位版本比其他基于整数的版本更快。这是由于转换的开销。由于标量循环和 add 指令对于 8 位、16 位和 32 位整数同样快（对于大多数 64 位主流平台来说这应该是正确的），前两个是同样快的.由于缺少（适当的）循环展开，整数实现比浮点数慢。
由于标量 SIMD 指令，浮点实现同样快。实际上，指令 vaddss（对于 np.float32）和 vaddsd（对于 np.float64）具有相同的延迟和吞吐量（至少在所有现代英特尔处理器上）。因此，这两个实现的吞吐量是相同的，因为这两个实现的循环是相似的（相同的展开）。
补充说明
很遗憾 np.sum 没有充分利用 SIMD 指令，因为这会大大加快使用它的计算速度（尤其是小整数）。
[更新] 查看 Numpy 代码，我发现代码没有向量化，因为 array stride 是一个运行时值并且编译器不会生成步幅为 1 的专用版本。实际上，这在之前的汇编代码中可以部分看出：编译器使用指令 add %r10, %rdi 因为 %r10（步幅目标数组）在编译时是未知的。目前还没有针对 Numpy 代码中这种特定情况下的缩减进行优化（函数相对通用）。这可能会在不久的将来改变。
除了步长问题，一个重要的问题是编译器很难自动向量化代码：浮点加法是not commutative, nor associative（除非使用-ffast-math这样的标志）。
我必须对相对较小的整数进行大量运算（加法），我开始考虑哪种数据类型在 64 位机器上性能最佳。
我确信将 4 uint16 相加所花费的时间与一个 uint64 相加所花费的时间相同，因为 ALU 仅使用 1 uint64 就可以进行 4 uint16 次相加加法器。 （进位传播意味着这对于单个 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 操作都花费相同的时间（比 float 操作多），没有加速。这是因为 numpy 的 C 实现没有完全优化，还是存在一些我不知道的硬件限制？
TL;DR: 我在Numpy 1.21.1上做了实验分析。实验结果表明 np.sum 没有（真正）使用 SIMD 指令：没有 SIMD 指令用于整数，标量 SIMD 指令用于浮点数！此外，对于较小的整数类型，Numpy 默认将整数转换为 64 位值，以避免溢出!
请注意，这可能无法反映所有 Numpy 版本，因为正在进行的工作是为常用功能提供 SIMD 支持（尚未发布的版本 Numpy 1.22.0rc1 继续这项长期工作）。此外，所使用的编译器或处理器可能会对结果产生重大影响。以下实验是使用从具有 i5-9600KF 处理器的 Debian Linux 上的 pip 检索的 Numpy 完成的。
在 np.sum
的幕后
对于浮点数，Numpy 使用 成对算法，该算法在相对较快的同时非常数值稳定。这可以在代码中看到，但也可以简单地使用分析器：TYPE_pairwise_sum 是在运行时调用以计算总和的 C 函数（其中 TYPE 是 DOUBLE 或 FLOAT ).
对于整数，Numpy 使用经典的朴素归约。在 AVX2 兼容机器上调用的 C 函数是 ULONG_add_avx2。如果类型不是 np.int64.
这里是DOUBLE_pairwise_sum函数
3,65 │ a0:┌─→add [=10=]x8,%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 编译代码清楚地使用了 scalar SIMD 指令 vaddsd（仅计算单个双精度项）尽管它成功地展开了循环 8 次。为 FLOAT_pairwise_sum 生成了相同的代码：vaddss 被调用了 8 次。
对于np.uint32，这里是生成的汇编代码的热点部分：
2,37 │160:┌─→add [=11=]x1,%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 指令。它甚至不展开循环。由于对累加器的数据依赖性，执行累加的 add (%rdi),%rdx 指令是这里的瓶颈。对于 np.uint64(despite the name of the function isULONG_add_avx2`).
但是，np.uint32 版本调用 C 函数 _aligned_contig_cast_uint_to_ulong 以便 将整数项转换为更宽的类型 。 Numpy 这样做是为了 避免整数溢出 。对于类型 np.uint8 和 np.uint16 可以看到同样的事情（尽管函数名称不同）。希望这个函数使用 SIMD 指令 (SSE)，但仍然占用执行时间的很大一部分（约 np.sum 时间的 30%）。
编辑：正如@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 位版本比其他基于整数的版本更快。这是由于转换的开销。由于标量循环和 add 指令对于 8 位、16 位和 32 位整数同样快（对于大多数 64 位主流平台来说这应该是正确的），前两个是同样快的.由于缺少（适当的）循环展开，整数实现比浮点数慢。
由于标量 SIMD 指令，浮点实现同样快。实际上，指令 vaddss（对于 np.float32）和 vaddsd（对于 np.float64）具有相同的延迟和吞吐量（至少在所有现代英特尔处理器上）。因此，这两个实现的吞吐量是相同的，因为这两个实现的循环是相似的（相同的展开）。
补充说明
很遗憾 np.sum 没有充分利用 SIMD 指令，因为这会大大加快使用它的计算速度（尤其是小整数）。
[更新] 查看 Numpy 代码，我发现代码没有向量化，因为 array stride 是一个运行时值并且编译器不会生成步幅为 1 的专用版本。实际上，这在之前的汇编代码中可以部分看出：编译器使用指令 add %r10, %rdi 因为 %r10（步幅目标数组）在编译时是未知的。目前还没有针对 Numpy 代码中这种特定情况下的缩减进行优化（函数相对通用）。这可能会在不久的将来改变。
除了步长问题，一个重要的问题是编译器很难自动向量化代码：浮点加法是not commutative, nor associative（除非使用-ffast-math这样的标志）。