在Python中对向量进行高性能的全对全比较

首先介绍一下背景:几种在聚类之间进行比较的方法依赖于所谓的对计数。我们在相同的n实体上有两个平面聚类向量a和b。在对所有可能的实体对进行配对计数时，我们检查它们是否在两个中属于同一簇，或者在a中相同但在b中不同，或者相反，或者在两个中都属于不同的簇。这样我们得到4个计数，我们称它们为n11, n10, n01, n00。这些是不同度量标准的输入。

当实体数量在10,000左右，并且要比较的聚类数量为数十或更多时，性能就会成为一个问题，因为每次比较的配对数量是10^8，而对于聚类的所有对所有比较，这需要执行10^4次。

对于一个简单的Python实现，它花了很长时间，所以我转向Cython和numpy。通过这种方式，我可以将一次比较的时间降低到0.9-3.0秒左右。对我来说，这仍然意味着半天的运行时间。我想知道你是否看到任何性能成就的可能性，用一些聪明的算法或C库，或其他什么。

以下是我的尝试:

1)在不为所有对分配巨大数组的情况下计数，取2个长度为n的成员向量a1, a2，并返回计数:

cimport cython
import numpy as np
cimport numpy as np
ctypedef np.int32_t DTYPE_t
@cython.boundscheck(False)
def pair_counts(
    np.ndarray[DTYPE_t, ndim = 1] a1,
    np.ndarray[DTYPE_t, ndim = 1] a2,
    ):
    cdef unsigned int a1s = a1.shape[0]
    cdef unsigned int a2s = a2.shape[0]
    cdef unsigned int n11, n10, n01, n00
    n11 = n10 = n01 = n00 = 0
    cdef unsigned int j0
    for i in range(0, a1s - 1):
        j0 = i + 1
        for j in range(j0, a2s):
            if a1[i] == a1[j] and a2[i] == a2[j]:
                n11 += 1
            elif a1[i] == a1[j]:
                n10 += 1
            elif a2[i] == a2[j]:
                n01 += 1
            else:
                n00 += 1
    return n11, n10, n01, n00

2)首先计算2个聚类中的每个聚类的成员向量(长度n * (n-1) / 2，每个实体对一个元素)，然后计算这些向量的计数。它在每次比较时分配~20-40M内存，但有趣的是，它比前一次比较快。注意:c是一个围绕集群的包装类，具有通常的成员向量，但也有一个c.members字典，其中包含numpy数组中每个集群的实体索引。

cimport cython
import numpy as np
cimport numpy as np
@cython.boundscheck(False)
def comembership(c):
    """
    Returns comembership vector, where each value tells if one pair
    of entites belong to the same group (1) or not (0).
    """
    cdef int n = len(c.memberships)
    cdef int cnum = c.cnum
    cdef int ri, i, ii, j, li
    cdef unsigned char[:] cmem = 
        np.zeros((int(n * (n - 1) / 2), ), dtype = np.uint8)
    for ci in xrange(cnum):
        # here we use the members dict to have the indices of entities
        # in cluster (ci), as a numpy array (mi)
        mi = c.members[ci]
        for i in xrange(len(mi) - 1):
            ii = mi[i]
            # this is only to convert the indices of an n x n matrix
            # to the indices of a 1 x (n x (n-1) / 2) vector:
            ri = n * ii - 3 * ii / 2 - ii ** 2 / 2 - 1
            for j in mi[i+1:]:
                # also here, adding j only for having the correct index
                li = ri + j
                cmem[li] = 1
    return np.array(cmem)
def pair_counts(c1, c2):
    p1 = comembership(c1)
    p2 = comembership(c2)
    n = len(c1.memberships)
    a11 = p1 * p2
    n11 = a11.sum()
    n10 = (p1 - a11).sum()
    n01 = (p2 - a11).sum()
    n00 = n - n10 - n01 - n11
    return n11, n10, n01, n00

3)这是一个纯粹基于numpy的解决方案，它创建了一个包含实体对成员关系的n x n布尔数组。输入是成员向量(a1, a2)。

def pair_counts(a1, a2):
    n = len(a1)
    cmem1 = a1.reshape([n,1]) == a1.reshape([1,n])
    cmem2 = a2.reshape([n,1]) == a2.reshape([1,n])
    n11 = int(((cmem1 == cmem2).sum() - n) / 2)
    n10 = int((cmem1.sum() - n) / 2) - n11
    n01 = int((cmem2.sum() - n) / 2) - n11
    n00 = n - n11 - n10 - n01
    return n11, n10, n01, n00

