广播 NumPy 算术 - 为什么一种方法的性能如此之高

Question

这个问题是我在 计算范德蒙德矩阵的有效方法。

这是设置：

x = np.arange(5000)  # an integer array
N = 4

现在，我将以两种不同的方式计算范德蒙德矩阵：

m1 = (x ** np.arange(N)[:, None]).T

和

m2 = x[:, None] ** np.arange(N)

健全性检查：

np.array_equal(m1, m2)
True

这些方法是相同的，但它们的性能不是：

%timeit m1 = (x ** np.arange(N)[:, None]).T
42.7 µs ± 271 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit m2 = x[:, None] ** np.arange(N)
150 µs ± 995 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

因此，第一种方法尽管最后需要换位，但仍比第二种方法快 3 倍以上。

唯一的区别是，在第一种情况下，广播的阵列较小，而在第二种情况下，它是较大的阵列。

因此，对numpy的工作原理有相当不错的了解，我可以猜测答案 将涉及缓存。第一种方法对缓存更加友好 比第二个。但是，我想要一个官方的话 比我更有经验。

时间上出现这种鲜明对比的原因可能是什么？

Answer 1

我也试着看broadcast_arrays：

In [121]: X,Y = np.broadcast_arrays(np.arange(4)[:,None], np.arange(1000))
In [122]: timeit X+Y
10.1 µs ± 31.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [123]: X,Y = np.broadcast_arrays(np.arange(1000)[:,None], np.arange(4))
In [124]: timeit X+Y
26.1 µs ± 30.6 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [125]: X.shape, X.strides
Out[125]: ((1000, 4), (4, 0))
In [126]: Y.shape, Y.strides
Out[126]: ((1000, 4), (0, 4))

np.ascontiguousarray将 0 步距尺寸转换为完整尺寸

In [132]: Y1 = np.ascontiguousarray(Y)
In [134]: Y1.strides
Out[134]: (16, 4)
In [135]: X1 = np.ascontiguousarray(X)
In [136]: X1.shape
Out[136]: (1000, 4)

使用完整阵列的操作速度更快：

In [137]: timeit X1+Y1
4.66 µs ± 161 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

因此，使用 0 步长数组会受到某种时间损失，即使它没有先显式扩展数组。 成本与形状有关，可能还与扩展的维度有关。

Answer 2

我不相信缓存真的是这里最有影响力的因素。

我也不是一个训练有素的计算机科学家，所以我很可能错了，但让我带你了解几个问题。 为简单起见，我使用 @hpaulj 的调用，即"+"显示与"**"基本相同的效果。

我的工作催眠是，我认为外循环的开销比连续的可矢量化最内层循环要贵得多。

因此，让我们首先尽量减少重复的数据量，因此缓存不太可能产生太大影响：

>>> from timeit import repeat
>>> import numpy as np
>>> 
>>> def mock_data(k, N, M):
...     x = list(np.random.randint(0, 10000, (k, N, M)))
...     y = list(np.random.randint(0, 10000, (k, M)))
...     z = list(np.random.randint(0, 10000, (k, N, 1)))
...     return x, y, z
...   
>>> k, N, M = 500, 5000, 4
>>>
>>> repeat('x.pop() + y.pop()', setup='x, y, z = mock_data(k, M, N)', globals=globals(), number=k)
[0.017986663966439664, 0.018148145987652242, 0.018077059998176992]
>>> repeat('x.pop() + y.pop()', setup='x, y, z = mock_data(k, N, M)', globals=globals(), number=k)
[0.026680009090341628, 0.026304758968763053, 0.02680662798229605]

在这里，两种方案都有连续的数据，添加次数相同，但具有 5000 次外部迭代的版本要慢得多。当我们恢复缓存时，尽管在试验中，差异大致相同，但比率变得更加明显：

>>> repeat('x[0] + y[0]', setup='x, y, z = mock_data(k, M, N)', globals=globals(), number=k)
[0.011324503924697638, 0.011121788993477821, 0.01106808998156339]
>>> repeat('x[0] + y[0]', setup='x, y, z = mock_data(k, N, M)', globals=globals(), number=k)
[0.020170683041214943, 0.0202067659702152, 0.020624138065613806]

