矩阵行差，输出一个布尔向量

Question

我有一个m x 3 矩阵A 及其行子集B (n x 3)。两者都是指向另一个大型 4D 矩阵的索引集；他们的数据类型是dtype('int64')。我想生成一个布尔向量x，其中x[i] = True 如果B 不包含行A[i,:]。

A 或 B 中没有重复行。

我想知道在 Numpy 中是否有一种有效的方法可以做到这一点？我找到了一个有点相关的答案：https://stackoverflow.com/a/11903368/265289;但是，它返回实际行（不是布尔向量）。

Answer 1

您可以遵循jterrace's answer 中所示的相同模式，除了使用np.in1d 而不是np.setdiff1d：

import numpy as np
np.random.seed(2015)

m, n = 10, 5
A = np.random.randint(10, size=(m,3))
B = A[np.random.choice(m, n, replace=False)]
print(A)
# [[2 2 9]
#  [6 8 5]
#  [7 8 0]
#  [6 7 8]
#  [3 8 6]
#  [9 2 3]
#  [1 2 6]
#  [2 9 8]
#  [5 8 4]
#  [8 9 1]]

print(B)
# [[2 2 9]
#  [1 2 6]
#  [2 9 8]
#  [3 8 6]
#  [9 2 3]]

def using_view(A, B, assume_unique=False):
    Ad = np.ascontiguousarray(A).view([('', A.dtype)] * A.shape[1])
    Bd = np.ascontiguousarray(B).view([('', B.dtype)] * B.shape[1])
    return ~np.in1d(Ad, Bd, assume_unique=assume_unique)

print(using_view(A, B, assume_unique=True))

产量

[False  True  True  True False False False False  True  True]

您可以使用assume_unique=True（可以加快计算速度）因为 A 或 B 中没有重复行。

---

注意A.view(...) 会引发

ValueError: new type not compatible with array.

如果A.flags['C_CONTIGUOUS'] 是False（即如果A 不是C 连续数组）。 因此，一般我们需要在调用view之前使用np.ascontiguous(A)。

---

作为 B.M.建议，您可以改为使用 "void" dtype 查看每一行：

def using_void(A, B):
    dtype = 'V{}'.format(A.dtype.itemsize * A.shape[-1])
    Ad = np.ascontiguousarray(A).view(dtype)
    Bd = np.ascontiguousarray(B).view(dtype)
    return ~np.in1d(Ad, Bd, assume_unique=True)

这对于整数 dtype 是安全的。但是，请注意

In [342]: np.array([-0.], dtype='float64').view('V8') == np.array([0.], dtype='float64').view('V8')
Out[342]: array([False], dtype=bool)

因此在查看为 void 后使用 np.in1d 可能会返回不正确的数组结果 带有浮点数据类型。

---

以下是一些建议方法的基准：

import numpy as np
np.random.seed(2015)

m, n = 10000, 5000
# Note A may contain duplicate rows, 
# so don't use assume_unique=True for these benchmarks. 
# In this case, using assume_unique=False does not improve the speed much anyway.
A = np.random.randint(10, size=(2*m,3))
# make A not C_CONTIGUOUS; the view methods fail for non-contiguous arrays
A = A[::2]  
B = A[np.random.choice(m, n, replace=False)]

def using_view(A, B, assume_unique=False):
    Ad = np.ascontiguousarray(A).view([('', A.dtype)] * A.shape[1])
    Bd = np.ascontiguousarray(B).view([('', B.dtype)] * B.shape[1])
    return ~np.in1d(Ad, Bd, assume_unique=assume_unique)

from scipy.spatial import distance
def using_distance(A, B):
    return ~np.any(distance.cdist(A,B)==0,1)

from functools import reduce 
def using_loop(A, B):
    pred = lambda i: A[:, i:i+1] == B[:, i]
    return ~reduce(np.logical_and, map(pred, range(A.shape[1]))).any(axis=1)

