Numpy矢量化及其算法复杂性

我读过很多关于numpy中的 * vectorized * 代码的文章。我知道一个事实，一个python for循环可能比一个等价的numpy操作慢大约100倍。然而，我认为numpy * vectorization * 的强大之处不仅仅是将一个for循环从Python转换为C/Fortran代码。
我到处读过SIMD（单指令多数据）、BLAS（基本线性代数子程序）和其他低级的东西，但并没有清楚地了解这些库的内部情况，我认为这些库由于在CPU级别进行了并行化，因此能够执行操作，以便以次线性方式进行扩展。
也许举个例子可以帮助我们澄清这一点。假设我们希望计算两个矩阵的乘积，我想检查增加第一个矩阵的行数将如何影响所用时间（这个操作和机器学习有很大关系，如果我们认为行数是 * 批量大小 *，左边的矩阵是数据，右边的矩阵包含模型的参数）.那么，我天真的理解是，在这种情况下，总的运行时间将以次线性的方式缩放（达到某个限制），因此，原则上，只要所有内容都适合RAM，我认为增加bath大小总是一个好主意。
我做了一些基准测试，但情况并不是我所期望的，看起来增长是线性的，而且操作的数量是行数的线性函数，看起来在引擎盖下运行的C/Fortran代码只是在做for循环。
这是密码：

k = 128
N = 100
time_info = {}
for size in list(range(100, 5000, 200)):
    A, B = np.random.random((size, k)), np.random.random((k, k))
    t = time()
    for i in range(N):
        np.dot(A, B)
    delta = time() - t
    time_info[size] = delta / N

x = np.array(list(time_info.keys()))
y = np.array(list(time_info.values())) * 1000
 
# Plotting the Graph
plt.plot(x, y)
plt.title("Elapsed time vs number of rows")
plt.xlabel("Number of rows of left matrix")
plt.ylabel("Time in miliseconds")
plt.show()

看起来趋势是线性的。顺便说一下，我已经检查了np.show_config()，它显示我有openblas。
所以我的问题是：

- 矢量化 * 的确切含义是什么？
基准测试是否合理，是否符合您的预期？
如果是这样，是否有任何值得注意的优化，感谢像SIMD这样的"低级别"的东西，应该会对"矢量化"操作产生影响？或者，也许，它只会在您从非常小的矢量/矩阵到中等大小的矢量/矩阵时产生影响？

最后一个选项在只有当大小非常小时CPU才没有被完全占用时才有意义。例如，如果我说了一些愚蠢的话，请纠正我，如果我们有一个能够并行执行8个数学运算的架构，那么您会期望将一个（1，）向量乘以a（1，）向量将与乘以a一样快（8，）向量乘以a（1，）向量。因此，对于非常小的尺寸，增益将是巨大的但是如果你有上千个元素的向量，那么这种影响将是可以忽略的，时间将线性缩放，因为你将总是CPU被完全占用。这种天真的解释有意义吗？

