列方向求和与行方向求和:为什么在NumPy中看不到区别?
我测试了在这个讲座中展示的一个例子,使用的是numpy(第20页/共57页)。
文中提到,a[:,1].sum()
的执行时间是9.3毫秒,而a[1,:].sum()
只需要72微秒。
我尝试复现这个结果,但没有成功。我是测量错了吗?还是说自2010年以来,NumPy发生了变化?
$ python2 -m timeit -n1000 --setup \
'import numpy as np; a = np.random.randn(4000,4000);' 'a[:,1].sum()'
1000 loops, best of 3: 16.5 usec per loop
$ python2 -m timeit -n1000 --setup \
'import numpy as np; a = np.random.randn(4000,4000);' 'a[1,:].sum()'
1000 loops, best of 3: 13.8 usec per loop
$ python2 --version
Python 2.7.7
$ python2 -c 'import numpy; print numpy.version.version'
1.8.1
虽然我能测到第二种写法的好处(可能是因为numpy使用C风格的行优先存储,导致缓存命中率更高),但我没有看到像pytables的贡献者所说的那种明显差异。
而且,我似乎在使用列求和和行求和时,并没有看到更多的缓存未命中。
编辑
到目前为止,我意识到我使用
timeit
模块的方式不对。对同一个数组(或数组的行/列)进行多次运行几乎肯定会被缓存(我有32KiB
的L1数据缓存,所以一行数据很容易就能放进去:4000 * 4字节 = 15k < 32k
)。我使用了@alim在回答中的脚本,设置单次循环(
nloop=1
)和十次试验(nrep=10
),并改变随机数组的大小(n x n
),我测量到n row/us col/us penalty col 1k 90 100 1 4k 100 210 2 10k* 110 350 3.5 20k* 120 1200 10
*
n=10k
及以上的大小已经无法放入L1d缓存了。
我仍然不确定造成这种情况的原因,因为perf
显示快速的行求和的缓存未命中率大致相同(有时甚至更高)。
Perf
数据:
nloop = 2
和nrep=2
,所以我预计第二次运行时一些数据仍在缓存中。
行求和 n=10k
perf stat -B -e cache-references,cache-misses,L1-dcache-loads,L1-dcache-load-misses,L1-dcache-stores,L1-dcache-store-misses,L1-dcache-prefetches,cycles,instructions,branches,faults,migrations ./answer1.py 2>&1 | sed 's/^/ /g'
row sum: 103.593 us
Performance counter stats for './answer1.py':
25850670 cache-references [30.04%]
1321945 cache-misses # 5.114 % of all cache refs [20.04%]
5706371393 L1-dcache-loads [20.00%]
11733777 L1-dcache-load-misses # 0.21% of all L1-dcache hits [19.97%]
2401264190 L1-dcache-stores [20.04%]
131964213 L1-dcache-store-misses [20.03%]
2007640 L1-dcache-prefetches [20.04%]
21894150686 cycles [20.02%]
24582770606 instructions # 1.12 insns per cycle [30.06%]
3534308182 branches [30.01%]
3767 faults
6 migrations
7.331092823 seconds time elapsed
列求和 n=10k
perf stat -B -e cache-references,cache-misses,L1-dcache-loads,L1-dcache-load-misses,L1-dcache-stores,L1-dcache-store-misses,L1-dcache-prefetches,cycles,instructions,branches,faults,migrations ./answer1.py 2>&1 | sed 's/^/ /g'
column sum: 377.059 us
Performance counter stats for './answer1.py':
26673628 cache-references [30.02%]
1409989 cache-misses # 5.286 % of all cache refs [20.07%]
5676222625 L1-dcache-loads [20.06%]
11050999 L1-dcache-load-misses # 0.19% of all L1-dcache hits [19.99%]
2405281776 L1-dcache-stores [20.01%]
126425747 L1-dcache-store-misses [20.02%]
2128076 L1-dcache-prefetches [20.04%]
21876671763 cycles [20.00%]
24607897857 instructions # 1.12 insns per cycle [30.00%]
3536753654 branches [29.98%]
3763 faults
9 migrations
7.327833360 seconds time elapsed
编辑2
我想我理解了一些方面,但我觉得问题还没有得到解答。目前我认为这个求和的例子并没有揭示任何关于CPU缓存的内容。为了消除numpy/python带来的不确定性,我尝试在C中进行求和,并在下面的回答中给出了结果。
4 个回答
我在使用Numpy 1.9.0.def-ff7d5f9时,发现你发的两个测试行的执行速度差了10倍。我觉得你的电脑和你用来编译Numpy的工具,对速度的影响可能和Numpy的版本一样重要。
不过实际上,我觉得像这样只对一行或一列进行计算的情况并不常见。我认为更好的测试方式是比较对所有行进行计算
a.sum(axis=0)
和对所有列进行计算
a.sum(axis=1)
对我来说,这两种操作的速度差别很小(对列进行计算大约只需要对行计算的95%的时间)。
补充:一般来说,我对比较微秒级别的操作速度非常谨慎。在安装Numpy时,确保有一个好的BLAS库是非常重要的,因为在大多数大型矩阵操作(比如矩阵乘法)中,这个库负责处理大部分的计算工作。在比较BLAS库时,最好用像矩阵乘法这样的复杂操作来做比较,因为这才是你花费时间最多的地方。我发现有时候,性能较差的BLAS库在向量乘法上反而比性能更好的BLAS库快一点。但更糟糕的是,像矩阵乘法和特征值分解这样的操作要慢上很多倍,而这些操作比简单的向量乘法重要得多。我认为这些差异出现的原因是,写一个相对快速的向量乘法在C语言中并不难,但写一个好的矩阵乘法需要更多的思考和优化,而矩阵乘法的成本更高,所以好的BLAS库会把精力放在这里。
在Numpy中也是如此:任何优化都会集中在较大的操作上,而不是小操作上,所以不要过于纠结小操作之间的速度差异。此外,很难判断小操作的速度差异是否真的是计算时间造成的,还是因为为了优化更复杂的操作而增加的额外开销。
我写了一个求和的例子,使用的是C语言:结果是关于CPU时间的测量,我总是用gcc -O1 using-c.c
来编译(gcc版本:gcc version 4.9.0 20140604)。下面是源代码。
我选择的矩阵大小是n x n
。当n<2k
时,行和列的求和没有明显的差别(对于n=2k
,每次运行大约6-7微秒)。
行求和
n first/us converged/us
1k 5 4
4k 19 12
10k 35 31
20k 70 61
30k 130 90
例如n=20k
Run 0 taken 70 cycles. 0 ms 70 us
Run 1 taken 61 cycles. 0 ms 60 us # this is the minimum I've seen in all tests
Run 1 taken 61 cycles. 0 ms 61 us
<snip> (always 60/61 cycles)
列求和
n first/us converged/us
1k 5 4
4k 112 14
10k 228 32
20k 550 246
30k 1000 300
例如n=20k
Run 0 taken 552 cycles. 0 ms 552 us
Run 1 taken 358 cycles. 0 ms 358 us
Run 2 taken 291 cycles. 0 ms 291 us
Run 3 taken 264 cycles. 0 ms 264 us
Run 4 taken 252 cycles. 0 ms 252 us
Run 5 taken 275 cycles. 0 ms 275 us
Run 6 taken 262 cycles. 0 ms 262 us
Run 7 taken 249 cycles. 0 ms 249 us
Run 8 taken 249 cycles. 0 ms 249 us
Run 9 taken 246 cycles. 0 ms 246 us
讨论
行求和更快。它对缓存的利用不大,也就是说,重复的求和并没有比第一次求和快多少。列求和则慢得多,但在进行5-8次迭代时速度会逐渐提高。特别是在n=4k
到n=10k
之间,缓存的作用使得速度提高了大约十倍。在更大的数组中,速度提升大约只有两倍。我还观察到,行求和收敛得非常快(经过一两次试验后),而列求和的收敛则需要更多的迭代(5次或更多)。
我得到的主要教训是:
- 对于大数组(超过2k个元素),求和速度是有差别的。我认为这是因为从RAM到L1d缓存获取数据时的协同效应。虽然我不知道一次读取的块/行大小,但我假设它大于8字节。因此,下一个要加的元素已经在缓存中了。
- 列求和的速度首先受到内存带宽的限制。当从RAM读取分散的数据块时,CPU似乎会缺乏数据。
- 在重复进行求和时,期望某些数据不需要从RAM中获取,而是已经存在于L2/L1d缓存中。对于行求和,这在
n>30k
时才明显,而对于列求和,在n>2k
时就已经显现出来。
使用perf
工具,我没有看到很大的差别。不过,C程序的大部分工作是用随机数据填充数组。我不知道如何消除这些“准备”数据...
下面是这个例子的C代码:
#include <stdio.h>
#include <stdlib.h> // see `man random`
#include <time.h> // man time.h, info clock
int
main (void)
{
// seed
srandom(62);
//printf ("test %g\n", (double)random()/(double)RAND_MAX);
const size_t SIZE = 20E3;
const size_t RUNS = 10;
double (*b)[SIZE];
printf ("Array size: %dx%d, each %d bytes. slice = %f KiB\n", SIZE, SIZE,
sizeof(double), ((double)SIZE)*sizeof(double)/1024);
b = malloc(sizeof *b * SIZE);
//double a[SIZE][SIZE]; // too large!
int i,j;
for (i = 0; i< SIZE; i++) {
for (j = 0; j < SIZE; j++) {
b[i][j] = (double)random()/(double)RAND_MAX;
}
}
double sum = 0;
int run = 0;
clock_t start, diff;
int usec;
for (run = 0; run < RUNS; run++) {
start = clock();
for (i = 0; i<SIZE; i++) {
// column wise (slower?)
sum += b[i][1];
// row wise (faster?)
//sum += b[1][i];
}
diff = clock() - start;
usec = ((double) diff*1e6) / CLOCKS_PER_SEC; // https://stackoverflow.com/a/459704/543411
printf("Run %d taken %d cycles. %d ms %d us\n",run, diff, usec/1000, usec%1000);
}
printf("Sum: %g\n", sum);
return 0;
}
有趣。我可以复现Sebastian的表现:
In [21]: np.__version__
Out[21]: '1.8.1'
In [22]: a = np.random.randn(4000, 4000)
In [23]: %timeit a[:, 1].sum()
100000 loops, best of 3: 12.4 µs per loop
In [24]: %timeit a[1, :].sum()
100000 loops, best of 3: 10.6 µs per loop
不过,如果我试着用一个更大的数组:
In [25]: a = np.random.randn(10000, 10000)
In [26]: %timeit a[:, 1].sum()
10000 loops, best of 3: 21.8 µs per loop
In [27]: %timeit a[1, :].sum()
100000 loops, best of 3: 15.8 µs per loop
但是,如果我再试一次:
In [28]: a = np.random.randn(10000, 10000)
In [29]: %timeit a[:, 1].sum()
10000 loops, best of 3: 64.4 µs per loop
In [30]: %timeit a[1, :].sum()
100000 loops, best of 3: 15.9 µs per loop
所以,我不太确定这里发生了什么,但这种波动可能是由于缓存的影响。也许新的架构在预测数据访问模式上更聪明,因此在预取数据时表现得更好?
无论如何,为了比较,我使用的是NumPy 1.8.1,操作系统是Linux Ubuntu 14.04,电脑是配有i5-3380M CPU @ 2.90GHz的笔记本。
编辑:经过一番思考,我可以说,第一次timeit执行求和时,数据是从内存中取出来的,但第二次运行这个操作时,数据已经在缓存中了(无论是按行还是按列),所以执行得很快。由于timeit会取运行时间的最小值,这就是为什么我们没有看到时间上有很大差别的原因。
另一个问题是,为什么我们有时会看到这种差异(使用timeit)。不过缓存的行为很奇怪,尤其是在多核机器同时执行多个进程时。
我觉得你复制的结果没有什么问题,但要记住那些幻灯片是2010年的,numpy自那时以来变化很大。根据numpy发布的日期,我猜Francesc可能在用的是v1.5版本。
这里有一个用来比较行和列求和速度的脚本:
#!python
import numpy as np
import timeit
print "numpy version == " + str(np.__version__)
setup = "import numpy as np; a = np.random.randn(4000, 4000)"
rsum = "a[1, :].sum()"
csum = "a[:, 1].sum()"
nloop = 1000
nrep = 3
print "row sum:\t%.3f us" % (
min(timeit.repeat(rsum, setup, repeat=nrep, number=nloop)) / nloop * 1E6)
print "column sum:\t%.3f us" % (
min(timeit.repeat(csum, setup, repeat=nrep, number=nloop)) / nloop * 1E6)
我发现使用numpy v1.5时,列求和的速度大约慢了50%:
$ python sum_benchmark.py
numpy version == 1.5.0
row sum: 8.472 us
column sum: 12.759 us
而你用的v1.8.1版本,速度大约慢了30%:
$ python sum_benchmark.py
numpy version == 1.8.1
row sum: 12.108 us
column sum: 15.768 us
有趣的是,最近的numpy版本中,这两种求和方式的速度其实都有点变慢。要想搞清楚具体原因,我得深入研究numpy的源代码。
更新
- 顺便说一下,我是在一台运行Ubuntu 14.04(内核v3.13.0-30)的四核i7-2630QM CPU @ 2.0GHz的电脑上进行测试的。两个版本的numpy都是通过pip安装的,并使用GCC-4.8.1编译的。
- 我意识到我最初的基准测试脚本不是特别清晰——你需要把总时间除以循环次数(1000),才能得到每次调用的时间。
- 而且取重复测试中的最小值可能更合理,而不是平均值,因为这样更能代表执行时间的下限(再加上背景进程等带来的波动)。
我已经相应地更新了我的脚本和结果
我们还可以通过为每次调用创建一个全新的随机数组来消除缓存的影响(时间局部性)——只需将nloop
设置为1,nrep
设置为一个比较小的数字(除非你真的喜欢看油漆干),比如10。
nloop=1
,nreps=10
在一个4000x4000的数组上:
numpy version == 1.5.0
row sum: 47.922 us
column sum: 103.235 us
numpy version == 1.8.1
row sum: 66.996 us
column sum: 125.885 us
这样就更接近了,但我仍然无法完全复制Francesc幻灯片上显示的巨大效果。不过,这也许并不奇怪——这个效果可能非常依赖于编译器、架构和/或内核。