Question

我们在例6-6写的Python代码还有一个问题，也是使用Python进行此类矢量计算的核心问题：原生Python并不支持矢量操作。这有两个原因：Python列表存储的是指向实际数据的指针，且Python字节码并没有针对矢量操作进行优化，所以for循环无法预测何时使用矢量操作能带来好处。

Python列表储存的是指针意味着列表保存的并不是实际的数据，而是这些数据的位置。在大多数情况下，这是一种优势，因为它允许我们在列表中存储任意类型的数据。但是对于矢量和矩阵操作，这会导致性能的大幅下降。

性能下降的原因在于每次我们需要获得grid矩阵中的元素时，我们都必须进行多次查找。比如，grid[5][2]这样的语句需要我们首先对列表grid查找索引5。这会返回一个指针指向数据所在的位置。然后我们需要第二次列表查找找到索引2的元素。查到之后，我们得到的是实际数据所储存的位置。

一次这样的查找开销并不大，在大多数情况下可以忽略。但是，如果我们想要定位的数据在内存的一个连续存储块内，我们本来可以一次性将所有的数据全部读取出来，而不是分别对每个元素进行两次操作。这是数据分片的一个主要问题：当你的数据被分成小片，你只能对每一片分别进行传输，而不是一次性传输整个块。这意味着你引入了更多的内存传输开销，且强制CPU在数据传输的过程中等待。我们可以用perf看到在缓存失效的情况下这个问题会有多严重。

这个在正确的时候将正确的数据传输给CPU的问题被称为“冯诺伊曼瓶颈”。意思是现代计算机所使用的层次化的内存架构会导致CPU和内存之间的带宽受到限制。如果我们数据传输的速度可以无限快，我们就不需要任何缓存，因为CPU可以立即获得任何它需要的数据。此时瓶颈就不再存在。

由于数据传输的速度不可能无限快，我们必须从RAM中预取数据并将其保存在一个更小但更快的CPU缓存中，并希望当CPU需要某个数据时，它可以从中更快读取到。虽然这已经是一个严重理想化了的场景，我们依然可以看到其中的一些问题——我们如何知道未来需要哪些数据？CPU内部的分支预测和流水线技术会试图在处理当前指令的同时预测其下一条指令并将相应的内存读进缓存。但是减少瓶颈最好的方法是让代码知道如何分配我们的内存以及如何使用我们的数据进行计算。

探测内存移动至CPU的性能相当困难，不过，Linux上的perf工具可以让我们洞察CPU如何处理运行中的程序。比如，我们可以对例6-6的纯Python代码运行perf并看到CPU运行我们代码的效率。结果见例6-8。注意该例以及之后的perf例子的输出都被截取以适应页面边界。被删除的数据包括各测量值的方差，用于表示测量值在几次测量中发生了多大的变化。这有助于看到测量值在多大程度上依赖于实际的程序特性以及多大程度上受到来自系统的其他干扰，比如其他正在使用系统资源的程序。

例6-8　减少内存分配后纯Python 2阶扩散函数的性能指标

$ 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_python_memory.py

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

 329,155,359,015 cycles           #  3.477 GHz
  76,800,457,550 stalled-cycles-frontend  #   23.33% frontend cycles idle
  46,556,100,820 stalled-cycles-backend   #   14.14% backend cycles idle
 598,135,111,009 instructions       #  1.82 insns per cycle
                      #  0.13 stalled cycles per insn
   35,497,196 cache-references      #  0.375 M/sec
   10,716,972 cache-misses        #   30.191 % of all cache refs
133,881,241,254 branches          # 1414.067 M/sec
  2,891,093,384 branch-misses       #  2.16% of all branches
   94678.127621 task-clock        #  0.999 CPUs utilized
      5,439 page-faults         #  0.057 K/sec
      5,439 minor-faults        #  0.057 K/sec
      125 context-switches      #  0.001 K/sec
        6 CPU-migrations      #  0.000 K/sec

   94.749389121 seconds time elapsed

6.3.1　理解perf

让我们花一秒钟来理解perf告诉我们的各种性能指标以及它们跟我们代码的关系。task-clock指标告诉我们的任务花了多少个时钟周期。这跟总体的运行时间不同，因为如果我们的程序花了一秒钟来运行但是使用了两个CPU，那么task-clock将是1000（task-clock的单位是毫秒）。方便的是，perf会帮我们计算并在该指标旁边告诉我们有多少个CPU被使用了。这个数字不完全等于1是因为进程有一段时间依赖于其他子系统的指令（比如分配内存时）。