回到最初的"外部和"情景，我们看到在不连续的长维情况下，情况变得更糟。由于我们不必读取比连续场景中更多的数据，因此这不能用未缓存的数据来解释。

>>> repeat('z.pop() + y.pop()', setup='x, y, z = mock_data(k, M, N)', globals=globals(), number=k)
[0.013918839977122843, 0.01390116906259209, 0.013737019035033882]
>>> repeat('z.pop() + y.pop()', setup='x, y, z = mock_data(k, N, M)', globals=globals(), number=k)
[0.0335254140663892, 0.03351909795310348, 0.0335453050211072]

此外，两者都从跨试验缓存中受益：

>>> repeat('z[0] + y[0]', setup='x, y, z = mock_data(k, M, N)', globals=globals(), number=k)
[0.012061356916092336, 0.012182610924355686, 0.012071475037373602]
>>> repeat('z[0] + y[0]', setup='x, y, z = mock_data(k, N, M)', globals=globals(), number=k)
[0.03265167598146945, 0.03277428599540144, 0.03247103898320347]

从卡奇主义者的角度来看，这充其量是不确定的。

因此，让我们看一下来源。 从压缩包构建当前的 NumPy 后，您会在树中的某个位置找到名为"loops.c"的文件中近 15000 行计算机生成的代码。这些循环是ufuncs的最内层循环，与我们的情况最相关的部分似乎是

#define BINARY_LOOP
char *ip1 = args[0], *ip2 = args[1], *op1 = args[2];
npy_intp is1 = steps[0], is2 = steps[1], os1 = steps[2];
npy_intp n = dimensions[0];
npy_intp i;
for(i = 0; i < n; i++, ip1 += is1, ip2 += is2, op1 += os1)
/*
* loop with contiguous specialization
* op should be the code working on `tin in1`, `tin in2` and
* storing the result in `tout * out`
* combine with NPY_GCC_OPT_3 to allow autovectorization
* should only be used where its worthwhile to avoid code bloat
*/
#define BASE_BINARY_LOOP(tin, tout, op) 
BINARY_LOOP { 
const tin in1 = *(tin *)ip1; 
const tin in2 = *(tin *)ip2; 
tout * out = (tout *)op1; 
op; 
}
etc.

在我们的例子中，有效负载似乎足够精简，特别是如果我正确解释有关连续专业化和自动矢量化的评论。现在，如果我们只进行 4 次迭代，开销与有效载荷的比率开始看起来有点麻烦，而且它并没有就此结束。

在文件ufunc_object.c中，我们找到以下代码片段

/*
* If no trivial loop matched, an iterator is required to
* resolve broadcasting, etc
*/
NPY_UF_DBG_PRINT("iterator loopn");
if (iterator_loop(ufunc, op, dtypes, order,
buffersize, arr_prep, arr_prep_args,
innerloop, innerloopdata) < 0) {
return -1;
}
return 0;

实际循环如下所示

NPY_BEGIN_THREADS_NDITER(iter);
/* Execute the loop */
do {
NPY_UF_DBG_PRINT1("iterator loop count %dn", (int)*count_ptr);
innerloop(dataptr, count_ptr, stride, innerloopdata);
} while (iternext(iter));
NPY_END_THREADS;

内环路是我们上面看到的内部环路。迭代下一个会带来多少开销？

为此，我们需要转到文件nditer_templ.c.src，在那里我们找到了

