numpy 删除具有唯一值的一维数组的公共值的最有效方法不是使用np.setDi1d和np.in1d

ubof19bj  于 2022-11-10  发布在  其他
关注(0)|答案(2)|浏览(143)

我需要一个更快的代码来删除一维数组(数组长度~10-15)的值,这些值与另一个一维数组(数组长度~1e5-5e5-->极少最高可达7e5)一样,是包含整数的索引数组。数组中没有重复,也没有排序,修改后值的顺序必须保持在主数组中。我知道可以使用这样的np.setdiff1dnp.in1d(这两个都不支持在no-python模式下运行的Numba)和其他类似的POST(例如this)来实现这一点,但性能在这里很重要,因为主索引数组中的所有值都将在循环中逐渐删除。

import numpy as np
import numba as nb

n = 500000
r = 10
arr1 = np.random.permutation(n)
arr2 = np.random.randint(0, n, r)

# @nb.jit

def setdif1d_np(a, b):
    return np.setdiff1d(a, b, assume_unique=True)

# @nb.jit

def setdif1d_in1d_np(a, b):
    return a[~np.in1d(a, b)]

还有一个由norok2提出的用于2D数组的相关帖子,即比那里描述的通常方法快约15倍(使用Numba的类似散列的方式)。如果可以为一维阵列做好准备,此解决方案可能是最好的:

@nb.njit
def mul_xor_hash(arr, init=65537, k=37):
    result = init
    for x in arr.view(np.uint64):
        result = (result * k) ^ x
    return result

@nb.njit
def setdiff2d_nb(arr1, arr2):
    # : build `delta` set using hashes
    delta = {mul_xor_hash(arr2[0])}
    for i in range(1, arr2.shape[0]):
        delta.add(mul_xor_hash(arr2[i]))
    # : compute the size of the result
    n = 0
    for i in range(arr1.shape[0]):
        if mul_xor_hash(arr1[i]) not in delta:
            n += 1
    # : build the result
    result = np.empty((n, arr1.shape[-1]), dtype=arr1.dtype)
    j = 0
    for i in range(arr1.shape[0]):
        if mul_xor_hash(arr1[i]) not in delta:
            result[j] = arr1[i]
            j += 1
    return result

我试着为一维阵列做准备,但我对此有一些问题。

  • 首先,IDU的mul_xor_hash到底是做什么的,如果initk是任意选择的
  • 为什么没有nb.njitmul_xor_hash将无法工作:
File "C:/Users/Ali/Desktop/test - Copy - Copy.py", line 21, in mul_xor_hash
    result = (result * k) ^ x
TypeError: ufunc 'bitwise_xor' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''
  • IDK如何在一维阵列上实现mul_xor_hash(如果可以),我猜这可能会使它比2DS更快,所以我通过[None, :]将输入数组广播到2D,这仅针对arr2得到以下错误:
print(mul_xor_hash(arr2[0]))
ValueError: new type not compatible with array
  • delta有什么作用

在这方面,我正在寻找最有效的方法。在没有比Norok2解决方案更好的方法的情况下,如何为一维阵列准备这个解决方案?

5fjcxozz

5fjcxozz1#

了解基于哈希的解决方案

首先,IDU的mul_xor_hash到底做了什么,以及是否任意选择了init和k
mul_xor_hash是一个定制的散列函数。已知混合XOR和MAXPLY(可能带有移位)的函数计算原始数据缓冲区的散列相对较快。乘法倾向于对位进行混洗,并且使用XOR以某种方式将结果组合/累加成固定大小的小值(即,最后的散列)。有许多不同的散列函数。一些比另一些更快,一些在给定的上下文中比另一些引起更多的冲突。导致太多冲突的快速散列函数在实践中可能是无用的,因为它将导致需要比较所有冲突值的病态情况。这就是为什么快速哈希函数很难实现的原因。
initk是参数,这肯定会导致散列非常平衡。这在这样的散列函数中很常见。k需要足够大,以便乘法对位进行混洗,并且它通常也应该是质数(由于模算术行为,像2的幂这样的值往往会增加冲突)。init仅在非常小的阵列(例如,对于1项):它通过将最终散列与一个非平凡常量进行异或来帮助减少冲突。事实上,如果arr.size = 1,那么result = (init * k) ^ arr[0],其中init * k是一个常量。众所周知,身份散列函数等于arr[0]是不好的,因为它往往会导致许多冲突(这是一个复杂的主题,但简而言之,arr[0]可以除以散列表中的桶的数量)。因此,init应该是一个相对较大的数字,而init * k也应该是一个很大的非平凡值(质数是一个很好的目标值)。
为什么没有nb.njit,mul_xor_hash将无法工作
这取决于输入。输入需要是一维数组,并具有可被8整除的原始字节大小(例如,64位项、2n x 32位项、4n x 16位项或8n个8位项)。以下是一些例子:

