Python 初学者。因此,我在尝试仅使用 numpy 库计算输入矩阵的行之间生成的二进制成对汉明顿距离矩阵时遇到麻烦。我应该避免循环并使用矢量化。例如,如果我有类似的东西:
[ 1, 0, 0, 1, 1, 0]
[ 1, 0, 0, 0, 0, 0]
[ 1, 1, 1, 1, 0, 0]
矩阵应该是这样的:
[ 0, 2, 3]
[ 2, 0, 3]
[ 3, 3, 0]
即,如果原始矩阵为 A,汉明距离矩阵为 B。B[0,1] = 汉明距离(A[0] 和 A[1])。在这种情况下,答案是 2,因为它们只有两个不同的元素。
所以我的代码是这样的
def compute_HammingDistance(X):
hammingDistanceMatrix = np.zeros(shape = (len(X), len(X)))
hammingDistanceMatrix = np.count_nonzero ((X[:,:,None] != X[:,:,None].T))
return hammingDistanceMatrix
但是它似乎只是返回一个标量值而不是预期的矩阵。我知道我可能在数组/矢量广播方面做错了什么,但我不知道如何修复它。我尝试过使用 np.sum 而不是 np.count_nonzero 但它们几乎都给了我类似的东西。
尝试这种方法,沿着
axis = 1
创建一个新轴,然后进行广播并使用 sum
: 计数 true 或非零值
(arr[:, None, :] != arr).sum(2)
# array([[0, 2, 3],
# [2, 0, 3],
# [3, 3, 0]])
def compute_HammingDistance(X):
return (X[:, None, :] != X).sum(2)
解释:
1)创建一个形状为 (3,1,6) 的 3d 数组
arr[:, None, :]
#array([[[1, 0, 0, 1, 1, 0]],
# [[1, 0, 0, 0, 0, 0]],
# [[1, 1, 1, 1, 0, 0]]])
2) 这是一个形状为 (3, 6) 的二维数组
arr
#array([[1, 0, 0, 1, 1, 0],
# [1, 0, 0, 0, 0, 0],
# [1, 1, 1, 1, 0, 0]])
3) 这会触发广播,因为它们的形状不匹配,并且 2d 数组 arr 首先沿着 3d 数组 arr[:, None, :] 的 0 轴广播,然后我们有形状数组 ( 1, 6) 针对 (3, 6) 进行广播。两个广播步骤一起对原始数组进行笛卡尔比较。
arr[:, None, :] != arr
#array([[[False, False, False, False, False, False],
# [False, False, False, True, True, False],
# [False, True, True, False, True, False]],
# [[False, False, False, True, True, False],
# [False, False, False, False, False, False],
# [False, True, True, True, False, False]],
# [[False, True, True, False, True, False],
# [False, True, True, True, False, False],
# [False, False, False, False, False, False]]], dtype=bool)
4) 沿第三轴的
sum
计算有多少个元素不相等,即给出汉明距离的 true。
出于我不明白的原因
(2 * np.inner(a-0.5, 0.5-a) + a.shape[1] / 2)
对于较大的数组来说,似乎比 @Psidom 快得多:
a = np.random.randint(0,2,(100,1000))
timeit(lambda: (a[:, None, :] != a).sum(2), number=100)
# 2.297890231013298
timeit(lambda: (2 * np.inner(a-0.5, 0.5-a) + a.shape[1] / 2), number=100)
# 0.10616962902713567
对于非常小的例子,Psidom 的速度要快一些:
a
# array([[1, 0, 0, 1, 1, 0],
# [1, 0, 0, 0, 0, 0],
# [1, 1, 1, 1, 0, 0]])
timeit(lambda: (a[:, None, :] != a).sum(2), number=100)
# 0.0004370050155557692
timeit(lambda: (2 * np.inner(a-0.5, 0.5-a) + a.shape[1] / 2), number=100)
# 0.00068191799800843
更新
部分原因似乎是浮点数比其他数据类型更快:
timeit(lambda: (0.5 * np.inner(2*a-1, 1-2*a) + a.shape[1] / 2), number=100)
# 0.7315902590053156
timeit(lambda: (0.5 * np.inner(2.0*a-1, 1-2.0*a) + a.shape[1] / 2), number=100)
# 0.12021801102673635
这是一个老话题,但这是我今天发现更正确的答案。
x = YourMatrix
Hd = x @ x.t()