性能追求第四部分：SIMD 力量与您同在

此时，人们可能想知道 numba（围绕 numpy python 代码的 python 编译器）如何提供比 numpy 更高的性能。为此，让我们单独检查所有三角函数的时序（是的，如果您还记得高中的复杂分析课程，指数和对数就是三角函数！）。但该测试不仅仅与那些想要做高中三角学的人相关：物理引擎，例如游戏中会用到这些函数，机器学习和统计计算大量使用log和exp。因此，正确使用三角函数是实现各种应用的一小步，但却是重要的一步。下表显示了 50m 就地转换的时间：

该表首先体现了性能优势：平方根，一个具有专用 simd 向量化指令的函数，在 numba 和 numpy 中执行的时间完全相同，而所有其他函数的速度都加快了 2 -2.5 时间，表示代码使用 simd 或自动线程自动矢量化。第二条线索是通过使用 ulp（最后一位单位）检查 numba 和 numpy 之间结果的差异来提供的。 ulp 是数值计算准确性的衡量标准，可以使用以下 python 函数轻松计算 numpy 数组：

def compute_ulp_error(array1, array2):
## maxulp set up to a very high number to avoid throwing an exception in the code
    return np.testing.assert_array_max_ulp(array1, array2, maxulp=100000)

这些数值基准表明，平方根函数在 numpy 和 numba 中使用几乎等效的代码，而对于所有其他三角函数，所有 50m 数字的平均值、中值、第 99.9 个百分位数和最大 ulp 值都不同。这是 simd 发挥作用的微妙暗示：矢量化稍微改变了浮点代码的语义以利用关联数学，而浮点数值运算不是关联的。

最后，我们可以使用下面的代码检查 numba 生成函数的代码以获取矢量化汇编指令，如此处详述：

@njit(nogil=true, fastmath=false, cache=true)
def compute_sqrt_with_numba(array):
    np.sqrt(array, array)


@njit(nogil=true, fastmath=false, cache=true)
def compute_sin_with_numba(array):
    np.sin(array, array)


@njit(nogil=true, fastmath=false, cache=true)
def compute_cos_with_numba(array):
    np.cos(array, array)


@njit(nogil=true, fastmath=false, cache=true)
def compute_exp_with_numba(array):
    np.exp(array, array)


@njit(nogil=true, fastmath=false, cache=true)
def compute_log_with_numba(array):
    np.log(array, array)

## check for vectorization
## code lifted from https://tbetcke.github.io/hpc_lecture_notes/simd.html
def find_instr(func, keyword, sig, limit=5):
    count = 0
    for l in func.inspect_asm(func.signatures[sig]).split("\n"):
        if keyword in l:
            count += 1
            print(l)
            if count >= limit:
                break
    if count == 0:
        print("no instructions found")

# compile the function to avoid the overhead of the first call
compute_sqrt_with_numba(np.array([np.random.rand() for _ in range(1, 6)]))
compute_sin_with_numba(np.array([np.random.rand() for _ in range(1, 6)]))
compute_exp_with_numba(np.array([np.random.rand() for _ in range(1, 6)]))
compute_cos_with_numba(np.array([np.random.rand() for _ in range(1, 6)]))
compute_exp_with_numba(np.array([np.random.rand() for _ in range(1, 6)]))

这就是我们如何探测 vmovups 指令是否存在，表明 ymm avx2 寄存器正在用于计算

print("sqrt")
find_instr(compute_sqrt_with_numba, keyword="vsqrtsd", sig=0)
print("\n\n")
print("sin")
find_instr(compute_sin_with_numba, keyword="vmovups", sig=0)

正如我们在下面的输出中看到的，平方根函数使用 x86 simd 矢量化平方根指令 vsqrtsd ，而 sin（以及所有三角函数）使用 simd 指令来访问内存。

sqrt
        vsqrtsd %xmm0, %xmm0, %xmm0
        vsqrtsd %xmm0, %xmm0, %xmm0


sin
        vmovups (%r12,%rsi,8), %ymm0
        vmovups 32(%r12,%rsi,8), %ymm8
        vmovups 64(%r12,%rsi,8), %ymm9
        vmovups 96(%r12,%rsi,8), %ymm10
        vmovups %ymm11, (%r12,%rsi,8)

现在让我们把注意力转移到 perl 和 c 上（毕竟这是一个 perl 博客！）。在第一部分
我们看到，在评估嵌套函数 cos(sin(sqrt(x))) 时，pdl 和 c 给出了相似的性能，并且在第二部分中，单线程 perl 代码与 numba 一样快。但是，如果没有 inline 模块，pdl 和 c 中的各个三角函数又如何呢？在这种情况下将被评估的 c 代码块是：

#pragma omp for simd
  for (int i = 0; i < array_size; i++) {
    double x = array[i];
    array[i] = foo(x);
  }

其中 foo 是 sqrt、sin、cos、log、exp 之一。我们将使用 simd omp pragma 与以下编译标志和 gcc 编译器结合使用，看看我们是否可以让代码拾取提示并自动向量化使用 simd 指令生成的机器代码。

cc = gcc
cflags = -o3 -ftree-vectorize  -march=native -mtune=native -wall -std=gnu11 -fopenmp -fstrict-aliasing -fopt-info-vec-optimized -fopt-info-vec-missed
ldflags = -fpie -fopenmp 
libs =  -lm

