在Cython与NumPy中对int与float求和时的巨大性能差异

在Cython与NumPy中对int与float求和时的巨大性能差异

我不会回答您所有的问题，而只是（在我看来）最有趣的问题。

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

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. 它执行成对求和（种类），从而导致较小的舍入误差。

2. 它以块为单位计算总和（python的代码有些难以理解，这里是相应的模板，其次是使用的函数列表

   pairwise_sum_DOUBLE

   ）

第二点是您要加快速度的原因，其计算类似于以下模式（至少从下面的源代码中可以理解）：

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

这种求和的优点：的结果

a2+a10

不依赖于

a1+a9

并且这两个值都可以在现代CPU上同时计算（例如，流水线化），从而提高了观察的速度。

---

就其价值而言，在我的计算机上cython-integer-sum比numpy的慢。

需要考虑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)

的确，在这种情况下，将内存视图用于双精度不会带来任何提速（不知道为什么），但通常会简化优化器的寿命。例如，将memory-view-
variant与

-ffast-math

compile选项结合使用，这将允许关联性，从而产生与numpy相当的性能：

%%cython -c=-ffast-mathcimport numpy as npdef 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

：

static npy_doublepairwise_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;                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]);        }                res = ((r[0] + r[1]) + (r[2] + r[3])) +   ((r[4] + r[5]) + (r[6] + r[7]));                for (; i < n; i++) { res += (a[i * stride]);        }        return res;    }    else {                npy_uintp n2 = n / 2;        n2 -= n2 % 8;        return pairwise_sum_DOUBLE(a, n2, stride) +    pairwise_sum_DOUBLE(a + n2 * stride, n - n2, stride);    }}