我试图找到使用 Cython 对 numpy 数组进行水平求和的最快方法。首先，假设我有一个 10 x 100,000 随机浮点数的 2D 数组。我可以创建一个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)
但是，我相信我可以使用data作为memoryview对象的属性访问底层 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 + 缓存未命中的开销。然而，这是一个相当多余的分析的结果，可能还有更微妙的原因。
第二个问题不太有趣：我们将看到，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-arrays，所以它必须使用Python-dispatch才能知道，这a_obj[i]+a_obj[i+1]实际上是两个numpy-arrays的相加。总结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 调度开销理论并不能解释一切。
让我们来看看以下数组的大小：
>>> 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)运行时需要逐行获取信息：每次从内存中将双精度数提取到缓存中时，会预取接下来的 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-cache-misses。
但还有更多乐趣，对象版本可能会更快！
>>> 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-function access_obj：
a[i]并假设它是一个通用的 Python 对象。a[i][j]它必须调用__get_item__的a[i]与Python基础设施的帮助。顺便提一句。因为j必须将其转换为成熟的 Python 整数。__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+i)*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)
它 3比迄今为止最快的版本（12 µs）更快。大型数组的表现如何：
%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"]在设置中）将改进生成的程序集并导致 shape 求和的 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()]
)))
| 归档时间： |
|
| 查看次数： |
1640 次 |
| 最近记录： |