有效地按元素比较numpy数组与自身答案

【问题标题】：Comparing numpy array with itself by element efficiently有效地按元素比较numpy数组与自身
【发布时间】：2018-01-16 10:34:04
【问题描述】：

我正在执行大量这样的计算：

A == A[np.newaxis].T

其中 A 是一个密集的 numpy 数组，它经常具有共同的值。

出于基准测试目的，我们可以使用：

n = 30000
A = np.random.randint(0, 1000, n)
A == A[np.newaxis].T

当我执行这个计算时，我遇到了内存问题。我相信这是因为输出不是更有效的 bitarray 或 np.packedbits 格式。第二个问题是我们执行的比较次数是必要的两倍，因为生成的布尔数组是对称的。

我的问题是：

是否可以在不牺牲速度的情况下以更节省内存的方式生成布尔 numpy 数组输出？我知道的选项是 bitarray 和 np.packedbits，但我只知道在创建大型布尔数组之后如何应用它们。
我们能否利用计算的对称性将处理的比较次数减半，同时又不牺牲速度？

我需要能够执行 & 和 |对布尔数组输出的操作。我尝试过 bitarray，这对于这些按位运算来说非常快。但是打包 np.ndarray -> bitarray 然后解包 bitarray -> np.ndarray 很慢。

[已编辑以提供说明。]

【问题讨论】：

x == x[:, None] 似乎也能做你想做的事。
你的输出应该非常稀疏，所以我想说你想保存True 值的索引，而不是一些自制的packbits 实现
您可以argsort 并比较连续的元素。即 O(n log n) 时间和 O(n) 空间。与直接方法的 O(n^2) 和 O(n^2) 相比。
@DanielF，对稀疏矩阵执行布尔运算是否可能/有效？我试过使用 & 和 |操作，但这些似乎没有实现。
您希望能够对两个 30k x 30k 布尔矩阵执行 and 和 or 操作？或者你想以某种方式广播它们？

标签： python numpy

【解决方案1】：

这里有一个 numba 给我们一个 NumPy 布尔数组作为输出 -

from numba import njit

@njit
def numba_app1(idx, n, s, out):
    for i,j in zip(idx[:-1],idx[1:]):
        s0 = s[i:j]
        c = 0
        for p1 in s0[c:]:
            for p2 in s0[c+1:]:
                out[p1,p2] = 1
                out[p2,p1] = 1
            c += 1
    return out

def app1(A):
    s = A.argsort()
    b = A[s]
    n = len(A)
    idx = np.flatnonzero(np.r_[True,b[1:] != b[:-1],True])
    out = np.zeros((n,n),dtype=bool)
    numba_app1(idx, n, s, out)
    out.ravel()[::out.shape[1]+1] = 1
    return out

时间安排 -

In [287]: np.random.seed(0)
     ...: n = 30000
     ...: A = np.random.randint(0, 1000, n)

# Original soln
In [288]: %timeit A == A[np.newaxis].T
1 loop, best of 3: 317 ms per loop

# @Daniel F's soln-1 that skips assigning lower diagonal in output
In [289]: %timeit sparse_outer_eq(A)
1 loop, best of 3: 450 ms per loop

# @Daniel F's soln-2 (complete one)
In [291]: %timeit sparse_outer_eq(A)
1 loop, best of 3: 634 ms per loop

# Solution from this post
In [292]: %timeit app1(A)
10 loops, best of 3: 66.9 ms per loop

【讨论】：

嗯，不知道为什么你的时间和我的不一样。我做同样的事情，广播方法获得约 800 毫秒，稀疏应用获得约 500 毫秒。
@DanielF 我正在使用您的 - Alternatively, you can define sparse_outer_eq as: 代码。对吗？
不，我用的是第一个。
@DanielF 好吧，第一个sparse_outer_eq 看起来并不完整，因为它跳过了较低的诊断。你有完整的解决方案吗？
最后，我喜欢这个解决方案，因为生成的数组可以利用 numpy 的所有“好”功能，例如〜A。如果稀疏矩阵的实现足够通用，这样默认值可以非零，那就太好了。

【解决方案2】：

这甚至不是一个简单的答案，但应该可以通过使用一些自制的稀疏表示法来降低您的数据需求

from numba import jit

@jit   # because this is gonna be loopy
def sparse_outer_eq(A):
    n = A.size
    c = []
    for i in range(n):
        for j in range(i + 1, n):
            if A[i] == A[j]:
                 c.append((i, j))
    return c

现在c 是坐标元组(i, j)、i < j 的列表，它们对应于布尔数组中为“真”的坐标。您可以轻松地对这些 setwise 进行 and 和 or 操作：

list(set(c1) & set(c2))
list(set(c1) | set(c2))

稍后，当您想将此掩码应用于数组时，您可以取消坐标并将它们用于花哨的索引：

i_, j_ = list(np.array(c).T)
i = np.r_[i_, j_, np.arange(n)]
j = np.r_[j_, i_, np.arange(n)]

