将压缩的1D numpy数组转换为2D汉明距离矩阵

Question

我正在寻找一种可靠的方法，将使用scipy.spatial.distance.pdist函数生成的压缩汉明距离数组转换为相应的2D汉明距离矩阵。我知道scipy.spatial.distance.squareform函数，但是我需要计算高达100,000 x 100,000矩阵的汉明距离，这在Python中会导致MemoryError。 我需要以逐行的方式将压缩矩阵转换成方阵形式。是否有人知道使用NumPy和/或相关软件包的可靠（可能快速）实现呢？ 我需要对每行执行numpy.sum计算，但无法在内存中存储完整的N x N矩阵。 目前，我正在使用嵌套循环遍历输入矩阵并手动计算距离。

identity = 0.7
hamming_sum = numpy.zeros(msa_mat.shape[0], dtype=numpy.float64)
hamming_dist = numpy.zeros(msa_mat.shape[0], dtype=numpy.float64)
for i, row1 in enumerate(msa_mat):
    hamming_dist.fill(0)
    for j, row2 in enumerate(msa_mat):
        if i != j:
            hamming_dist[j] = scipy.spatial.distance.hamming(row1, row2)
    hamming_sum[i] = numpy.sum(numpy.where(hamming_dist < (1 - identity), 1, 0), axis=0)

---

编辑1

我的数据看起来像以下矩阵：

>>> a = numpy.array([1, 2, 3, 4, 5, 4, 5, 4, 2, 7, 9, 4, 1, 5, 6, 2, 3, 6], dtype=float).reshape(3, 6)
>>> a
array([[ 1.,  2.,  3.,  4.,  5.,  4.],
       [ 5.,  4.,  2.,  7.,  9.,  4.],
       [ 1.,  5.,  6.,  2.,  3.,  6.]])

我希望计算这个矩阵的海明距离。对于小矩阵，可以使用SciPy中的cdist命令轻松完成，并返回如下结果：

>>> cdist(a, a, 'hamming')
array([[ 0.        ,  0.83333333,  0.83333333],
       [ 0.83333333,  0.        ,  1.        ],
       [ 0.83333333,  1.        ,  0.        ]])

然而，在矩阵较大的情况下，这会在Python中引发MemoryError。

我知道可以使用“pdist”命令计算这些情况下的汉明距离。它返回一个1D数组中上三角的距离。

>>> pdist(a, 'hamming')
array([ 0.83333333,  0.83333333,  1.        ])

我的问题与我不知道如何在每一行上从结果重构矩阵有关。

我知道函数，但对于大型矩阵，它会引发MemoryErrors的问题。

- fsimkovic

---

pdist和cdist的输入是点数组。也就是说，如果输入是形状为(m, n)的数组，则表示n维空间中的m个点。在您的示例中，a的形状为(4, 4)：它表示4维空间中的4个点。所以我想知道为什么a恰好是对称的。这只是您选择的示例的偶然吗？ - Warren Weckesser

@WarrenWeckesser 您是正确的，这个例子选择得不好。我已经更新了帖子。 - fsimkovic

标题显示为"欧几里得距离"，而在答案的评论中您提到了欧几里得距离，但在问题正文中您说的是"汉明距离"，而示例代码计算的是汉明距离。那么您需要计算欧几里得距离的地方在哪里？ - Warren Weckesser

根据行hamming_sum[i] = numpy.sum(numpy.where(hamming_dist < (1 - identity), 1, 0), axis=0)，看起来你的最终目标是：对于每个点，计算在该点半径为1-identity内有多少其他点，使用汉明距离计算点之间的距离。正确吗？ - Warren Weckesser

Answer 1

以下是使用基于ID求和的方法，结合np.bincount函数的一种方法 -

def getdists_v1(a):
    n = a.shape[0]
    r,c = np.triu_indices(n,1)
    vals = pdist(a, 'hamming') < (1 - identity)
    return np.bincount(r,vals,minlength=n) + np.bincount(c,vals,minlength=n) + 1

这是另一个基于二进制的内容，注重内存效率，使用 np.add.reduceat。

