python - 在 Cython 与 NumPy 中对整数与 float 求和时性能差异很大

Question

我正在使用 Cython 或 NumPy 对一维数组中的每个元素求和。对整数求和时，Cython 快了约 20%。对 float 求和时，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

时间

创建两个数组，每个数组包含 100 万个元素:

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 在 float 方面的表现优于(相当大的优势)，甚至超过了它自己的整数和。
• sum_float 的性能差异与boundscheck相同和 wraparound缺少指令。为什么？
• 转换 sum_int 中的整数 numpy 数组指向 C 指针 (np.int64_t *arr = <np.int64_t *> a.data) 可将性能提高 25%。对花车这样做没有任何作用

主要问题

如何在 Cython 中使用 float 获得与使用整数相同的性能？

编辑 - 只是计数很慢？!？

我写了一个更简单的函数，它只计算迭代次数。第一个将计数存储为 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)

Answer 1

我不会回答您所有的问题，只会回答(在我看来)最有趣的问题。

让我们从您的计数示例开始:

1. 编译器能够在整数情况下优化 for 循环 - 生成的二进制文件不计算任何东西 - 它只需要返回在编译阶段预先计算的值。
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. 它执行 pairwise summation(某种程度上)，这会导致更小的舍入误差。

2. 它以 block 的形式计算总和(python 的代码有点难以理解，这里是 corresponding template 以及使用函数 pairwise_sum_DOUBLE 的列表)

   <

第二点是您观察到的加速的原因，计算类似于以下模式(至少我从下面的源代码中理解):

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

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

这种求和的优点:a2+a10的结果不依赖于a1+a9，这两个值可以同时计算(例如pipelining ) 在现代 CPU 上，这会导致您观察到的加速。

---

就其值(value)而言，在我的机器上，cython-integer-sum 比 numpy 的慢。

需要考虑 numpy 数组的步幅(这仅在运行时已知，另请参阅 this question 关于矢量化)阻止了一些优化。一种解决方法是使用内存 View ，您可以明确指出数据是连续的，即:

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)

的确，在这种情况下，对 double 使用内存 View 不会带来任何加速(不知道为什么)，但总的来说它简化了优化器的生命周期。例如，将 memory-view-variant 与 -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);
    }
}