python - 为什么循环在这里胜过索引？

Question

几年前，有人posted在 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))

第三个使用单个循环创建距离矩阵(考虑到只有 1,000 个二维点的距离矩阵有 100 万个条目，这显然是很多循环)。乍一看，这个函数就像我在学习 NumPy 时编写的代码，我会先编写 Python 代码，然后逐行翻译它来编写 NumPy 代码。

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

带有循环的函数实际上显着优于另外两个:

from numpy import mat, zeros, newaxis

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

线程中的一位参与者 (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 并且 calcDistanceMatrixFastEuclidean2 将是 fast2。

fast1 基于 Matlab 的处理方式，repmap 函数证明了这一点。 repmap 函数在这种情况下创建一个数组，它只是一遍又一遍重复的原始数据。但是，如果您查看该函数的代码，它的效率非常低。它使用许多 Numpy 函数(3 个 reshape 和 2 个 repeat)来执行此操作。 repeat 函数还用于创建一个数组，其中包含每个数据项重复多次的原始数据。如果我们的输入数据是 [1,2,3] 那么我们从中减去 [1,2,3,1,2,3,1,2,3] [1,1,1,2,2,2,3,3,3]。 Numpy 不得不在运行 Numpy 的 C 代码之间创建很多额外的矩阵，这本可以避免。

fast2 使用 Numpy 的更多繁重工作，而无需在 Numpy 调用之间创建尽可能多的矩阵。 fast2 遍历点的每个维度，进行减法运算并保留每个维度之间平方差的总和。只有在最后才完成平方根。到目前为止，这听起来可能不如 fast1 高效，但是 fast2 通过使用 Numpy 的索引避免了执行 repmat 操作。为了简单起见，让我们看一下一维情况。 fast2 创建一维数据数组，然后从二维 (N x 1) 数据数组中减去它。这会在每个点和所有其他点之间创建差异矩阵，而无需使用 repmat 和 repeat，从而绕过创建大量额外数组。在我看来，这就是真正的速度差异所在。 fast1 在矩阵之间创建了很多额外的东西(并且它们的计算成本很高)以找到点之间的差异，而 fast2 更好地利用 Numpy 的强大功能来避免这些差异。

顺便说一句，这里是 fast2 的一个稍微快一点的版本:

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 创建为零矩阵。