mul_xor_hash(np.random.rand(10))
mul_xor_hash(np.arange(10)) # Do not work with 9

达美航空是做什么的?
它是包含arr2行的散列的set,因此查找匹配行比不使用散列进行比较更快。
如何为一维阵列准备此解决方案?
AFAIK,哈希仅用于避免行的比较,但这是因为输入是2D数组。在1D中,没有这样的问题。
这种方法有一个很大的问题:它只有在没有散列冲突的情况下才起作用。否则,实现会错误地假设值是相等的,即使它们不是!@norok在评论中明确提到了这一点:
请注意,还应实现散列的冲突处理

更快实施

对于1D,使用@norok2的2D解决方案不是一个好主意,因为散列不会使其使用方式更快。事实上,无论如何,set已经在内部使用了哈希函数。更不用说需要正确实现冲突(这是由set完成的)。
使用set是一个相对较好的想法,因为它导致复杂性为O(n + m),其中n = len(arr1)m = len(arr2)。也就是说,如果将arr1转换为set,那么它将太大而无法放入一级缓存(由于您的情况是arr1的大小),从而导致较慢的缓存未命中。此外,set的不断增长的大小将导致重新散列值,这是非常低效的。如果将arr2转换为set,那么许多哈希表提取将不会非常高效,因为在您的例子中,arr2非常小。这就是为什么这个解决方案不是最优的。
一种解决方案是arr1拆分成块,然后基于目标块构建set。然后,您可以有效地检查某个值是否在集合中。由于规模不断扩大,建造布景的效率仍然不是很高。这个问题是由Python本身造成的,它不像其他语言那样提供为数据结构保留一些空间的方法(例如。C++)。避免这个问题的一种解决方案是简单地重新实现哈希表,这并不是微不足道和繁琐的。实际上,可以使用Bloom filters来加快这一过程,因为它们可以快速发现两个集合arr1arr2之间是否平均没有冲突(尽管它们的实现并不容易)。
另一种优化是使用多线程并行计算块,因为它们是独立的。也就是说,以并行方式高效地完成对最终数组的追加并非易事,尤其是因为您不希望修改顺序。一种解决方案是将副本从并行循环中移走并顺序执行,但这很慢,而且AFAIK目前在Numba中没有简单的方法来做到这一点(因为并行层非常有限)。考虑使用像C/C++这样的本地语言来实现高效的并行实现。

最后,与带有两个嵌套循环的简单实现相比,散列可能非常复杂,并且速度可能非常小,因为arr2只有很少的项,并且现代处理器可以使用SIMD指令快速比较值(而基于散列的方法在主流处理器上很难从这些指令中受益)。展开可以帮助编写一个非常简单快速的实现。同样不幸的是,Numba在内部使用了LLVM-Jit,这似乎无法向量化如此简单的代码(当然是由于LLVM-Jit甚至LLVM本身中的“缺少优化”)。结果,未矢量化的代码最终变慢了一些(而不是在现代主流处理器上快4到10倍)。一种解决方案是使用C/C++代码来实现这一点(或者可能使用Cython)。
以下是使用基本布隆过滤器的一系列实现:

@nb.njit('uint32(int32)')
def hash_32bit_4k(value):
    return (np.uint32(value) * np.uint32(27_644_437)) & np.uint32(0x0FFF)

@nb.njit(['int32[:](int32[:], int32[:])', 'int32[:](int32[::1], int32[::1])'])
def setdiff1d_nb_faster(arr1, arr2):
    out = np.empty_like(arr1)
    bloomFilter = np.zeros(4096, dtype=np.uint8)
    for j in range(arr2.size):
        bloomFilter[hash_32bit_4k(arr2[j])] = True
    cur = 0
    for i in range(arr1.size):
        # If the bloom-filter value is true, we know arr1[i] is not in arr2.
        # Otherwise, there is maybe a false positive (conflict) and we need to check to be sure.
        if bloomFilter[hash_32bit_4k(arr1[i])] and arr1[i] in arr2:
            continue
        out[cur] = arr1[i]
        cur += 1
    return out[:cur]

以下是一个未经测试的变体,它应该适用于64位整数(浮点数需要内存视图,可能还需要素数常量):

