在Cython中对数组求和
我试图找到最快的方式来水平总结使用Cython的数组的一个numpy数组。 为了让事情开始，我们假设我有一个随机浮点数为10 x 100,000的二维数组。 我可以创建一个object数组，每列都作为数组中的一个值，如下所示：
n = 10 ** 5
a = np.random.rand(10, n)
a_obj = np.empty(n, dtype='O')
for i in range(n):
a_obj[i] = a[:, i]
我想要做的是找到每一行的总和。 它们都可以这样简单计算：
%timeit a.sum(1)
414 µs ± 24.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit a_obj.sum()
113 ms ± 7.01 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
对象数组速度慢了250倍。
在尝试总结之前，我想要时间访问每个元素。 当直接遍历每个项目时，Cython无法加速访问对象数组中的每个成员：
def access_obj(ndarray[object] a):
cdef int i
cdef int nc = len(a)
cdef int nr = len(a[0])
for i in range(nc):
for j in range(nr):
a[i][j]
%timeit access_obj(a_obj)
42.1 ms ± 665 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
但是，我相信我可以使用作为memoryview对象的data属性访问底层C数组，并获得更快的访问权限：
def access_obj_data(ndarray[object] a):
cdef int i
cdef int nc = len(a)
cdef int nr = len(a[0])
cdef double **data = <double**>a.data
for i in range(nc):
for j in range(nr):
data[i][j]
我认为这在计时时缓存结果，所以我不得不这样做一次。
%timeit -n 1 -r 1 access_obj_data(a_obj)
8.17 µs ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)
问题是你不能只将一个a.data投入到一个C数组中。 如果我打印出第一列，我会得到以下结果。
5e-324
2.1821467023e-314
2.2428810855e-314
2.1219957915e-314
6.94809615162997e-310
6.94809615163037e-310
2.2772194067e-314
2.182150145e-314
2.1219964234e-314
0.0
也有类似的问题，但没有任何明确的答案可以用代码运行快速完成此操作。
我想有两个不同的问题：
a_obj.sum()这么慢？ 我们将看到a_obj.sum()慢度可以追踪到Python-dispatch + cache未命中的开销。 然而，这是一个相当多余的分析的结果，可能会有更微妙的原因。
第二个问题不太有趣：我们将看到，Cython对阵列中的对象不够了解，无法加速该功能。
还有一个问题：“什么是最快的方式？”。 答案是“这取决于您的数据和硬件”，因为几乎总是出现这种问题。
不过，我们将利用已获得的见解对现有版本进行略微改进。
我们先从“ a_obj.sum()慢”开始。
作为你的例子在我的机器上的基准（是的，因子200 ）：
>>> %timeit a.sum(1)
510 µs ± 14.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
>>> %timeit a_obj.sum()
90.3 ms ± 494 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
首先在a_obj.sum()发生了什么？ 简化算法如下所示：
result=a_obj[0]+a_obj[1]+a_obj[2]+...
然而，该算法并不知道a_obj[i]是numpy数组，所以它必须使用Python-dispatch来知道， a_obj[i]+a_obj[i+1]实际上是两个numpy数组的加法。 这是总共10双打的开销，我们必须支付10**5次。 如果您对更多细节感兴趣，我前一段时间回答了有关Python-dispatch成本的问题。
如果我们的数组的形状为10**5x10 ，那么我们将只支付10倍的开销，与10**5增加相比，它不会很大：
>>> b=np.random.rand(n,10)
>>> b_obj=to_obj_array(b) # see listing at the end for code
>>> %timeit b.sum(1)
2.43 ms ± 13.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
>>> %timeit b_obj.sum()
5.53 ms ± 17.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
正如你所看到的，现在的因素（大约2）要小得多。 然而，Python-dispatch-overhead理论无法解释所有的事情。
我们来看看下面这些数组的大小：
>>> c=np.random.rand(10**4,10**4)
>>> c_obj=to_obj_array(c)
>>> %timeit c.sum(1)
52.5 ms ± 354 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
>>> %timeit c_obj.sum()
369 ms ± 2.53 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
这一次的因素大约是8！ 原因是obj版本中的缓存未命中。 这是一项内存限制任务，所以基本上内存速度是瓶颈。
在numpy数组中，元素连续存储在内存中，一行一行。 每次从内存中获取双精度数据时，都会获取7个双精度数据，并将这8个精度数据保存在高速缓存中。
因此，当c.sum(1)运行时，信息需要按行进行：每次从内存中将双c.sum(1)提取到高速缓存中时，接下来的7个数字都会被预取，因此可以进行总结。 这是最有效的内存访问模式。
不同的是对象版本。 这些信息是需要列的。 因此，每当我们向缓存中获取额外的7双打时，我们并不需要它，并且在我们可以处理这些双倍时，它们已经从缓存中被逐出（因为缓存对于列中的所有数字来说都不够大）以及需要再次从慢速存储器中取出。
这意味着采用第二种方法时，每个双精度数据将从内存中读取8次，或者从内存/高速缓存未命中中读取8次。
我们可以很容易地验证使用valgrind的cachegrind来衡量这两种变体的缓存未命中情况（请参见test_cache.py`末尾列表）：
valgrind --tool=cachegrind python test_cache.py array
...
D1 misses: 144,246,083
...
valgrind --tool=cachegrind python test_cache.py object
...
D1 misses: 1,285,664,578
...
即c_obj -version约有8倍以上的D1缓存未命中。
但更有趣的是，对象版本可能会更快！
>>> d=np.random.rand(800,6)
>>> d_obj=to_obj_array(d)
>>> %timeit d.sum(1)
20.2 µs ± 217 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
>>> %timeit d_obj.sum()
12.4 µs ± 155 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
那怎么可能？ 他们可能会很快，因为
但速度较慢的版本快两倍？
我不是100％肯定的，但这是我的理论：对于数组版本，我们计算第i行的总和：
...
sum_i=sum_i+a[i][j]
sum_i=sum_i+a[i][j+1]
sum_i=sum_i+a[i][j+2]
...
所以硬件不能并行执行指令，因为前一个操作的结果需要计算下一个。
但是对于对象版本，计算如下所示：
...
sum_i=sum_i+a[i][j]
sum_k=sum_k+a[k][j]
sum_m=sum_m+a[i][j]
...
这些步骤可以由CPU并行完成，因为它们之间没有依赖关系。
那么，为什么cython如此缓慢呢？ 因为cython对数组中的对象一无所知，所以它并不比通常很慢的Python代码更快！
让我们分解一下cython函数access_obj ：
a[i]一无所知，并且认为它是一个通用的Python对象。 a[i][j]它必须在Python基础结构的帮助下调用a[i] __get_item__ 。 顺便说一句。 因为那个j必须被转换成一个完整的Python-integer。 __get_item__从存储的c-double构造一个完整的Python-float并将其返回给cython。 在底线之下，我们使用Python-dispatch并创建了两个临时的Python对象 - 您可以看到所需的成本。
加快总和
我们已经看到，总结是一个内存限制任务（至少在我的机器上），所以最重要的是避免缓存未命中（这最重要的还可能意味着，例如，我们需要不同的C- order和F-order numpy-arrays）。
但是，如果我们的任务必须共享一个具有CPU高负荷任务的CPU，那么它可能会再次成为事实上的CPU限制，因此它有一些优点可以专注于不仅仅是现金缺失。
在后面将介绍的C版本中，我们希望将两种方法的优良特性结合起来：
这是一个简化的版本，它处理8行的行：
//fast.c
#define MAX 8
void sum_rows_blocks(double *orig, double *out, int n_rows, int n_cols){
int n_blocks=n_rows/MAX;
for (int b=0; b<n_blocks; b++){ //for every block
double res[MAX]={0.0}; //initialize to 0.0
for(int j=0;j<n_cols;j++){
for(int i=0;i<MAX;i++){ //calculate sum for MAX-rows simultaniously
res[i]+=orig[(b*MAX+i)*n_cols+j];
}
}
for(int i=0;i<MAX;i++){
out[b*MAX+i]=res[i];
}
}
//left_overs:
int left=n_rows-n_blocks*MAX;
double res[MAX]={0.0}; //initialize to 0.0
for(int j=0;j<n_cols;j++){
for(int i=0;i<left;i++){ //calculate sum for left rows simultaniously
res[i]+=orig[(n_blocks*MAX)*n_cols+j];
}
}
for(int i=0;i<left;i++){
out[n_blocks*MAX+i]=res[i];
}
}
现在我们可以使用cython来包装这个函数：
//c_sum.pyx
cimport numpy as np
import numpy as np
cdef extern from "fast.c":
void sum_rows_blocks(double *orig, double *out, int n_rows, int n_cols)
def sum_rows_using_blocks(double[:,:] arr):
cdef int n_rows
cdef int n_cols
cdef np.ndarray[np.double_t] res
n_rows=arr.shape[0]
n_cols=arr.shape[1]
res=np.empty((n_rows, ),'d')
sum_rows_blocks(&arr[0,0], &res[0], n_rows, n_cols)
return res
在构建扩展后（请参阅setup.py列表），我们可以进行比较：
import c_sum
%timeit c_sum.sum_rows_using_blocks(d) #shape=(800,6)
4.26 µs ± 95.8 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
它比现在最快的版本（12微秒）要快3倍。 它如何为大型阵列提供支持：
%timeit c_sum.sum_rows_using_blocks(c) #shape=(10**4,10**4)
49.7 ms ± 600 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
好，它比略快c.sum(1)这需要52毫秒，这意味着我们没有额外的缓存未命中。
可以改进吗？ 我确实这样认为：例如，如果不使用默认编译标志-fwrapv （通过在设置中包含extra_compile_args = ["-fno-wrapv"] ）将会改进所得到的组件并且导致形状总和的10％加速(800, 6) 。
人数：
def to_obj_array(a):
n=a.shape[1]
obj=np.empty((n,),dtype='O')
for i in range(n):
obj[i]=a[:,i]
return obj
check_cache.py ：
import sys
import numpy as np
def to_obj_array(a):
n=a.shape[1]
obj=np.empty((n,),dtype='O')
for i in range(n):
obj[i]=a[:,i]
return obj
c=np.random.rand(10**4,10**4)
c_obj=to_obj_array(c)
for i in range(10):
if sys.argv[1]=="array":
c.sum(1)
elif sys.argv[1]=="object":
c_obj.sum()
else:
raise Exception("unknown parameter")
setup.py ：
#setup.py
from distutils.core import setup, Extension
from Cython.Build import cythonize
import numpy as np
setup(ext_modules=cythonize(Extension(
name='c_sum',
sources = ["c_sum.pyx"],
#extra_link_args = ["fast.o"],
include_dirs=[np.get_include()]
)))