Python中相似度矩阵的高效计算（NumPy）

Question

设X 为Bxn numpy 矩阵，即，

import numpy as np
B = 10
n = 2
X = np.random.random((B, n))

现在，我有兴趣计算所谓的核（甚至相似度）矩阵K，其形状为BxB，其{i,j}-th 元素如下所示：

K(i,j) = fun(x_i, x_j)

其中x_t 表示矩阵X 的第t 行，fun 是x_i、x_j 的某个函数。例如，这个函数可以是所谓的 RBF 函数，即

K(i,j) = exp(-|x_i - x_j|^2)。

为此，一种天真的方法如下：

K = np.zeros((B, B))
for i in range(X.shape[0]):
    x_i = X[i, :]
    for j in range(X.shape[0]):
        x_j = X[j, :]
        K[i, j] = np.exp(-np.linalg.norm(x_i - x_j, 2) ** 2)

为了提高效率，我想要以矢量化的方式进行上述操作。你能帮忙吗？

Answer 1

如果您利用 broadcasting 的力量，这在 numpy 中当然是可能的。

您只需以矢量化的方式编写内部距离范数计算：

X1 = X[:, np.newaxis, :]
X2 = X[np.newaxis, :, :]
K = np.exp(-np.sum((X1 - X2)**2, axis=-1))

Answer 2

不要向量化，直接编译即可

这几乎每次都更快，并且代码更易于阅读。 由于可以使用像 Numba 这样的良好 jit 编译器，因此这是一件非常简单的事情。

在你的情况下：

import numpy as np
import numba as nb
@nb.njit(fastmath=True)
def Test_1(X):
  K = np.zeros((B, B))
  for i in range(X.shape[0]):
      x_i = X[i, :]
      for j in range(X.shape[0]):
          x_j = X[j, :]
          K[i, j] = np.exp(-np.linalg.norm(x_i - x_j, 2) ** 2)

  return K

函数的并行化也很容易：

import numpy as np
import numba as nb
@nb.njit(fastmath=True,parallel=True)
def Test_1(X):
  K = np.zeros((B, B))
  for i in nb.prange(X.shape[0]):
      x_i = X[i, :]
      for j in range(X.shape[0]):
          x_j = X[j, :]
          K[i, j] = np.exp(-np.linalg.norm(x_i - x_j, 2) ** 2)

  return K

这很容易胜过迄今为止提供的所有其他解决方案。第一个函数调用大约需要 0.5 秒，因为这里你的代码被编译，但我猜你想多次调用这个函数。

如果使用单线程版本，还可以缓存编译结果。多线程代码的缓存可能很快就会实现。

Answer 3

我不确定您是否可以仅使用 numpy.我会使用 scipy 库中的 cdist 方法，如下所示：

import numpy as np 
from scipy.spatial.distance import cdist
B=5
X=np.random.rand(B*B).reshape((B,B))
dist = cdist(X, X, metric='euclidean')
K = np.exp(dist)

dist
array([[ 0.        ,  1.2659804 ,  0.98231231,  0.80089176,  1.19326493],
       [ 1.2659804 ,  0.        ,  0.72658078,  0.80618767,  0.3776364 ],
       [ 0.98231231,  0.72658078,  0.        ,  0.70205336,  0.81352455],
       [ 0.80089176,  0.80618767,  0.70205336,  0.        ,  0.60025858],
       [ 1.19326493,  0.3776364 ,  0.81352455,  0.60025858,  0.        ]])
K
array([[ 1.        ,  3.5465681 ,  2.67062441,  2.22752646,  3.29783084],
       [ 3.5465681 ,  1.        ,  2.06799756,  2.23935453,  1.45883242],
       [ 2.67062441,  2.06799756,  1.        ,  2.01789192,  2.25584482],
       [ 2.22752646,  2.23935453,  2.01789192,  1.        ,  1.82259002],
       [ 3.29783084,  1.45883242,  2.25584482,  1.82259002,  1.        ]])

希望对您有所帮助。干得好

编辑 您也可以只使用 numpy 数组，用于 theano 实现：

dist = (X ** 2).sum(1).reshape((X.shape[0], 1)) + (X ** 2).sum(1).reshape((1, X.shape[0])) - 2 * X.dot(X.T)

应该可以了！