在编译过程中，gcc 会通知我们所有错过的优化循环的绝佳机会。下面的性能表也证明了这一点；请注意，标准 c 实现和 pdl 与 numba 等效且性能相同。 perl 通过 pdl 可以在我们的数据科学世界中提供性能。

为了让编译器使用 simd，我们将标志 -o3 替换为 -ofast，这样就不再错过所有这些绝佳的性能机会，并且 c 代码现在可以提供（但有适用于 -ofast 标志的常见警告） 。

WowTo

用AI建立视频知识库

下载

通过这些基准测试，让我们回到最初的 perl 基准测试，并对比通过内联 c 代码的非 simd 感知调用获得的时序：

use inline (
    c    => 'data',
    build_noisy => 1,
    with => qw/alien::openmp/,
    optimize => '-o3 -march=native -mtune=native',
    libs => '-lm'
);

以及 simd 感知型（在下面的代码中，必须包含数学库的矢量化版本才能编译代码）：

use inline (
    c    => 'data',
    build_noisy => 1,
    with => qw/alien::openmp/,
    optimize => '-ofast -march=native -mtune=native',
    libs => '-lmvec'
);

代码的非矢量化版本产生下表

inplace  in  base python took 11.9 seconds
inplace  in pythonjoblib took 4.42 seconds
inplace  in         perl took 2.88 seconds
inplace  in perl/mapcseq took 1.60 seconds
inplace  in    perl/mapc took 1.50 seconds
c array  in            c took 1.42 seconds
vector   in       base r took 1.30 seconds
c array  in   perl/c/seq took 1.17 seconds
inplace  in     pdl - st took 0.94 seconds
inplace in  python numpy took 0.93 seconds
inplace  in python numba took 0.49 seconds
inplace  in   perl/c/omp took 0.24 seconds
c array  in   c with omp took 0.22 seconds
c array  in    c/omp/seq took 0.18 seconds
inplace  in     pdl - mt took 0.16 seconds

而矢量化的这个：

Inplace  in  base Python took 11.9 seconds
Inplace  in PythonJoblib took 4.42 seconds
Inplace  in         Perl took 2.94 seconds
Inplace  in Perl/mapCseq took 1.59 seconds
Inplace  in    Perl/mapC took 1.48 seconds
Vector   in       Base R took 1.30 seconds
Inplace  in     PDL - ST took 0.96 seconds
Inplace in  Python Numpy took 0.93 seconds
Inplace  in Python Numba took 0.49 seconds
C array  in   Perl/C/seq took 0.30 seconds
C array  in            C took 0.26 seconds
Inplace  in   Perl/C/OMP took 0.24 seconds
C array  in   C with OMP took 0.23 seconds
C array  in    C/OMP/seq took 0.19 seconds
Inplace  in     PDL - MT took 0.17 seconds

为了便于与 python 和 r 的各种风格进行比较，我们将之前提供的结果插入到这两个表中。

带回家的要点（其中一些可能有点令人惊讶）：

1. 基础 python 的性能非常糟糕 - 比基础 perl 慢了近 4 倍。即使在并行处理（joblib）下，8 线程 python 代码也比 perl 慢
2. r 在考虑的 3 种动态类型语言中表现最好：比基本 python 快近 9 倍，比基本 perl 快 2.3 倍。
3. 通过访问 c (perl/mapc) 中的数组容器对 perl 数组进行就地修改，缩小了 perl 和 r 之间的差距
4. 考虑到时序的可变性，三个代码、base r、perl/mapc 和 c 数组的等效就地修改大致相同
5. 单线程 pdl 代码提供与 numpy 相同的性能
6. simd 感知代码，例如通过 numba 或（对于本机 c 代码编译器编译指示的情况，例如向量化和非向量化版本的表之间的 c 时序中的 cntract c 数组）提供了最佳的单线程性能。
7. openmp pragmas 确实可以加速 perl 容器的操作（例如映射）。
8. 添加线程（omp 代码）或多线程版本的 pdl 可提供最佳性能。

这些观察产生了以下大问题：

• 观察 1 和 2 令人惊讶的是，基础 python 赢得了数据科学语言的声誉。事实上，凭​​借相当简单的数字代码的这些性能特征，人们应该在数据分析代码中用基本 python 来替代 r（或者对于我们这些使用 perl 的人来说）。
• 由于涉及 c 和 perl 的混合实现可以使用 simd 感知代码提供真正的性能优势，即使使用单个线程，我们是否应该升级 inline::c 代码以使用 simd 友好的编译器标志？或者（对于那些害怕-ofast的人来说，也许是精心准备的内在函数混合学，甚至是汇编？
• 我们是否应该考虑升级 perl 的各种 c/xs 模块，以便它们使用 openmp？
• 为什么没有更多的人使用 pdl 软件这一瑰宝？仅自动线程属性就应该让人们考虑将其用于要求苛刻的数据密集型任务。
• 是否可以创建依赖 pdl 和 openmp/inline::c 进行不同计算的混合应用程序？
• 我们是否应该研究 pdl 的特殊构建，利用 simd 属性和编译器编译指令来允许用户鱼与熊掌兼得（simd 矢量化）和吃掉它（自动线程）？