如果您关心订单，那么您可以np.lexsort i nd j

或者，您可以将sparse_outer_eq 定义为：

@jit
def sparse_outer_eq(A):
    n = A.size
    c = []
    for i in range(n):
        for j in range(n):
            if A[i] == A[j]:
                 c.append((i, j))
    return c

它保留了 >2x 的数据，但是坐标很简单：

 i, j = list(np.array(c).T)

如果你已经完成了任何set 操作，如果你想要一个合理的顺序，这仍然需要lexsorted。

如果您的坐标是每个 n 位整数，那么这应该比布尔格式更节省空间，只要您的稀疏度小于 1/n -> 32 位的 3% 左右。

至于时间，感谢numba 比广播还要快：

n = 3000
A = np.random.randint(0, 1000, n)

%timeit sparse_outer_eq(A)
100 loops, best of 3: 4.86 ms per loop

%timeit A == A[:, None]
100 loops, best of 3: 11.8 ms per loop

和比较：

a = A == A[:, None]

b = B == B[:, None]

a_ = sparse_outer_eq(A)

b_ = sparse_outer_eq(B)

%timeit a & b
100 loops, best of 3: 5.9 ms per loop

%timeit list(set(a_) & set(b_))
1000 loops, best of 3: 641 µs per loop

%timeit a | b
100 loops, best of 3: 5.52 ms per loop

%timeit list(set(a_) | set(b_))
1000 loops, best of 3: 955 µs per loop

编辑：如果您想做&~（根据您的评论），请使用第二个sparse_outer_eq 方法（这样您就不必跟踪对角线），然后就这样做：

list(set(a_) - set(b_))

【讨论】：

这很有用，但有几个问题：(1) 性能对唯一值的 n 和 # 非常敏感，例如对于 n = 15,000 和 200 个唯一值，我看到您的方法慢了 3 倍。 (2) 诸如~A 之类的操作并不那么容易。我可能必须存储一整套所有坐标 [或循环遍历它们] 并排除 set(a_) 中的项目。但是参见github.com/scipy/scipy/issues/1166.
是的，这种方法很大程度上依赖于这样一个事实，即所述的原始问题非常稀疏（1k 唯一值）。如果您的输入和输出矩阵相对密集，则布尔矩阵的性能会更高。
但是我要指出，即使在 1bit/boolean 下，只要稀疏度小于1/log2(n)，这种方法也比布尔矩阵更 memory 有效，即n~ 会破坏任何形式的稀疏性并使布尔值更有效。

【解决方案3】：

这里或多或少是规范的argsort 解决方案：

import numpy as np

def f_argsort(A):
    idx = np.argsort(A)
    As = A[idx]
    ne_ = np.r_[True, As[:-1] != As[1:], True]
    bnds = np.flatnonzero(ne_)
    valid = np.diff(bnds) != 1
    return [idx[bnds[i]:bnds[i+1]] for i in np.flatnonzero(valid)]

n = 30000
A = np.random.randint(0, 1000, n)
groups = f_argsort(A)

for grp in groups:
    print(len(grp), set(A[grp]), end=' ')
print()

【讨论】：

这很好，但就像 DanielF 的解决方案一样，它不容易允许 ~A 的计算/存储。公平地说，我应该在我的原始查询中明确说明这一点。但是，Divakar 的解决方案也使用了您的规范 argsort 函数，这就是我要在这里解决的问题。
@jp_data_analysis 别担心，选择最适合您的。

【解决方案4】：

我正在为我的问题添加一个解决方案，因为它满足以下 3 个属性：

低、固定、内存要求
快速按位运算（&、|、~ 等）
低存储空间，每个布尔值 1 位，通过打包整数

缺点是它以 np.packbits 格式存储。它比其他方法（尤其是 argsort）要慢得多，但如果速度不是问题，该算法应该运行良好。如果有人想出进一步优化的方法，这将非常有帮助。

更新：可以在此处找到以下算法的更高效版本：Improving performance on comparison algorithm np.packbits(A==A[:, None], axis=1)。

import numpy as np
from numba import jit

@jit(nopython=True)
def bool2int(x):
    y = 0
    for i, j in enumerate(x):
        if j: y += int(j)<<(7-i)
    return y

@jit(nopython=True)
def compare_elementwise(arr, result, section):
    n = len(arr)

    for row in range(n):
        for col in range(n):

            section[col%8] = arr[row] == arr[col]

            if ((col + 1) % 8 == 0) or (col == (n-1)):
                result[row, col // 8] = bool2int(section)
                section[:] = 0

    return result

A = np.random.randint(0, 10, 100)
n = len(A)
result_arr = np.zeros((n, n // 8 if n % 8 == 0 else n // 8 + 1)).astype(np.uint8)
selection_arr = np.zeros(8).astype(np.uint8)

packed = compare_elementwise(A, result_arr, selection_arr)

【讨论】：