我想有两个不同的问题：
- Why is
a_obj.sum()太慢了？
- 为什么 cython 函数代码很慢？
我们将看到缓慢的a_obj.sum()可以跟踪Python调度+缓存未命中的开销。然而，这是一个相当多余的分析结果，可能还有更微妙的原因。
第二个问题不太有趣：我们会看到，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 调度来知道，a_obj[i]+a_obj[i+1]实际上是两个 numpy 数组的加法。总结起来这是相当大的开销10我们必须支付的费用加倍10**5次。我回答了一个question前段时间关于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 个双精度值保存在缓存中。
So when c.sum(1)运行时需要逐行获取信息：每次从内存中将双精度值提取到缓存中时，都会预取接下来的 7 个数字，以便对它们进行求和。这是最有效的内存访问模式。
对象版本不同。这些信息需要按列进行。所以每次我们获取额外的7将双精度数据放入缓存中，我们还不需要它，当我们可以处理这些双精度数据时，它们已经从缓存中逐出（因为缓存不足以容纳列中的所有数字），并且需要从慢速数据中获取再次记忆。
这意味着使用第二种方法，每个 double 将从内存中读取 8 次，或者从内存/缓存未命中的读取次数将增加 8 倍。
我们可以使用 valgrind 轻松验证缓存研磨测量两个变体的缓存未命中（请参阅 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
...
i.e c_obj-version 的 D1 缓存未命中次数大约是原来的 8 倍。
但还有更多乐趣，对象版本可以更快！
>>> 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)
怎么可能呢？它们可能很快相似，因为
- 整个数组可以保存在缓存中，因此不会出现额外的缓存未命中情况。
- python-dispatch 的开销是可以承受的。
但对于其他较慢的版本来说，速度要快两倍？
我不是 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:
- Cython 一无所知
a[i]并假设它是一个通用的 Python 对象。
- 这意味着，对于
a[i][j]它必须调用__get_item__ of a[i]在Python基础设施的帮助下。顺便提一句。为了那个原因j必须转换为完整的 Python 整数。
-
__get_item__从存储的 c-double 构造一个成熟的 Python-float 并将其返回给 cython。
- 强制转换将 Python-float 返回为 c-float。
在底线下方，我们使用了 Python 调度并创建了两个临时 Python 对象 - 您可以看到所需的成本。
加快求和速度
我们已经看到，求和是一个内存限制任务（至少在我的机器上），所以最重要的是避免缓存未命中（这一点也可能意味着，例如，我们需要不同的 C 算法）阶和 F 阶 numpy 数组）。
然而，如果我们的任务必须与 CPU 密集型任务共享一个 CPU，它可能会再次成为事实上的 CPU 密集型任务，因此它有一些优点，可以不仅仅关注现金缺失。
在稍后将介绍的 C 版本中，我们希望结合这两种方法的良好特性：
- 与对象版本类似，我们希望并行计算多个行总和，以便能够使用 CPU 固有的并行性。
- 与 numpy-array-sum 类似，我们希望缓存未命中次数最少，因此我们不会同时计算所有行，而只计算一小块（例如 8）。这种按块计算的优点是，一旦缓存的数据在真正需要之前就不会从缓存中逐出。
下面是一个简化版本，它以 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)这需要52ms，这意味着我们没有额外的缓存未命中。
可以改进吗？我确实这么认为：例如不使用默认编译标志-fwrapv（通过包括extra_compile_args = ["-fno-wrapv"]在设置中）将改善结果assembly并导致形状求和加速 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()]
)))
