Большая разница в производительности при суммировании чисел и чисел в Cython против NumPy
Я суммирую каждый элемент в одномерном массиве, используя либо Cython, либо NumPy. При суммировании целых чисел Cython работает на ~20% быстрее. При суммировании чисел с запятой Cython работает примерно в 2,5 раза медленнее. Ниже приведены две простые функции.
#cython: boundscheck=False
#cython: wraparound=False
def sum_int(ndarray[np.int64_t] a):
cdef:
Py_ssize_t i, n = len(a)
np.int64_t total = 0
for i in range(n):
total += a[i]
return total
def sum_float(ndarray[np.float64_t] a):
cdef:
Py_ssize_t i, n = len(a)
np.float64_t total = 0
for i in range(n):
total += a[i]
return total
Задержки
Создайте два массива по 1 миллиону элементов каждый:
a_int = np.random.randint(0, 100, 10**6)
a_float = np.random.rand(10**6)
%timeit sum_int(a_int)
394 µs ± 30 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit a_int.sum()
490 µs ± 34.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit sum_float(a_float)
982 µs ± 10.8 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit a_float.sum()
383 µs ± 4.42 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
Дополнительные баллы
- NumPy превосходит (с большим отрывом) результаты с плавающей точкой и даже превосходит собственную целочисленную сумму.
- Разница в производительности для
sum_floatто же самое сboundscheckа такжеwraparoundдирективы отсутствуют. Зачем? - Преобразование массива целых чисел в
sum_intна указатель C (np.int64_t *arr = <np.int64_t *> a.data) повышает производительность еще на 25%. Делать это для поплавков ничего не делает
Главный вопрос
Как я могу получить ту же производительность в Cython с плавающей точкой, что и с целыми числами?
РЕДАКТИРОВАТЬ - Просто считать медленный?!?
Я написал еще более простую функцию, которая просто подсчитывает количество итераций. Первый хранит счет как int, последний как double.
def count_int():
cdef:
Py_ssize_t i, n = 1000000
int ct=0
for i in range(n):
ct += 1
return ct
def count_double():
cdef:
Py_ssize_t i, n = 1000000
double ct=0
for i in range(n):
ct += 1
return ct
Сроки подсчета
Я запускал их только один раз (боюсь кеширования). Не знаю, если цикл на самом деле выполняется для целого числа, но count_double имеет ту же производительность, что и sum_float сверху. Это безумие...
%timeit -n 1 -r 1 count_int()
1.1 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)
%timeit -n 1 -r 1 count_double()
971 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)
1 ответ
Я не собираюсь отвечать на все ваши вопросы, а только (на мой взгляд) самые интересные.
Давайте начнем с вашего примера подсчета:
- компилятор может оптимизировать цикл for в целочисленном случае - результирующий двоичный файл ничего не вычисляет - он только должен вернуть значение, предварительно рассчитанное на этапе компиляции.
- Это не относится к двойному случаю, потому что из-за ошибок округления результат не будет
1.0*10**6и потому что Cython компилируется в IEEE 754 (не-ffast-math) режим по умолчанию.
Об этом следует помнить при взгляде на свой код Cython: компилятору не разрешается переставлять суммирования (IEEE 754), и поскольку результат первых суммирований необходим для следующего, есть только одна длинная строка, в которой все операции Подождите.
Но самое важное понимание: numpy делает не то же самое, что ваш код на Cython:
>>> sum_float(a_float)-a_float.sum()
2.9103830456733704e-08
Да, никто не сказал numpy (в отличие от вашего кода Cython), что сумма должна быть рассчитана следующим образом
((((a_1+a2)+a3)+a4)+...
И Numpy воспользоваться этим двумя способами:
он выполняет парное суммирование (вид), что приводит к меньшей ошибке округления.
он вычисляет сумму в блоках (код на python немного сложен для понимания, вот соответствующий шаблон и далее список используемых функций
pairwise_sum_DOUBLE)
Второй момент - причина ускорения, которое вы наблюдаете, расчет происходит аналогично следующей схеме (по крайней мере, как я понимаю из исходного кода ниже):
a1 + a9 + ..... = r1
a2 + a10 + ..... = r2
..
a8 + a16 + = r8
----> sum=r1+....+r8
Преимущество такого рода суммирования: результат a2+a10 не зависит от a1+a9 и эти оба значения могут быть рассчитаны одновременно (например, конвейерная обработка) на современных процессорах, что приводит к ускорению, которое вы наблюдаете.
Что бы это ни стоило, на моей машине cython-integer-sum медленнее, чем numpy's.
Необходимость учета шага массива numpy (который известен только во время выполнения, см. Также этот вопрос о векторизации) предотвращает некоторые оптимизации. Обходной путь должен использовать представления памяти, для которых вы можете ясно дать понять, что данные являются непрерывными, то есть:
def sum_int_cont(np.int64_t[::1] a):
Что приводит к значительному ускорению на моей машине (фактор 2):
%timeit sum_int(a_int)
2.64 ms ± 46.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit sum_int_cont(a_int)
1.31 ms ± 19 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit a_int.sum()
2.1 ms ± 105 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Это правда, что в этом случае использование просмотров памяти для удвоений не приводит к ускорению (не знаю почему), но в целом упрощает жизнь оптимизатора. Например, объединение варианта просмотра памяти с -ffast-math Параметры компиляции, которые допускают ассоциативность, приводят к производительности, сравнимой с numpy:
%%cython -c=-ffast-math
cimport numpy as np
def sum_float_cont(np.float64_t[::1] a):
cdef:
Py_ssize_t i, n = len(a)
np.float64_t total = 0
for i in range(n):
total += a[i]
return total
И сейчас:
>>> %timeit sum_float(a_float)
3.46 ms ± 226 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
>>> %timeit sum_float_cont(a_float)
1.87 ms ± 44 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
>>> %timeit a_float.sum()
1.41 ms ± 88.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
Листинг pairwise_sum_DOUBLE:
/*
* Pairwise summation, rounding error O(lg n) instead of O(n).
* The recursion depth is O(lg n) as well.
* when updating also update similar complex floats summation
*/
static npy_double
pairwise_sum_DOUBLE(npy_double *a, npy_uintp n, npy_intp stride)
{
if (n < 8) {
npy_intp i;
npy_double res = 0.;
for (i = 0; i < n; i++) {
res += (a[i * stride]);
}
return res;
}
else if (n <= PW_BLOCKSIZE) {
npy_intp i;
npy_double r[8], res;
/*
* sum a block with 8 accumulators
* 8 times unroll reduces blocksize to 16 and allows vectorization with
* avx without changing summation ordering
*/
r[0] = (a[0 * stride]);
r[1] = (a[1 * stride]);
r[2] = (a[2 * stride]);
r[3] = (a[3 * stride]);
r[4] = (a[4 * stride]);
r[5] = (a[5 * stride]);
r[6] = (a[6 * stride]);
r[7] = (a[7 * stride]);
for (i = 8; i < n - (n % 8); i += 8) {
r[0] += (a[(i + 0) * stride]);
r[1] += (a[(i + 1) * stride]);
r[2] += (a[(i + 2) * stride]);
r[3] += (a[(i + 3) * stride]);
r[4] += (a[(i + 4) * stride]);
r[5] += (a[(i + 5) * stride]);
r[6] += (a[(i + 6) * stride]);
r[7] += (a[(i + 7) * stride]);
}
/* accumulate now to avoid stack spills for single peel loop */
res = ((r[0] + r[1]) + (r[2] + r[3])) +
((r[4] + r[5]) + (r[6] + r[7]));
/* do non multiple of 8 rest */
for (; i < n; i++) {
res += (a[i * stride]);
}
return res;
}
else {
/* divide by two but avoid non-multiples of unroll factor */
npy_uintp n2 = n / 2;
n2 -= n2 % 8;
return pairwise_sum_DOUBLE(a, n2, stride) +
pairwise_sum_DOUBLE(a + n2 * stride, n - n2, stride);
}
}