简单向量运算的优化（python）答案

【问题标题】：Optimisation of a simple vector operation (python)简单向量运算的优化（python）
【发布时间】：2014-10-17 17:17:28
【问题描述】：

在 python 中是否有一种快速的方法来执行一个简单的操作，从而产生一个矩阵，例如A[i,j] = a[i] - b[j] 给定两个数组 a 和 b（长度相同，但这可能不相关）？

更准确地说，我拥有的是二维空间中的 N 个点，其位置存储在两个数组 dx 和 dy 中，还有 N 个点的位置存储在 tx 和 ty 中。我需要一个矩阵

A[i,j] = (dx[j]-tx[i])**2+(dy[j]-ty[i])**2

我想到的唯一方法就是做

A = np.empty([nData,nData])
for i in range(nData):
        A[i] = (dx-tx[i])**2+(dy-ty[i])**2
return A

问题是这太慢了（nData 会很大）。如果可以更快地更改表示法，我们欢迎任何更改。

（顺便问一下，x**2 是不是比 x*x 慢？）

【问题讨论】：

请展示一个最小的数据示例和预期结果。

标签： python optimization numpy

【解决方案1】：

您想要计算点之间的所有成对平方欧几里得距离。最快的是使用scipy.distance.cdist:

>>> import numpy as np
>>> from scipy.spatial.distance import cdist
>>> x = np.random.rand(10, 2)
>>> t = np.random.rand(8, 2)

>>> cdist(x, t, 'sqeuclidean')
array([[ 0.61048982,  0.04379578,  0.30763149],
       [ 0.02709455,  0.30235292,  0.25135934],
       [ 0.21249888,  0.14024951,  0.28441688],
       [ 0.39221412,  0.01994213,  0.17699239]])

如果你想自己在 numpy.这样的事情应该可以解决问题：

>>> np.sum((x[:, None] - t)**2, axis=-1)
array([[ 0.61048982,  0.04379578,  0.30763149],
       [ 0.02709455,  0.30235292,  0.25135934],
       [ 0.21249888,  0.14024951,  0.28441688],
       [ 0.39221412,  0.01994213,  0.17699239]])

或者，为 x 和 y 坐标使用单独的数组：

>>> dx, dy = x.T
>>> tx, ty = t.T

>>> (dx[:, None] - tx)**2 + (dy[:, None] - ty)**2
array([[ 0.61048982,  0.04379578,  0.30763149],
       [ 0.02709455,  0.30235292,  0.25135934],
       [ 0.21249888,  0.14024951,  0.28441688],
       [ 0.39221412,  0.01994213,  0.17699239]])

【讨论】：

谢谢！像魅力一样工作。

【解决方案2】：

试试

>>> a = arange(1, 10)
>>> b = arange(1, 10)
>>> a.reshape(9, 1) - b.reshape(1, 9)
array([[ 0, -1, -2, -3, -4, -5, -6, -7, -8],
       [ 1,  0, -1, -2, -3, -4, -5, -6, -7],
       [ 2,  1,  0, -1, -2, -3, -4, -5, -6],
       [ 3,  2,  1,  0, -1, -2, -3, -4, -5],
       [ 4,  3,  2,  1,  0, -1, -2, -3, -4],
       [ 5,  4,  3,  2,  1,  0, -1, -2, -3],
       [ 6,  5,  4,  3,  2,  1,  0, -1, -2],
       [ 7,  6,  5,  4,  3,  2,  1,  0, -1],
       [ 8,  7,  6,  5,  4,  3,  2,  1,  0]])

该片段中发生的事情称为broadcasting，请查看该页面以获取解释。如果您不惜一切代价避免显式循环，Numpy 通常是最快的。谷歌搜索“numpy 矢量化”应该会为您提供详细信息。

翻译成你的例子，完整的公式是

(dx.reshape(len(dx), 1) - tx.reshape(1, len(tx)))**2 + \
(dy.reshape(len(dy), 1) - ty.reshape(1, len(ty)))**2

【讨论】：