Python - 如何生成成对汉明距离矩阵

6

我是Python的初学者。我正在尝试使用仅使用numpy库来计算输入矩阵行之间的二进制对比Hamming距离矩阵,但遇到了困难。我需要避免使用循环,并进行向量化处理。例如,如果我有以下内容:

   [ 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]=hammingdistance(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,但它们基本上都给我类似的结果。

2个回答

9
尝试这个方法,沿着axis = 1创建一个新的轴,然后使用广播和计算sum中的真值或非零值:
(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)的三维数组。

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的数量,这就是汉明距离。


我用 hammingDistanceMatrix = np.count_nonzero((X[:, None, :] != X).sum(2)) 替换了代码中原来的 i,但结果似乎仍然是一个标量值。 - user2444400
不需要使用np.count_nonzero。使用sum函数即可完成操作。只需返回 hammingDistanceMatrix = (arr[:, None, :] != arr).sum(2) 即可。 - Psidom
哦,好的!非常感谢。您能否解释一下为什么广播是这样完成的arr[:,None,:],以及为什么要使用.sum(2)? - user2444400

5

由于我不理解原因

(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

这是一种聪明的比较元素并一次性计数的方法。我认为浮点数数据类型在这种情况下更快的主要原因是BLAS例程仅适用于这些数据类型。所有其他数据类型都由不太优化的numpy C代码处理。 - user7138814
哇,我不知道BLAS会有这么大的差异! - Paul Panzer

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接