@nb.njit('uint64(int64)')
def hash_32bit_4k(value):
    return (np.uint64(value) * np.uint64(67_280_421_310_721)) & np.uint64(0x0FFF)

注意,如果小数组中的所有值都包含在每个循环的主数组中,那么我们可以通过在找到值时从arr2中删除这些值来加速arr1[i] in arr2部分。尽管如此,冲突和发现应该非常罕见,所以我预计这不会显著加快(更不用说它增加了一些开销和复杂性)。如果项目是按块计算的,那么最后的块可以直接复制,而不需要任何检查,但好处应该仍然相对较小。请注意,该策略对于前面提到的朴素(C/C++)SIMD实现是有效的(它可以快2倍左右)。

泛化与并行实现

本节重点介绍关于输入大小要使用的算法。它特别详细地介绍了基于SIMD的实现,并讨论了多线程的使用。
首先,对于值r,使用的最佳算法可以不同。更具体地说:

  • r为0时,最好的做法是原封不动地返回输入数组arr1(可能是副本,以避免就地算法问题);
  • r为1时,我们可以使用一个基本循环迭代数组,但最好的实现可能是使用经过高度优化的np.where的Numpy
  • r很小(如<10)时,使用基于SIMD的实现应该特别有效,特别是在基于arr2的循环的迭代范围在编译时已知并且展开的情况下
  • 对于仍然相对较小的较大r值(例如r < 1000r << n),提供的基于哈希的解决方案应该是最好的之一;
  • 对于r << n的较大r值,可以通过将布尔值打包为bloomFilter中的位,并使用多个散列函数而不是一个散列函数来优化基于散列的解决方案,以便更好地处理冲突,同时更有利于缓存(事实上,这就是实际的Bloom过滤器所做的);请注意,可以使用多线程,因此当r很大且r << n时,可以加快查找速度;
  • r很大并且不比n小很多时,这个问题很难有效解决,最好的解决方案当然是对两个数组进行排序(通常使用基数排序),并使用基于合并的方法来删除重复项,当rn都很大(难以实现)时,可能会使用多个线程。

让我们从基于SIMD的解决方案开始。以下是一个实现:

@nb.njit('int32[:](int32[::1], int32[::1])')
def setdiff1d_nb_simd(arr1, arr2):
    out = np.empty_like(arr1)
    limit = arr1.size // 4 * 4
    limit2 = arr2.size // 2 * 2
    cur = 0
    z32 = np.int32(0)

    # Tile (x4) based computation
    for i in range(0, limit, 4):
        f0, f1, f2, f3 = z32, z32, z32, z32
        v0, v1, v2, v3 = arr1[i], arr1[i+1], arr1[i+2], arr1[i+3]
        # Unrolled (x2) loop searching for a match in `arr2`
        for j in range(0, limit2, 2):
            val1 = arr2[j]
            val2 = arr2[j+1]
            f0 += (v0 == val1) + (v0 == val2)
            f1 += (v1 == val1) + (v1 == val2)
            f2 += (v2 == val1) + (v2 == val2)
            f3 += (v3 == val1) + (v3 == val2)
        # Remainder of the previous loop
        if limit2 != arr2.size:
            val = arr2[arr2.size-1]
            f0 += v0 == val
            f1 += v1 == val
            f2 += v2 == val
            f3 += v3 == val
        if f0 == 0: out[cur] = arr1[i+0]; cur += 1
        if f1 == 0: out[cur] = arr1[i+1]; cur += 1
        if f2 == 0: out[cur] = arr1[i+2]; cur += 1
        if f3 == 0: out[cur] = arr1[i+3]; cur += 1

    # Remainder
    for i in range(limit, arr1.size):
        if arr1[i] not in arr2:
            out[cur] = arr1[i]
            cur += 1

    return out[:cur]

事实证明,在我的机器上,这种实现总是比基于散列的实现慢,因为Numba显然为基于arr2的内部循环生成了一个低效的循环,这似乎来自与==相关的中断优化:Numba只是将SIMD指令用于此操作(没有明显的原因)。这阻止了许多替代SIMD相关代码的快速,只要它们使用Numba即可。
Numba的另一个问题是np.where速度很慢,因为它使用了一个简单的实现,而Numpy的实现已经进行了大量优化。由于之前的问题,在Numpy中进行的优化很难应用于Numba实现。这防止了在Numba代码中使用np.where的任何加速。
在实践中,基于散列的实现相当快,并且在我的机器上复制已经花费了相当多的时间。计算部分可以使用多线程来加速。这并不容易,因为Numba的并行模型非常有限。使用Numba(可以使用非临时存储,但Numba还不支持这一点)无法轻松优化副本,除非计算可能就地完成。
要使用多个线程,一种策略是首先将范围分成块,然后:

  • 构建布尔数组,确定arr1的每一项是否在arr2中找到(完全并行)

  • 统计按区块找到的项目数(完全并行)

  • 计算目标块的偏移量(很难并行,特别是用Numba,但由于块,速度很快)

  • 将区块复制到目标位置,而不复制找到的项目(完全并行)

