Широковещательная арифметика NumPy - почему один метод намного эффективнее?
Этот вопрос является продолжением моего ответа в Эффективном способе вычисления матрицы Вандермонда.
Вот настройки:
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, я могу предположить, что ответом будет кеш. Первый метод гораздо более дружественен к кешу, чем второй. Тем не менее, я хотел бы получить официальное слово от кого-то с большим опытом, чем я.
Что может быть причиной такого резкого контраста во времени?
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, y1 а также x2, y2соответственно соответствуют вашим оригинальным тестовым кейсам. ??_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. Более конкретно, для массивов с основными рядами ("C-порядок") конечные измерения являются смежными, в то время как для массивов с основными столбцами ("Fortran-порядок") ведущие измерения являются смежными. Для достаточно больших размеров arr.sum(axis=-1) должно быть быстрее, чем arr.sum(axis=0) для больших рядных массивов дайте или возьмите мелкий шрифт.
Здесь происходит то, что существует огромная разница между этими двумя измерениями (размер 4 и 5000 соответственно), но огромная асимметрия производительности между двумя транспонированными случаями имеет место только для случая широковещания. По общему признанию, у меня сложилось впечатление, что радиовещание использует 0 шагов для построения представлений соответствующего размера. Эти 0 шагов означают, что в более быстром случае доступ к памяти выглядит так долго x массив:
[mem0,mem1,mem2,...,mem4999, mem0,mem1,mem2,...,mem4999, ...] # and so on
где mem* просто обозначает float64 ценность x сидит где-то в оперативке. Сравните это с более медленным случаем, когда мы работаем с формой (5000,4):
[mem0,mem0,mem0,mem0, mem1,mem1,mem1,mem1, mem2,mem2,mem2,mem2, ...]
Мое наивное представление заключается в том, что работа с первым позволяет процессору кэшировать большие фрагменты отдельных значений x в то время, так что производительность отличная. В последнем случае 0 шагов заставляют процессор перескочить на тот же адрес памяти x четыре раза каждый, делая это 5000 раз. Я считаю разумным полагать, что эта установка работает против кэширования, что приводит к общей плохой производительности. Это также согласуется с тем фактом, что смежные случаи не показывают такого снижения производительности: там ЦП должен работать со всеми 5000*4 уникальными float64 значениям и кешированию могут не мешать эти странные операции чтения.
Я тоже пытался смотреть на 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 шагами, даже если оно вначале явно не расширяет массивы. И стоимость связана с формами, и, возможно, какое измерение расширяется.
Я не уверен, что кеширование - это действительно самый влиятельный фактор.
Я также не обученный специалист по информатике, так что, может быть, я ошибаюсь, но позвольте мне провести вас через пару аберраций. Для простоты я использую вызов @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 из архива вы найдете где-нибудь в дереве почти 15000 строк компьютерного кода в файле с именем loops.c. Эти циклы являются самыми внутренними циклами 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 loop\n");
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 %d\n", (int)*count_ptr);
innerloop(dataptr, count_ptr, stride, innerloopdata);
} while (iternext(iter));
NPY_END_THREADS;
innerloop - это внутренний цикл, который мы рассмотрели выше. Сколько накладных расходов идет с последующим?
Для этого нам нужно обратиться к файлу 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 итерациями тела тривиального цикла это будет суверенным.
Я, вероятно, должен профилировать это, но мои навыки недостаточны.