Такое поведение обусловлено тем, каким образом numy получает доступ к памяти во время операции Reduce («добавление» является лишь особым случаем) с целью улучшения использования кэша.
Для некоторых случаев (какодин выше), можно применить парное суммирование без большого влияния на производительность.Но в целом принудительное применение этого приведет к огромной потере производительности - возможно, будет проще использовать двойную точность, что в большинстве случаев уменьшит вышеуказанную проблему.
Парное суммирование можно рассматривать как оченьспециальная оптимизация для операции «add», которая выполняется, если соблюдаются некоторые ограничения (подробнее об этом позже).
Суммирование (и многие другие операции сокращения) ограничено пропускной способностью памяти.Жизнь хороша, если мы суммируем по смежной оси: память, извлеченная в кэш для индекса i, будет непосредственно повторно использована для расчета с индексом i+1, i+2, ... без удаления из кэша, док использованию.
Ситуация иная, когда суммирование происходит не вдоль непрерывной оси: для добавления элемента float32 16-float32 извлекаются в кеш, но 15 из них выселяются до того, как они могутбыть использованным и должен быть извлечен снова - что за пустая трата времени.
По этой причине numpy производит суммирование по строкам в этом случае: суммирование первой и второй строк, затем добавление третьей строки к результату, затемчетвертый и тд.Однако парное суммирование осуществляется только для одномерного суммирования и не может использоваться здесь.
Парное суммирование выполняется, когда:
sum вызывается на одно-размерный массив numpy sum вызывается вдоль непрерывной оси
numpy не (пока?) не предлагает способ принудительного парного суммирования без существенного негативного влияния на производительность,
Мой вывод: цель должна состоять в том, чтобы выполнить суммирование вдоль смежной оси, что не только более точно, но и может быть намного быстрее:
A=np.ones((N,2), dtype=np.float32, order="C") #non-contiguous
%timeit A.sum(axis=0)
# 326 ms ± 9.17 ms
B=np.ones((N,2), dtype=np.float32, order="F") # contiguous
%timeit B.sum(axis=0)
# 15.6 ms ± 898 µs
В этом особом случае, когда в строке всего 2 элемента, накладные расходы слишком велики (см. Также объяснение аналогичного поведения здесь ).
Это можно сделать, например, лучшепо-прежнему неточно einsum:
%timeit np.einsum("i...->...", A)
# 74.5 ms ± 1.47 ms
np.einsum("i...->...", A)
# array([16777216., 16777216.], dtype=float32)
или даже:
%timeit np.array([A[:,0].sum(), A[:,1].sum()], dtype=np.float32)
# 17.8 ms ± 333 µs
np.array([A[:,0].sum(), A[:,1].sum()], dtype=np.float32)
# array([17000000., 17000000.], dtype=float32)
, что не только почти так же быстро, как непрерывная версия (наказание за загрузку памяти дважды нетакой же высокой, как загрузка памяти в 16 раз), но также и точный, потому что sum используется для одномерных numpy-массивов.
Для большего числа столбцов разница в смежном случае намного меньше для numpy и einsum-способы:
B=np.ones((N,16), dtype=np.float32, order="F")
%timeit B.sum(axis=0)
# 121 ms ± 3.66 ms
A=np.ones((N,16), dtype=np.float32, order="C")
%timeit A.sum(axis=0)
# 457 ms ± 12.1 ms
%timeit np.einsum("i...->...", A)
# 139 ms ± 651 µs per loop
Но производительность для "точного" трюка очень плохая, вероятно, потому, что латентность больше не может быть скрыта вычислениями:
def do(A):
N=A.shape[1]
res=np.zeros(N, dtype=np.float32)
for i in range(N):
res[i]=A[:,i].sum()
return res
%timeit do(A)
# 1.39 s ± 47.8 ms
Воткровавые детали реализации numpy.
Разницу можно увидеть в коде FLOAT_add с определениями из здесь :
#define IS_BINARY_REDUCE ((args[0] == args[2])\
&& (steps[0] == steps[2])\
&& (steps[0] == 0))
#define BINARY_REDUCE_LOOP(TYPE)\
char *iop1 = args[0]; \
TYPE io1 = *(TYPE *)iop1; \
/** (ip1, ip2) -> (op1) */
#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)
/**begin repeat
* Float types
* #type = npy_float, npy_double, npy_longdouble#
* #TYPE = FLOAT, DOUBLE, LONGDOUBLE#
* #c = f, , l#
* #C = F, , L#
*/
/**begin repeat1
* Arithmetic
* # kind = add, subtract, multiply, divide#
* # OP = +, -, *, /#
* # PW = 1, 0, 0, 0#
*/
NPY_NO_EXPORT void
@TYPE@_@kind@(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func))
{
if (IS_BINARY_REDUCE) {
#if @PW@
@type@ * iop1 = (@type@ *)args[0];
npy_intp n = dimensions[0];
*iop1 @OP@= pairwise_sum_@TYPE@(args[1], n, steps[1]);
#else
BINARY_REDUCE_LOOP(@type@) {
io1 @OP@= *(@type@ *)ip2;
}
*((@type@ *)iop1) = io1;
#endif
}
else if (!run_binary_simd_@kind@_@TYPE@(args, dimensions, steps)) {
BINARY_LOOP {
const @type@ in1 = *(@type@ *)ip1;
const @type@ in2 = *(@type@ *)ip2;
*((@type@ *)op1) = in1 @OP@ in2;
}
}
}
, который когда-то был сгенерирован, выглядит следующим образом:
NPY_NO_EXPORT void
FLOAT_add(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func))
{
if (IS_BINARY_REDUCE) {
#if 1
npy_float * iop1 = (npy_float *)args[0];
npy_intp n = dimensions[0];
*iop1 += pairwise_sum_FLOAT((npy_float *)args[1], n,
steps[1] / (npy_intp)sizeof(npy_float));
#else
BINARY_REDUCE_LOOP(npy_float) {
io1 += *(npy_float *)ip2;
}
*((npy_float *)iop1) = io1;
#endif
}
else if (!run_binary_simd_add_FLOAT(args, dimensions, steps)) {
BINARY_LOOP {
const npy_float in1 = *(npy_float *)ip1;
const npy_float in2 = *(npy_float *)ip2;
*((npy_float *)op1) = in1 + in2;
}
}
}
FLOAT_add может использоваться для одномерного сокращения, в этом случае:
args[0] - указатель на результат / начальное значение (аналогично args[2]) args[1] - входной массив steps[0] и steps[2] равны 0, то есть указатели на скаляр.
и затем можно использовать попарное суммирование (проверяется с помощью IS_BINARY_REDUCE).
FLOAT_add может использоваться для добавления двух векторов, в этом случае:
args[0] первый входмассив args[1] второй входной массив args[2] выходной массив steps - переход от одного элемента к другому в массиве для указанных выше массивов.
Параметр @PW@ равен 1 только для суммирования - для всех других операций парное суммирование не используется.