矢量化的确切含义是什么？
在Python中，特别是Numpy代码中，矢量化是使用原生（C）代码来加速计算，避免Python的开销，主要是CPython解释器的开销。更具体地说，Numpy文档将其定义如下：
NumPy将数组处理交给了C，C中的循环和计算比Python中快得多。为了利用这一点，使用NumPy的程序员消除了Python循环，转而使用数组到数组的操作。向量化可以指C卸载和NumPy代码结构化来利用它。
这并不意味着函数是特别优化的，尽管大多数Numpy函数经常是这样。
基准测试是否合理，是否符合您的预期？
总的来说，是的，这是相当意料之中的。
矩阵乘法是受计算限制的，因此输入的大小不会对性能产生太大影响，因为BLAS使用平铺，因此计算可以有效地使用多级CPU缓存。
SIMD指令适用于所有提供的输入大小。它不仅适用于非常小的矩阵，但您无法使用Numpy来测量，因为使用Numpy从CPython调用BLAS函数的开销远远大于8x8矩阵乘法的计算时间。
OpenBLAS对大于预定义阈值的矩阵使用多线程。该阈值相对较小，因此128 x128矩阵乘法应该已经使用了多线程（至少，这是我在Windows上安装OpenBLAS的机器上的情况）。问题是OpenBLAS不能完美地并行化小矩阵和非常长的矩阵的计算，就像您的基准测试一样。这会导致一些内核在一小部分时间内无法工作，并且计算持续的时间比预期的要长。也就是说，这种开销是一个恒定因素，因此在您的特定基准测试中，它不会受到输入大小的太大影响。您应该会看到，使用更大的方阵时，性能会有小幅提升（相对于输入大小）。或者，您可以尝试其他BLAS实现，如BLIS或MKL。
因此，最后，您应该看到的大多是线性曲线，而在一些具有其他BLAS实现的机器上可能不是这样，或者是因为一些低级别的开销（特别是当线程数随着输入大小而增长时）。
如果是这样，是否有任何值得注意的优化得益于像SIMD这样的低级别的东西，应该对矢量化操作有影响？或者，也许，它只会有影响，当你从非常小的向量/矩阵到中等大小的向量/矩阵？
参见上文。SIMD优化在此处不可见，因为计算是在相对较大的矩阵上进行的，而基准测试是针对不同矩阵大小测量 * 相对 * 时间。SIMD优化确实显著提高了大多数Numpy函数的速度。大多数优化是由编译器执行的（即不是手动的），尽管我们需要帮助编译器这样做。一些特定的函数是手动优化的（正在进行的工作）。
如果您想了解SIMD优化的影响（以及Numpy开销和展开），那么您可以使用不同的矩阵形状来测量整数和所花费的时间：

dtype = np.uint64

a = np.full((2, 16*1024*1024), 1, dtype=dtype)
%timeit -n 10 np.sum(a, axis=0, dtype=dtype)       # 55.7 ms

a = np.full((4, 8*1024*1024), 1, dtype=dtype)
%timeit -n 10 np.sum(a, axis=0, dtype=dtype)       # 41.2 ms

a = np.full((8, 4*1024*1024), 1, dtype=dtype)
%timeit -n 10 np.sum(a, axis=0, dtype=dtype)       # 33.7 ms

a = np.full((16, 2*1024*1024), 1, dtype=dtype)
%timeit -n 10 np.sum(a, axis=0, dtype=dtype)       # 29.2 ms

a = np.full((32, 1*1024*1024), 1, dtype=dtype)
%timeit -n 10 np.sum(a, axis=0, dtype=dtype)       # 26.2 ms

a = np.full((64, 512*1024), 1, dtype=dtype)
%timeit -n 10 np.sum(a, axis=0, dtype=dtype)       # 19.8 ms

a = np.full((128, 256*1024), 1, dtype=dtype)
%timeit -n 10 np.sum(a, axis=0, dtype=dtype)       # 19.9 ms  <-- no speed-up

虽然总的输入大小相同，但性能却有很大的不同。SIMD指令的影响只占开销的一小部分。问题是，您无法独立于其他优化（如从Numpy代码展开）来查看SIMD优化的影响。因此，此基准测试不是非常严格。请注意，这是假设目标函数是矢量化的。这应该只是最近版本的Numpy的情况。
为了获得更好的基准测试，您需要在不使用SIMD指令的情况下重新构建Numpy（实际上，在主流x86/x86-64处理器上，编译器仍然会生成SIMD指令，但只会生成标量指令）。
请注意，大多数Numpy函数都没有使用多线程，事实上，矩阵乘法和其他BLAS函数是AFAIK唯一的多线程函数（由于内部BLAS函数调用，LAPACK函数也可以是多线程的）。

Numpy矢量化及其算法复杂性

1条答案

相关问题

热门标签

最新问答