def getdists_v2(a):
    n = a.shape[0]
    nr = (n*(n-1))//2
    vals = pdist(a, 'hamming') < (1 - identity)

    sfidx = n*np.arange(0,n-1) - np.arange(n-1).cumsum()
    id_arr = np.ones(nr,dtype=int)
    id_arr[sfidx[1:]] = -np.arange(n-3,-1,-1)
    c = id_arr.cumsum()

    out = np.bincount(c,vals)+1
    out[:n-1] += np.add.reduceat(vals,sfidx)
    return out

这是另一个循环计算下三角区域行求和的示例 -

def getdists_v3(a):
    n = a.shape[0]
    r_arr = np.arange(n-1)
    cr_arr = r_arr.cumsum()
    sfidx_c = (n-1)*r_arr - cr_arr
    vals = pdist(a, 'hamming') < (1 - identity)
    out = np.zeros(n)
    for i in range(n-1):
        out[i+1] = np.count_nonzero(vals[sfidx_c[:i+1] + i])
    out[:n-1] += np.add.reduceat(vals, n*r_arr - cr_arr)
    out[:] += 1
    return out

- Divakar

---

在之前的一个SO问题中，triu函数中的where对于大数组产生了内存错误。我会尝试查找一下。 - hpaulj

@Felix_Sim，你不需要使用np.where。看一下修改后的内容。 - Divakar

@hpaulj 对于上面那个来说很简单，但是对于下面那个可能需要一些额外的工作，并且考虑到内存限制，我不确定是否值得这样做。如果r和c可以在没有内存错误的情况下生成，我就不会担心它了。 - Divakar

@Felix_Sim 在哪个步骤？ - Divakar

在调用“hamming = numpy.bincount(c, vals) + 1”时，会引发“MemoryError”。 - fsimkovic

Answer 2

为了避免内存问题，一种方法是分批使用cdist：

import numpy as np
from scipy.spatial.distance import cdist


def count_hamming_neighbors(points, radius, batch_size=None):
    n = len(points)

    if batch_size is None:
        batch_size = min(n, 1000)

    hamming_sum = np.zeros(n, dtype=int)

    num_full_batches, last_batch = divmod(n, batch_size)
    batches = [batch_size]*num_full_batches
    if  last_batch != 0:
        batches.append(last_batch)
    for k, batch in enumerate(batches):
        i = batch_size*k
        dists = cdist(points[i:i+batch], points, metric='hamming')
        hamming_sum[i:i+batch] = (dists < radius).sum(axis=1)

    return hamming_sum

这里是与Divakar的getdists_v3(a)进行比较，以确保我们获得相同的结果：

In [102]: np.random.seed(12345)

In [103]: a = np.random.randint(0, 4, size=(16, 4))

In [104]: count_hamming_neighbors(a, 0.3)
Out[104]: array([1, 1, 3, 2, 2, 1, 2, 1, 3, 2, 3, 2, 2, 1, 2, 2])

In [105]: identity = 0.7

In [106]: getdists_v3(a)
Out[106]: 
array([ 1.,  1.,  3.,  2.,  2.,  1.,  2.,  1.,  3.,  2.,  3.,  2.,  2.,
        1.,  2.,  2.])

比较更大数组的时间：

In [113]: np.random.seed(12345)

In [114]: a = np.random.randint(0, 4, size=(10000, 4))

In [115]: %timeit hamming_sum = count_hamming_neighbors(a, 0.3)
1 loop, best of 3: 714 ms per loop

In [116]: %timeit v3result = getdists_v3(a)
1 loop, best of 3: 1.05 s per loop

因此，速度会稍微快一些。改变批次大小会影响性能，有时以出人意料的方式：

In [117]: %timeit hamming_sum = count_hamming_neighbors(a, 0.3, batch_size=250)
1 loop, best of 3: 643 ms per loop

In [118]: %timeit hamming_sum = count_hamming_neighbors(a, 0.3, batch_size=2000)
1 loop, best of 3: 875 ms per loop

In [119]: %timeit hamming_sum = count_hamming_neighbors(a, 0.3, batch_size=125)
1 loop, best of 3: 664 ms per loop

- Warren Weckesser

---

有没有办法自动确定最佳的“batch_size”？ - fsimkovic