context-switches和CPU-migrations告诉我们程序在等待内核操作（如I/O操作）完成时，为了让其他进程得以运行，被挂起或迁移到另一个CPU核心上执行的次数。当一个context-switch发生时程序的执行会被挂起，让另一个程序得以执行。这是一个对时间要求非常精细的任务，也是我们需要尽量避免的，但是我们对它的发生无能为力。只要进程允许切换，内核就会接手；不过，我们可以做一些事来抑制内核切换我们的程序。总的来说，内核会在程序进行I/O（比如读取内存、磁盘或网络）时将其挂起。我们在后续章节会看到，我们可以用异步操作来确保我们的程序在等待I/O时继续使用CPU，这会让我们的进程继续运行而不被切换出去。另外，我们还可以设置程序的nice值来给我们的程序更高的优先级以防止内核将它切换出去。类似的，CPU-migrations会发生在进程被挂起并迁移到另一个CPU上继续执行的情况，这是为了让所有的CPU都有同样程度的利用率。这可以被认为是一个特别糟糕的进程切换，因为我们的程序不仅被暂时挂起，而且丢失了L1缓存内所有的数据（每个CPU都有它自己的L1缓存）。

page-fault是现代UNIX内存分配机制的一部分。分配内存时，内核除了告诉程序一个内存的引用地址以外没做任何事。但是，之后在这块内存第一次被使用时，操作系统会抛出一个缺页小中断，这将暂停程序的运行并正确分配内存。这被称为延迟分配系统。虽然这种手段相比以前的内存分配系统是一个很大的优化，缺页小中断本身依然是一个相当昂贵的操作，因为大多数操作都发生在你的程序外部。另外还有一种缺页大中断，发生于当你的程序需要从设备（磁盘、网络等）上请求还未被读取的数据时。这些操作更加昂贵，因为他们不仅中断了你的程序，还需要读取数据所在的设备。这种缺页不总是影响CPU密集的工作，但是，它会给任何需要读写磁盘或网络的程序带来痛苦。

一旦我们引用内存中的数据，数据就需要通过内存的多个层级（对此的讨论见1.1.3节“通信层”）。每当我们引用缓存中的数据，cache-references指标就会增加。如果我们的缓存中还没有数据并需要从RAM获取，cache-miss指标的计数就会增加。如果我们读取一个最近刚读过的数据（那个数据还在缓存里）或者读取的数据在最近刚读过的数据附近（缓存从RAM中按块读取数据），缓存失效就不会发生。缓存失效会拖慢CPU密集型工作，因为我们不仅需要等待数据从RAM被读取，而且中断了我们的执行流水线（马上会说到这个）。本章后面将讨论如何通过优化内存中的数据布局降低缓存失效。

Instructions指标告诉我们的代码让CPU总共执行了多少条指令。流水线技术让CPU能在同一时间执行多条指令，insns per cycle指标会告诉我们一个时钟周期内具体了几条。为了更好地优化流水线，stalled-cycles-frontend和stalled-cycles-backend告诉我们的程序在等待流水线的前端或后端填满指令时等待了多少时钟周期。这可能是由于一次缓存失效，一次错误的分支预测或一次资源冲突。流水线前端负责从内存中获取下一条指令并解码成一个合法的操作，而后端则负责实际执行操作。有了流水线，CPU就能在运行当前指令的同时获取并准备下一条指令。

branch是代码执行流程发生改变的地方。想象一条if..then语句——基于条件语句的结果，我们会执行不同的代码段。这就是代码执行的分支——下一条指令将是两者之一。为了优化，特别是在使用了流水线的情况下，CPU会试图猜测分支的走向并预取与其相关的指令。当预测出错时，就会发生stalled-cycles和branch-miss。分支失效是相当令人困惑的，并可能导致很多奇怪的现象（比如，某些循环在排序列表上跑得远比乱序列表快，仅仅是因为发生了更少的分支失效）。

　备忘　

如果你想要一个在CPU层面各种性能指标更彻底的解释，请参考Gurpur M. Prabhu的“计算机架构导论”。它讲述在很底层发生的故事，让你对代码运行时底层发生的一些事有一个很好的理解。

6.3.2　根据perf输出做出抉择

例6-8的性能指标告诉我们在运行我们的代码时，CPU需要访问L1/L2缓存35 497 196次，其中10 716 972次（30.191%）是由于请求的数据不在内存中而需要被获取。另外，我们可以看到每个CPU周期我们能平均执行1.82条指令，这告诉了我们通过流水线，乱序执行以及超线程技术（或者任何其他能够让你的CPU在一个时钟周期内执行多条指令的技术）给予我们的性能增幅。

数据的分片不仅增加了从内存传输数据到CPU的次数，还让你无法进行矢量计算，因为在计算时你的CPU缓存中并没有多个数据。我们在1.1.3节中解释过，矢量计算（或者说让CPU在同一时间进行多个计算）仅能发生在我们能够将相关数据填满CPU缓存的情况下。由于总线只能移动连续的内存数据，也就是说只有当格子中的数据在RAM中顺序储存时才有可能。由于我们的列表存储的是数据的指针而不是实际数据，格子中实际的值被分散在内存里，无法一次性全部复制。

我们可以通过使用array模块代替列表来减轻这个问题。array对象可以在内存中顺序储存数据，一段array实际代表了一段连续的内存。但是这并没有完全解决问题——现在我们有了内存中的连续存储的数据，但Python依然不知道如何对我们的循环进行矢量计算。我们想要的是让原本一次只能处理一个数组元素的循环现在一次能够处理多个数据，但是之前说过，Python没有这样的字节码优化（部分是由于Python语言极端的动态特性）。

