Question

numpy对矢量操作优化的一个缺陷是它一次只能处理一个操作。这意味着，当我们对numpy矢量进行A * B + C这样的操作时，首先要等待A * B操作完成，数据保存在一个临时矢量中，然后将这个新的矢量和C相加。正如例6-14中使用就地操作的扩散代码所示。

然而，有许多模块可以对这点进行优化。numexpr模块可以将整个矢量表达式　 编译成非常高效的代码，可以将缓存失效以及临时变量的数量最小化。另外，它还能利用多个CPU核心（更多信息见第9章）以及Intel芯片专用的指令集来将速度最大化。

很容易修改代码来使用numexpr：我们只需要将表达式重写为使用本地变量的字符串即可。表达式会在后台被编译成优化过的代码（并被缓存来确保相同的表达式不会导致同样的编译过程发生多次）并运行。例6-19显示了evolve函数改用numexpr是多么的简单。我们在evaluate函数中使用out参数，这样numexpr就不会为返回一个新矢量而分配内存。

例6-19　使用numexpr来进一步优化大矩阵操作

from numexpr import evaluate

def evolve(grid, dt, next_grid, D=1):
  laplacian(grid, next_grid)
  evaluate("next_grid*D*dt+grid", out=next_grid)

numexpr的一个重要特点是它考虑到了CPU缓存。它特地移动数据来让各级　CPU缓存能够拥有正确的数据让缓存失效最小化。当我们对更新后的代码运行perf（例6-20）时，我们看到速度的确是提升了。但是，如果我们针对一个较小的512×512的矩阵（见本章最后的图6-4）比较速度，我们会看到大约15%的速度下降。这是为什么？

例6-20　使用了numpy就地内存操作特化laplacian函数以及numexpr的性能指标

$ perf stat -e cycles,stalled-cycles-frontend,stalled-cycles-backend,instructions,\
  cache-references,cache-misses,branches,branch-misses,task-clock,faults,\
  minor-faults,cs,migrations -r 3 python diffusion_numpy_memory2_numexpr.py

Performance counter stats for 'python diffusion_numpy_memory2_numexpr.py' (3 runs):

   5,940,414,581 cycles           #   1.447 GHz
   3,706,635,857 stalled-cycles-frontend  #  62.40% frontend cycles idle
   2,321,606,960 stalled-cycles-backend   #  39.08% backend cycles idle
   6,909,546,082 instructions       #  1.16 insns per cycle
                      #  0.54 stalled cycles per insn
   261,136,786 cache-references     #  63.628 M/sec
    11,623,783 cache-misses       #   4.451 % of all cache refs
   627,319,686 branches         # 152.851 M/sec
     8,443,876 branch-misses      #   1.35% of all branches
   4104.127507 task-clock         #   1.364 CPUs utilized
       9,786 page-faults        #   0.002 M/sec
       9,786 minor-faults       #   0.002 M/sec
       8,701 context-switches     #   0.002 M/sec
        60 CPU-migrations       #   0.015 K/sec

   3.009811418 seconds time elapsed

numexpr引入的大多数额外的机制都跟缓存相关。当我们的矩阵较小且计算所需的所有数据都能被放入缓存时，这些额外的机制只是白白增加了更多的指令而不能对性能有所帮助。另外，将字符串编译成矢量操作也会有很大的开销。当程序运行的整体时间较少时，这个开销就会变得相当引人注意。不过，当我们增加矩阵的大小时，我们会发现numexpr比原生numpy更好地利用了我们的缓存。而且，numexpr利用了多核来进行计算并尝试填满每个核心的缓存。当矩阵较小时，管理多核的开销盖过了任何可能的性能提升。

我们用来跑测试的电脑有20480 KB的缓存（Intel Xeon E5-2680）。因为有两个数组需要处理，一个作为输入，另一个作为输出，所以我们可以轻易计算出需要足以填满缓存的矩阵大小。矩阵元素的总数是20480 KB/64 bit = 2560000。因为我们有两个矩阵，这一总数被平分到两个对象中（所以每个对象最多可以有2560000 / 2 = 1280000个元素）。最后，对这个数字求平方根可以让我们知道能存放这么多元素的矩阵的大小。总的来说，这意味着大概两个大小为1131的2维数组就会填满缓存（ ）。但实际上，我们自己没办法完全填满缓存（其他程序会占用部分缓存），所以现实来说大概能填入两个800×800的数组。见表6-1和表6-2，我们可以看到当矩阵大小从512×512跳到1024×1024时，numexpr代码的性能就开始超越纯numpy。