/*NUMPY_API
* Compute the specialized iteration function for an iterator
*
* If errmsg is non-NULL, it should point to a variable which will
* receive the error message, and no Python exception will be set.
* This is so that the function can be called from code not holding
* the GIL.
*/
NPY_NO_EXPORT NpyIter_IterNextFunc *
NpyIter_GetIterNext(NpyIter *iter, char **errmsg)
{
etc.

此函数返回一个函数指针，指向预处理所做的事情之一

/* Specialized iternext (@const_itflags@,@tag_ndim@,@tag_nop@) */
static int
npyiter_iternext_itflags@tag_itflags@_dims@tag_ndim@_iters@tag_nop@(
NpyIter *iter)
{
etc.

解析这超出了我的范围，但无论如何，它是一个函数指针，必须在外循环的每次迭代中调用，据我所知，函数指针不能内联，所以与 4 次迭代相比，一个琐碎的循环体这将是悬而未决的。

我可能应该介绍这一点，但我的技能不足。

Answer 3

虽然我担心我的结论不会比你的结论更基本("可能是缓存")，但我相信我可以通过一组更本地化的测试来帮助集中我们的注意力。

考虑您的示例问题：

M,N = 5000,4
x1 = np.arange(M)
y1 = np.arange(N)[:,None]
x2 = np.arange(M)[:,None]
y2 = np.arange(N)
x1_bc,y1_bc = np.broadcast_arrays(x1,y1)
x2_bc,y2_bc = np.broadcast_arrays(x2,y2)
x1_cont,y1_cont,x2_cont,y2_cont = map(np.ascontiguousarray,
[x1_bc,y1_bc,x2_bc,y2_bc])

如您所见，我定义了一堆数组来进行比较。x1、y1y2和x2分别对应于原始测试用例。??_bc对应于这些数组的显式广播版本。它们与原始数据共享数据，但它们具有明确的 0 步长以获得适当的形状。最后，??_cont是这些广播阵列的连续版本，就像用np.tile构造一样。

所以x1_bc、y1_bc、x1_cont和y1_cont都有形状(4, 5000)，但前两者是零步长，后两者是连续数组。出于所有意图和目的，获取这些相应数组对中的任何一个的幂都应该给我们相同的连续结果(正如 hpaulj 在评论中指出的那样，换置本身本质上是免费的，所以我将在下面忽略最外层的转置)。

以下是与您的原始支票相对应的时间：

In [143]: %timeit x1 ** y1
...: %timeit x2 ** y2
...: 
52.2 µs ± 707 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
96 µs ± 858 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

以下是显式广播数组的计时：

In [144]: %timeit x1_bc ** y1_bc
...: %timeit x2_bc ** y2_bc
...: 
54.1 µs ± 906 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
99.1 µs ± 1.51 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each

同样的事情。这告诉我，这种差异并不是由于从索引表达式到广播数组的过渡造成的。这主要是意料之中的，但检查一下也无妨。

最后，连续数组：

In [146]: %timeit x1_cont ** y1_cont
...: %timeit x2_cont ** y2_cont
...: 
38.9 µs ± 529 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
45.6 µs ± 390 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

很大一部分差异消失了！

那我为什么要检查这个呢？有一个一般的经验法则，如果您在 python 中对大型尾随维度使用矢量化操作，则可以使用 CPU 缓存。更具体地说，对于行主("C 顺序")数组，尾随维度是连续的，而对于列主("Fortran-order")数组，前导维度是连续的。对于足够大的尺寸arr.sum(axis=-1)应该比arr.sum(axis=0)行主 numpy 数组给出或采取一些小印件的速度快。

这里发生的事情是两个维度(分别为尺寸 4 和 5000)之间存在巨大差异，但两个转置情况之间的巨大性能不对称只发生在广播案例中。我诚然挥手的印象是广播使用 0 步幅来构建适当大小的视图。这些 0 步长意味着在较快的情况下，长x数组的内存访问如下所示：

[mem0,mem1,mem2,...,mem4999, mem0,mem1,mem2,...,mem4999, ...] # and so on

其中mem*只是表示位于 RAM 中某处的xfloat64值。将此与我们正在处理形状(5000,4)的较慢情况进行比较：

[mem0,mem0,mem0,mem0, mem1,mem1,mem1,mem1, mem2,mem2,mem2,mem2, ...]

我天真的想法是，使用前者允许CPU一次缓存x单个值的更大块，因此性能非常好。在后一种情况下，0 步幅使 CPU 在相同的内存地址上跳来跳去，每次x四次，执行 5000 次。我发现有理由相信此设置对缓存有效，从而导致整体性能不佳。这也与连续的情况没有显示这种性能影响的事实一致：CPU必须使用所有5000 * 4的唯一float64值，并且缓存可能不会受到这些奇怪读取的阻碍。