　备忘　

为什么我们顺序存储在内存中的数据不能自动矢量化？如果我们看看CPU运行的机器指令，就会发现矢量操作（比如两个数组相乘）和非矢量操作使用的是不同的CPU计算单元和指令集。为了让Python能够使用这些特殊指令，我们必须有一个模块专门用来使用这些指令集。我们马上会看到numpy如何让我们能够访问这些特殊指令。

另外，由于实现细节，使用array类型来创建数据列表实际上比使用list要慢。这是因为array对象存储的数据是一个非常底层的抽象，在用户使用它们时必须被转换成一个Python兼容的版本。这个额外的开销在你每次索引一个array类型时都会发生。这一实现细节导致了array对象不适用于数学计算而更适用于在内存中存放类型固定的数据。

6.3.3　使用numpy

为了解决perf发现的分片问题，我们必须找到一个可以进行高效矢量操作的模块。幸运的是，numpy拥有我们需要的所有特性——它能将数据连续存储在内存中并支持数据的矢量操作。结果就是，任何我们对numpy数组的数学操作都能自动矢量化而无须我们显式遍历每一个元素。这不仅仅让矩阵计算更简单，而且也更快。让我们看这样一个例子：

from array import array import numpy def norm_square_list(vector): """ >>> vector = range(1000000) >>> %timeit norm_square_list(vector_list) 1000 loops, best of 3: 1.16 ms per loop """ norm = 0 for v in vector: norm += v*v return norm def norm_square_list_comprehension(vector): """ >>> vector = range(1000000) >>> %timeit norm_square_list_comprehension(vector_list) 1000 loops, best of 3: 913 μs per loop """ return sum([v*v for v in vector]) def norm_squared_generator_comprehension(vector): """ >>> vector = range(1000000) >>> %timeit norm_square_generator_comprehension(vector_list) 1000 loops, best of 3: 747 μs per loop """ return sum(v*v for v in vector) def norm_square_array(vector): """ >>> vector = array('l', range(1000000)) >>> %timeit norm_square_array(vector_array) 1000 loops, best of 3: 1.44 ms per loop """ norm = 0 for v in vector: norm += v*v return norm def norm_square_numpy(vector): """ >>> vector = numpy.arange(1000000) >>> %timeit norm_square_numpy(vector_numpy) 10000 loops, best of 3: 30.9 μs per loop """ return numpy.sum(vector * vector) # ❶ def norm_square_numpy_dot(vector): """ >>> vector = numpy.arange(1000000) >>> %timeit norm_square_numpy_dot(vector_numpy) 10000 loops, best of 3: 21.8 μs per loop """ return numpy.dot(vector, vector) # ❷

❶ 这创建了vector上的两个隐式循环，一个用于做乘法，另一个求和。这两个循环类似norm_squre_list_comprehension中的循环，只是使用了numpy优化的数学代码。

❷ 这才是numpy中求矢量范数的推荐做法，通过矢量化的numpy.dot操作。前面提供的低效的norm_squre_numpy代码仅作说明用途。

这段简单的numpy代码运行速度是norm_square_list的37.54倍，比使用Python列表“优化”后的版本也要快29.5倍。纯Python循环和列表优化版本在速度上的差距体现了让后台进行更多计算相比于Python代码显式计算的优势。 Python是用C语言写的，利用Python内建机制进行运算，我们就获得了原生C代码的速度。我们在numpy代码中有如此大的速度提升也是基于同样的原因：我们使用了高度优化且特殊构建的对象取代了通用的列表结构来处理数组。

除了更轻量级且特化的机制，numpy对象还给了我们内存的本地性和矢量操作，它们在处理数值计算时尤其重要。CPU是超快的，大多数时候，仅仅提升其获取所需的数据的速度就是最好的优化手段。我们之前使用perf工具运行每个函数显示了使用array的版本和纯Python版本的函数花费了大约1×10¹¹条指令，而numpy版本花费了大约3×10⁹条。另外，array和纯Python版本大约有80%缓存失效而numpy只有大约55%。

我们的norm_square_numpy代码中，在运算vector * vector时，numpy会为我们隐式调用一个循环。该循环和我们在其他例子中显式调用的一样：遍历vector内的每一项，乘以其自身。但是，由于我们让numpy来干这件事而不是自己用Python代码实现，numpy就可以进行任何它想要的优化。numpy在后台有极其优化的C代码来专门使用CPU的矢量操作。另外，numpy数组在内存中是连续储存的底层数字类型，和（array模块中的）array对象有同样的空间需求。

另外还有一个额外好处是，numpy支持我们将问题重新描述为一个点积。这让我们可以用一个单一操作来计算我们想要的值，而不是计算两个矢量的乘积并求和。在图6-3中可以看到，由于该函数的特化，norm_numpy_dot的效率大大超越了其他所有函数，这也是因为我们不需要存储vector * vector的中间结果。

图6-3　对不同长度的矢量求范数平方的函数运行时间