编辑:示例数据

import numpy as np
a1 = np.random.randint(0, 1868, 14440, dtype = np.int32)
a2 = np.random.randint(0, 484, 14440, dtype = np.int32)
# to have the members dicts used in example 2:
def get_cnum(a):
    """
    Returns number of clusters.
    """
    return len(np.unique(a))
def get_members(a):
    """
    Returns a dict with cluster numbers as keys and member entities
    as sorted numpy arrays.
    """
    members = dict(map(lambda i: (i, []), range(max(a) + 1)))
    list(map(lambda m: members[m[1]].append(m[0]),
        enumerate(a)))
    members = dict(map(lambda m:
       (m[0], np.array(sorted(m[1]), dtype = np.int)),   
       members.items()))
    return members
members1 = get_members(a1)
members2 = get_members(a2)
cnum1 = get_cnum(a1)
cnum2 = get_cnum(a2)

基于排序的方法有优点，但可以更简单地完成:

def pair_counts(a, b):
    n = a.shape[0]  # also b.shape[0]
    counts_a = np.bincount(a)
    counts_b = np.bincount(b)
    sorter_a = np.argsort(a)
    n11 = 0
    same_a_offset = np.cumsum(counts_a)
    for indices in np.split(sorter_a, same_a_offset):
        b_check = b[indices]
        n11 += np.count_nonzero(b_check == b_check[:,None])
    n11 = (n11-n) // 2
    n10 = (np.sum(counts_a**2) - n) // 2 - n11
    n01 = (np.sum(counts_b**2) - n) // 2 - n11
    n00 = n**2 - n - n11 - n10 - n01
    return n11, n10, n01, n00

如果这个方法在Cython中被有效地编码，那么就会获得另一个加速(可能是20倍)。

编辑:

我找到了一种方法来完全矢量化这个过程并降低时间复杂度:

def sizes2count(a, n):
    return (np.inner(a, a) - n) // 2
def pair_counts_vec_nlogn(a, b):
    # Size of "11" clusters (a[i]==a[j] & b[i]==b[j])
    ab = a * b.max() + b  # make sure max(a)*max(b) fits the dtype!
    _, sizes = np.unique(ab, return_counts=True)
    # Calculate the counts for each type of pairing
    n = len(a)  # also len(b)
    n11 = sizes2count(sizes, n)
    n10 = sizes2count(np.bincount(a), n) - n11
    n01 = sizes2count(np.bincount(b), n) - n11
    n00 = n**2 - n - n11 - n10 - n01
    return n11, n10, n01, n00
def pair_counts_vec_linear(a, b):
    # Label "11" clusters (a[i]==a[j] & b[i]==b[j])
    ab = a * b.max() + b
    # Calculate the counts for each type of pairing
    n = len(a)  # also len(b)
    n11 = sizes2count(np.bincount(ab), n)
    n10 = sizes2count(np.bincount(a), n) - n11
    n01 = sizes2count(np.bincount(b), n) - n11
    n00 = n**2 - n - n11 - n10 - n01
    return n11, n10, n01, n00

有时O(n log(n))算法比线性算法快，因为线性算法使用max(a)*max(b)存储。命名可能可以改进，我对术语不太熟悉。

比较A和B两个聚类在线性时间上的差异:

遍历A中的集群。设每个集群的大小为a_i。A中同一集群的pair总数为所有a_i*(a_i-1)/2的总和。
根据B中的集群对每个a集群进行分区。设每个分区的大小为b_j。A和B中同一集群的配对总数为所有b_j *(b_j-1)/2的总和。
两者之间的差值是A中属于同一簇但不属于B的对的总数
遍历B中的聚类，得到B中同一聚类的配对总数，从(2)中减去结果，得到B中同一聚类的配对，但不包括A。
以上3个结果的和是在A或B中相同的对的数量，减去n*(n-1)/2，得到在A和B中不同簇的对的总数

步骤(2)中的分区是通过为B创建一个映射item -> cluster的字典，然后查找每个a -cluster中的每个item来完成的。如果交叉比较很多聚类，可以为每个聚类只计算一次这些映射并保留它们，这样可以节省很多时间。

您不需要枚举和计数对。

相反，计算一个混淆矩阵，其中包含第一个聚类的每个聚类与第二个聚类的每个聚类的相交大小(这是对所有对象的一个循环)，然后使用公式n*(n-1)/2从该矩阵计算成对的数量。

这将你的运行时间从O(n^2)减少到O(n)，所以它应该给你一个可观的加速。

相关内容

最新更新

热门标签：