为什么循环优于索引？

Question

几年前，有人在Active State Recipes上使用了三个python/NumPy函数，每个函数接受相同的参数并返回相同的结果，即距离矩阵。在

其中两个是从已发布的源代码中获取的；它们都是——或者在我看来是——惯用的numpy代码。创建距离矩阵所需的重复计算由numpy优雅的索引语法驱动。以下是其中之一：

from numpy.matlib import repmat, repeat

def calcDistanceMatrixFastEuclidean(points):
  numPoints = len(points)
  distMat = sqrt(sum((repmat(points, numPoints, 1) - 
             repeat(points, numPoints, axis=0))**2, axis=1))
  return distMat.reshape((numPoints,numPoints))

第三个使用一个循环创建了距离矩阵（显然，这是一个很大的循环，因为距离矩阵只有1000个二维点，有一百万个条目）。乍一看，这个函数看起来像我学习NumPy时编写的代码，我会先编写Python代码，然后逐行翻译，来编写NumPy代码。在

在activestate发布几个月后，比较这三者的性能测试结果发布在NumPy邮件列表上的thread中并进行了讨论。在

实际上，带有循环的函数的性能明显优于其他两个函数：

^{pr2}$

线程中的一位参与者（Keir Mierle）给出了一个可能是真的原因：

The reason that I suspect this will be faster is that it has better locality, completely finishing a computation on a relatively small working set before moving onto the next one. The one liners have to pull the potentially large MxN array into the processor repeatedly.

据这名发帖人自己说，他的话只是一种怀疑，似乎没有进一步讨论过。在

关于如何解释这些结果还有什么想法吗？在

特别是，是否有一个有用的规则——关于何时循环和何时索引——可以从这个例子中提取出来作为编写numpy代码的指导？在

对于那些不熟悉NumPy的人，或者没有看过代码的人来说，这种比较并不是基于边缘情况——如果是的话，我肯定不会那么感兴趣。相反，这种比较涉及到一个函数，该函数在矩阵计算中执行一个公共任务（即，创建一个给定两个前件的结果数组）；此外，每个函数依次由最常见的numpy内置函数组成。在

Answer 1

TL；DR上面的第二个代码只在点的维数上循环（对于3D点，循环次数是for循环的3倍），因此循环次数不多。上面第二个代码的真正加速是它更好地利用Numpy的能力，以避免在找到点之间的差异时创建一些额外的矩阵。这减少了内存使用和计算工作量。在

更长的解释 我认为calcDistanceMatrixFastEuclidean2函数的循环可能在欺骗你。它只是在点的维数上循环。对于1D点，循环只执行一次，对于2D，执行两次，对于3D，执行三次。这真的不是什么循环。在

让我们分析一下代码，看看为什么其中一个比另一个快。calcDistanceMatrixFastEuclidean我将调用fast1，而{}将是{}。在

fast1是基于Matlab的做事方式的，repmap函数就证明了这一点。在本例中，repmap函数创建了一个数组，它只是反复重复的原始数据。但是，如果查看函数的代码，则效率非常低。它使用许多Numpy函数（3reshapes和2repeats）来实现这一点。repeat函数还用于创建一个数组，该数组包含每个数据项重复多次的原始数据。如果我们的输入数据是[1,2,3]，那么我们将从[1,1,1,2,2,2,3,3,3]中减去{}。Numpy不得不在运行Numpy的C代码之间创建大量额外的矩阵，这是可以避免的。在

fast2使用了更多Numpy的重担，而没有在Numpy调用之间创建那么多的矩阵。fast2遍历点的每个维度，进行减法，并保持每个维度之间的平方差的累计。只有到最后才算平方根。到目前为止，这听起来不像fast1那么有效，但是fast2通过使用Numpy的索引避免了repmat的工作。为了简单起见，让我们看一下1D案例。fast2生成数据的1D数组，并从2D（nx1）数组中减去它。这将在每个点和所有其他点之间创建差分矩阵，而不必使用repmat和repeat，从而避免了创建大量额外的数组。这就是我认为真正的速度差异所在。fast1在矩阵之间创建了很多额外的元素（它们的创建是非常昂贵的计算）来发现点之间的差异，而{}则更好地利用Numpy的能力来避免这些差异。在

{cd5>稍微快一点

def calcDistanceMatrixFastEuclidean3(nDimPoints):
  nDimPoints = array(nDimPoints)
  n,m = nDimPoints.shape
  data = nDimPoints[:,0]
  delta = (data - data[:,newaxis])**2
  for d in xrange(1,m):
    data = nDimPoints[:,d]
    delta += (data - data[:,newaxis])**2
  return sqrt(delta)

不同的是，我们不再创建delta作为零矩阵。在

Answer 2

dis为了好玩：

dis.dis(calcDistanceMatrixFastEuclidean)

  2           0 LOAD_GLOBAL              0 (len)
              3 LOAD_FAST                0 (points)
              6 CALL_FUNCTION            1
              9 STORE_FAST               1 (numPoints)

  3          12 LOAD_GLOBAL              1 (sqrt)
             15 LOAD_GLOBAL              2 (sum)
             18 LOAD_GLOBAL              3 (repmat)
             21 LOAD_FAST                0 (points)
             24 LOAD_FAST                1 (numPoints)
             27 LOAD_CONST               1 (1)
             30 CALL_FUNCTION            3

  4          33 LOAD_GLOBAL              4 (repeat)
             36 LOAD_FAST                0 (points)
             39 LOAD_FAST                1 (numPoints)
             42 LOAD_CONST               2 ('axis')
             45 LOAD_CONST               3 (0)
             48 CALL_FUNCTION          258
             51 BINARY_SUBTRACT
             52 LOAD_CONST               4 (2)
             55 BINARY_POWER
             56 LOAD_CONST               2 ('axis')
             59 LOAD_CONST               1 (1)
             62 CALL_FUNCTION          257
             65 CALL_FUNCTION            1
             68 STORE_FAST               2 (distMat)

  5          71 LOAD_FAST                2 (distMat)
             74 LOAD_ATTR                5 (reshape)
             77 LOAD_FAST                1 (numPoints)
             80 LOAD_FAST                1 (numPoints)
             83 BUILD_TUPLE              2
             86 CALL_FUNCTION            1
             89 RETURN_VALUE

dis.dis(calcDistanceMatrixFastEuclidean2)

^{pr2}$

我不是dis方面的专家，但似乎您需要更多地研究第一个调用的函数，以了解它们为什么要花费一段时间。Python还有一个性能分析器工具，cProfile。在