矩阵行差，输出布尔值

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

可以将A和B视为两组XYZ数组，并使用^{}计算它们之间的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 2

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

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)

这对于整数数据类型是安全的。但是，请注意

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

{as}查看后可能返回不正确的结果 使用float dtype。你知道吗

---

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

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 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中间数组，但它可能不具有内存效率。你知道吗

或者，如果值是索引，即非负整数，则可以使用^{}将其缩减为一维数组。这是一种有效的算法，熊猫在内部用于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)