以下是一种高效的基于散列的并行实现:

@nb.njit('int32[:](int32[:], int32[:])', parallel=True)
def setdiff1d_nb_faster_par(arr1, arr2):
    # Pre-computation of the bloom-filter
    bloomFilter = np.zeros(4096, dtype=np.uint8)
    for j in range(arr2.size):
        bloomFilter[hash_32bit_4k(arr2[j])] = True

    chunkSize = 1024 # To tune regarding the kind of input
    chunkCount = (arr1.size + chunkSize - 1) // chunkSize

    # Find for each item of `arr1` if the value is in `arr2` (parallel)
    # and count the number of item found for each chunk on the fly.
    # Note: thanks to page fault, big parts of `found` are not even written in memory if `arr2` is small
    found = np.zeros(arr1.size, dtype=nb.bool_)
    foundCountByChunk = np.empty(chunkCount, dtype=nb.uint16)
    for i in nb.prange(chunkCount):
        start, end = i * chunkSize, min((i + 1) * chunkSize, arr1.size)
        foundCountInChunk = 0
        for j in range(start, end):
            val = arr1[j]
            if bloomFilter[hash_32bit_4k(val)] and val in arr2:
                found[j] = True
                foundCountInChunk += 1
        foundCountByChunk[i] = foundCountInChunk

    # Compute the location of the destination chunks (sequential)
    outChunkOffsets = np.empty(chunkCount, dtype=nb.uint32)
    foundCount = 0
    for i in range(chunkCount):
        outChunkOffsets[i] = i * chunkSize - foundCount
        foundCount += foundCountByChunk[i]

    # Parallel chunk-based copy
    out = np.empty(arr1.size-foundCount, dtype=arr1.dtype)
    for i in nb.prange(chunkCount):
        srcStart, srcEnd = i * chunkSize, min((i + 1) * chunkSize, arr1.size)
        cur = outChunkOffsets[i]
        # Optimization: we can copy the whole chunk if there is nothing found in it 
        if foundCountByChunk[i] == 0:
            out[cur:cur+(srcEnd-srcStart)] = arr1[srcStart:srcEnd]
        else:
            for j in range(srcStart, srcEnd):
                if not found[j]:
                    out[cur] = arr1[j]
                    cur += 1

    return out

对于我的机器上的目标输入,这个实现是最快的。当n非常大并且在目标平台上创建线程的开销相对较小时(例如,在PC上,但通常不是具有许多核心的计算服务器)。并行实现的开销很大,因此目标机器上的核心数量需要至少为4个,这样实现才能显着快于顺序实现。
为目标输入调优chunkSize变量可能很有用。如果是r << n,最好使用较大的块大小。也就是说,块的数量需要足够大,以便多个线程在多个块上操作。因此,chunkSize应该明显小于n / numberOfThreads
在我的机器上,大部分时间(65%-70%)都花在最终副本上,这主要是内存受限的,很难用Numba进一步优化。

结果

以下是我的基于i5-9600KF的机器(6核)上的结果:

setdif1d_np:               2.65 ms
setdif1d_in1d_np:          2.61 ms
setdiff1d_nb:              2.33 ms
setdiff1d_nb_simd:         1.85 ms
setdiff1d_nb_faster:       0.73 ms
setdiff1d_nb_faster_par:   0.49 ms

最好的实现比其他实现快4~5倍。

yvgpqqbh

yvgpqqbh2#

我发现散列并没有帮助。对于2D情况,这只是一个技巧,将一维数组转换为单个数字并将其放入一个集合中。
下面是我将Norok2I转换为一维数组的方法(并添加了注解以加快编译速度)。请注意,这只比您已有的方法快一点(20%-30%)。当然,在第二次函数调用之后,由于编译的原因,第一次调用会稍微慢一些。

@nb.njit('int32[:](int32[:], int32[:])')
def setdiff1d_nb(arr1, arr2):
    delta = set(arr2)

    # : build the result
    result = np.empty(len(arr1), dtype=arr1.dtype)
    j = 0
    for i in range(arr1.shape[0]):
        if arr1[i] not in delta:
            result[j] = arr1[i]
            j += 1
    return result[:j]

相关问题