from pandas.core.groupby import get_group_index, _int64_overflow_possible
from functools import partial
def using_pandas(A, B):
    shape = [1 + max(A[:, i].max(), B[:, i].max()) for i in range(A.shape[1])]
    assert not _int64_overflow_possible(shape)

    encode = partial(get_group_index, shape=shape, sort=False, xnull=False)
    a1, b1 = map(encode, (A.T, B.T))
    return ~np.in1d(a1, b1)

def using_void(A, B):
    dtype = 'V{}'.format(A.dtype.itemsize * A.shape[-1])
    Ad = np.ascontiguousarray(A).view(dtype)
    Bd = np.ascontiguousarray(B).view(dtype)
    return ~np.in1d(Ad, Bd)

# Sanity check: make sure all the functions return the same result
for func in (using_distance, using_loop, using_pandas, using_void):
    assert (func(A, B) == using_view(A, B)).all()

---

In [384]: %timeit using_pandas(A, B)
100 loops, best of 3: 1.99 ms per loop

In [381]: %timeit using_void(A, B)
100 loops, best of 3: 6.72 ms per loop

In [378]: %timeit using_view(A, B)
10 loops, best of 3: 35.6 ms per loop

In [383]: %timeit using_loop(A, B)
1 loops, best of 3: 342 ms per loop

In [379]: %timeit using_distance(A, B)
1 loops, best of 3: 502 ms per loop

Answer 2

您可以将A 和B 视为两组XYZ 数组，并使用scipy.spatial.distance.cdist 计算它们之间的euclidean distances。我们会对零距离感兴趣。这个距离计算应该是一个非常有效的实现，所以希望我们能有一个有效的解决方案来解决我们的问题。因此，查找此类布尔输出的实现将如下所示 -

from scipy.spatial import distance

out = ~np.any(distance.cdist(A,B)==0,1)
# OR np.all(distance.cdist(A,B)!=0,1)

示例运行 -

In [582]: A
Out[582]: 
array([[0, 2, 2],
       [1, 0, 3],
       [3, 3, 3],
       [2, 0, 3],
       [2, 0, 1],
       [1, 1, 1]])

In [583]: B
Out[583]: 
array([[2, 0, 3],
       [2, 3, 3],
       [1, 1, 3],
       [2, 0, 1],
       [0, 2, 2],
       [2, 2, 2],
       [1, 2, 3]])

In [584]: out
Out[584]: array([False,  True,  True, False, False,  True], dtype=bool)

Answer 3

由于只有 3 列，一种解决方案是跨列减少：

>>> a
array([[2, 2, 9],
       [6, 8, 5],
       [7, 8, 0],
       [6, 7, 8],
       [3, 8, 6],
       [9, 2, 3],
       [1, 2, 6],
       [2, 9, 8],
       [5, 8, 4],
       [8, 9, 1]])
>>> b
array([[2, 2, 9],
       [1, 2, 6],
       [2, 9, 8],
       [3, 8, 6],
       [9, 2, 3]])

>>> from functools import reduce
>>> pred = lambda i: a[:, i:i+1] == b[:,i]
>>> reduce(np.logical_and, map(pred, range(a.shape[1]))).any(axis=1)
array([ True, False, False, False,  True,  True,  True,  True, False, False], dtype=bool)

虽然这会创建一个m x n 中间数组，但内存效率可能不高。

或者，如果值是索引，即非负整数，您可以使用pandas.groupby.get_group_index 减少到一维数组。这是 pandas 在内部用于groupby 操作的高效算法；唯一需要注意的是，您可能需要验证不会有任何整数溢出：

>>> from pandas.core.groupby import get_group_index, _int64_overflow_possible
>>> from functools import partial

>>> shape = [1 + max(a[:, i].max(), b[:, i].max()) for i in range(a.shape[1])]
>>> assert not _int64_overflow_possible(shape)

>>> encode = partial(get_group_index, shape=shape, sort=False, xnull=False)
>>> a1, b1 = map(encode, (a.T, b.T))
>>> np.in1d(a1, b1)
array([ True, False, False, False,  True,  True,  True,  True, False, False], dtype=bool)