列和与行和：为什么我看不到使用NumPy的区别？

Question

我用numpy测试了这个talk[pytables]中演示的一个例子（第20/57页）。在

结果表明，a[:,1].sum()只需9.3ms，而{}只需72 us。在

我试图复制它，但没有成功。我量错了吗？或者从2010年开始，纽比的情况发生了变化？在

$ 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贡献者所说的那种巨大差异。在

另外，在使用列V行求和时，似乎看不到更多的缓存未命中。在

---

编辑

• 到目前为止，我的洞察力是我用错了timeit模块。使用同一个数组（或数组的行/列）重复运行几乎肯定会被缓存（我有一级数据缓存的32KiB，因此其中有一行很适合：4000 * 4 byte = 15k < 32k）。

• 使用@alim的answer中的脚本和一个单循环（nloop=1）和十次试验nrep=10，并改变我测量的随机数组（n x n）的大小

  ^{pr2}$

  *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中使用perf进行求和，结果如下所示。在

Answer 1

有趣。我可以重现塞巴斯蒂安的表演：

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

但是，如果我尝试使用更大的阵列：

^{pr2}$

但是，如果我再试一次：

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

所以，不确定这是怎么回事，但这种抖动可能是由于缓存效应。也许新的体系结构在预测模式访问方面更明智了，因此可以更好地进行预取？在

无论如何，为了便于比较，我使用的是Numpy1.8.1、LinuxUbuntu14.04和一台i5-3380mCPU@2.90GHz的笔记本电脑。在

编辑：在考虑过这一点之后，是的，我会说第一次执行sum时，列（或行）是从RAM中获取的，但是第二次操作运行时，数据在缓存中（对于行和列两种版本），所以它执行的速度很快。由于时间花费最少的运行，这就是为什么我们看不到时间上的巨大差异。在

另一个问题是为什么我们有时会看到差异（使用timeit）。但是缓存是一种奇怪的野兽，尤其是在一次执行多个进程的多核机器中。在

Answer 2

我看不出你的复制尝试有什么错，但请记住，这些幻灯片都是2010年的，从那时起numpy已经发生了很大的变化。根据dates of numpy releases，我猜弗朗西斯可能在用v1.5。在

使用此脚本对第v行列和进行基准测试：

#!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%：

^{pr2}$

相比之下，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的源代码 去理解为什么会这样。在

更新

• 我在2.0GHz的四核i7-2630qmcpu上运行ubuntu14.04（kernelv3.13.0-30）。两个版本的numpy都是pip安装的，并使用GCC-4.8.1进行编译。在
• 我意识到我最初的基准测试脚本并不是完全自解释的-您需要将总时间除以循环数（1000）才能得到每次调用的时间。在
• 它也是probably makes more sense to take the minimum across repeats rather than the average，因为这更可能代表执行时间的下限（在这个下限上，您将得到由于后台进程等原因而产生的变化）。在

我已经相应地更新了我的脚本和上面的结果

我们还可以通过为每个调用创建一个全新的随机数组来消除跨调用缓存的任何效果（临时局部性）——只需将nloop设置为1，将nrep设置为一个相当小的数字（除非您真的很喜欢看油漆干燥），比如10。在

在4000x4000阵列上的nloop=1，nreps=10：

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

这有点像，但我仍然无法真正复制弗朗西斯的幻灯片所显示的巨大效果。也许这并不令人惊讶，但是-效果可能非常依赖于编译器、体系结构和/或内核。在

Answer 3

我在C中编写了求和示例：结果显示为CPU time度量，我总是使用gcc -O1 using-c.c来编译（gcc版本：gcc版本4.9.0 20140604）。源代码如下。在

我选择矩阵大小为n x n。对于n<2k，行和列的总和没有任何可测量的差异（对于n=2k，每次运行6-7 us）。在

行总和

n     first/us      converged/us
 1k       5                 4
 4k      19                12
10k      35                31
20k      70                61
30k     130                90

e、 gn=20k ^{pr2}$

列

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到{}之间，这种增长最为明显，其中缓存有助于将速度提高约10倍。在较大的阵列中，加速仅为因子2。我还观察到，虽然行和收敛非常快（经过一次或两次试验），列求和收敛需要更多的迭代（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;
}