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

Question

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

A 或 B 中没有重复的行。

我想知道在 Numpy 中是否有一种有效的方法来做到这一点？我找到了一个有点相关的答案：；但是，它 returns 实际行（不是布尔向量）。

Answer 1

您可以遵循 jterrace's answer, except use 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 连续数组）。所以，一般我们需要先用np.ascontiguous(A)再调用view.

如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)

所以在查看为 void 后使用 np.in1d 可能 return 不正确的数组结果使用 float 数据类型。

这里是一些建议方法的基准：

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 中间数组，它可能不会有效利用内存。

或者，如果值是 indices，即非负整数，您可以使用 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)

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

Matrix row difference, output a boolean vector

python

numpy

matrix

set-difference