在 Python 中构建矩阵的快速方法

Question

我一直在浏览这些问题，并且可以找到一些帮助，但我更喜欢通过直接询问来确认。所以这是我的问题。

我有一个维度为 N 的 (numpy) 数组 u，我想从中构建一个维度为 N^2 的方阵 k。基本上，每个矩阵元素k(i,j) 定义为k(i,j)=exp(-|u_i-u_j|^2)。

我的第一个天真方法是这样的，我相信这类似于 Fortran：

for i in range(N):
  for j in range(N):
    k[i][j]=np.exp(np.sum(-(u[i]-u[j])**2))

但是，这非常慢。例如，对于 N=1000，大约需要 15 秒。 我的其他方法如下（受其他问题/答案的启发）：

i, j = np.ogrid[:N,:N]
k = np.exp(np.sum(-(u[i]-u[j])**2,axis=2))

这样更快，对于 N=1000，结果几乎是瞬时的。 所以我有两个问题。

1) 为什么第一种方法这么慢，为什么第二种方法那么快？

2) 有更快的方法吗？对于 N=10000，它已经开始需要相当长的时间了，所以我真的不知道这是否是“正确”的做法。

提前谢谢你！

P.S：矩阵是对称的，所以必须有一种方法可以通过只计算矩阵的上半部分来使过程更快，但我的问题更多与操作数组的方式等有关。

Answer 1

首先，一个小说明，如果u可以改写为u = np.arange(N)，则不需要使用np.sum。这似乎是这种情况，因为您写道它的维度为N。

1) 第一个问题： 在 Python 中访问索引很慢，所以如果有办法不使用它，最好不要使用 []。另外，您可以多次调用np.exp 和np.sum，而它们可以调用向量和矩阵。因此，您的第二个建议更好，因为您一次性计算了 k，而不是逐个元素地计算。

2) 第二个问题： 就在这里。您应该考虑仅使用 numpy 函数而不使用索引（大约快 3 倍）：

k = np.exp(-np.power(np.subtract.outer(u,u),2))

（注意：您可以保留**2 而不是np.power，这会更快但精度更小）

---

编辑（考虑到u 是一个元组数组）

使用元组数据，就有点复杂了：

ma = np.subtract.outer(u[:,0],u[:,0])**2
mb = np.subtract.outer(u[:,1],u[:,1])**2
k = np.exp(-np.add(ma, mb))

你必须使用两次np.substract.outer，因为如果你一次性使用它会返回一个 4 维数组（并计算大量无用数据），而 u[i]-u[j] 返回一个 3 维数组。

我使用np.add 而不是np.sum，因为它保留了数组维度。

注意：我检查过

N = 10000
u = np.random.random_sample((N,2))

我返回的结果与您的建议相同。 （但快 1.7 倍）