Большая разница в производительности при суммировании ints vs floats в Cython vs NumPy

Я суммирую каждый элемент в массиве 1D, используя либо Cython, либо NumPy. При суммировании целых чисел Cython на 20% быстрее. При суммировании поплавков Cython ~ 2.5x медленнее. Ниже приведены две простые функции.

#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 директивы отсутствуют. Зачем?
  • Преобразование целочисленного массива numpy в sum_int в C-указатель (np.int64_t *arr = <np.int64_t *> a.data) повышает производительность на 25%. Делать это для поплавков ничего не делает

Основной вопрос

Как я могу получить такую же производительность в Cython с поплавками, которые я делаю с целыми числами?

EDIT - только подсчет медленный?!?

Я написал еще более простую функцию, которая просто подсчитывает количество итераций. Первый хранит счет как int, второй - как двойной.

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

Я не буду отвечать на все ваши вопросы, но только (в моих глазах) самые интересные.

Начните с вашего примера подсчета:

  1. компилятор способен оптимизировать for-loop в целочисленном случае - результирующий двоичный код не вычисляет ничего - он должен только вернуть значение, предварительно рассчитанное на этапе компиляции.
  2. Это не так для двойного случая, потому что из-за ошибок округления результат будет не 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 использовать его двумя способами:

  1. он выполняет парное суммирование (вид), что приводит к меньшей ошибке округления.

  2. он вычисляет сумму в кусках (код python несколько сложно понять, вот соответствующий шаблон и далее вниз список используемой функции pairwise_sum_DOUBLE)

Второй момент - причина ускорения вашего наблюдения, расчет происходит аналогично следующей схеме (по крайней мере, то, что я понимаю из исходного кода ниже):

a1  + a9 + .....  = r1 
a2  + a10 + ..... = r2
..
a8  + a16 +       = r8

----> sum=r1+....+r8

Преимущество такого рода суммирования: результат a2+a10 не зависит от a1+a9 и эти оба значения могут быть вычислены одновременно (например, конвейерная обработка) на современных процессорах, что приводит к ускорению, которое вы наблюдаете.

Для чего это стоит, на моей машине сумма cython-integer меньше, чем numpy.

Необходимость принятия шага множества numpy-array (который известен только во время выполнения, см. Также этот вопрос об векторизации) предотвращает некоторые оптимизации. Обходным путем является использование представлений памяти, для которых вы можете четко указать, что данные непрерывны, то есть:

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 компиляции -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);
    }
}