Ali*_*_Sh 3 python optimization performance numpy numba
我需要更快的代码来删除与另一个一维数组（数组长度〜1e5-5e5 -->很少达到7e5）常见的一维数组（数组长度〜10-15）的值，它们是包含整数的索引数组。数组中没有重复项，并且它们没有排序，并且修改后值的顺序必须保留在主数组中。我知道可以使用这样的np.setdiff1dor来实现np.in1d（两者都不支持在 no-python 模式下进行 numba jitted），而其他类似的帖子（例如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 数组提出了另一篇相关文章，该解决方案（使用 numba 的类似散列的方式）比那里描述的通常方法快约 15 倍。如果可以为一维数组准备，这个解决方案可能是最好的：
@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
我尝试为一维数组做准备，但我对此有一些问题/疑问。
mul_xor_hash到底做什么，以及 和init是否k是任意选择的mul_xor_hash没有以下内容就无法工作nb.njit： 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''
mul_xor_hash在 1D 数组上实现（如果可以的话），我想这可能比 2D 更快，所以我将输入数组广播到 2D by [None, :]，这仅导致以下错误arr2： print(mul_xor_hash(arr2[0]))
ValueError: new type not compatible with array
delta以及做什么我正在寻找这方面最有效的方法。在没有比norok2解决方案更好的方法的情况下，如何为一维数组准备该解决方案？
首先，IDU mul_xor_hash 到底做了什么，以及 init 和 k 是否任意选择
mul_xor_hash是一个自定义的哈希函数。众所周知，混合异或和乘法（可能带有移位）的函数计算原始数据缓冲区的哈希值相对较快。乘法往往会打乱位，而异或用于以某种方式将结果组合/累加为固定大小的小值（即最终散列）。有许多不同的散列函数。在给定的上下文中，有些比其他更快，有些比其他造成更多的冲突。导致太多冲突的快速散列函数在实践中可能是无用的，因为它会导致需要比较所有冲突值的病态情况。这就是快速哈希函数难以实现的原因。
init并且k参数肯定会导致哈希相当平衡。这在哈希函数中很常见。k需要足够大才能进行乘法来洗牌位，并且通常也应该是素数（由于模算术行为，像 2 的幂这样的值往往会增加冲突）。init仅对于非常小的数组（例如，只有 1 项）发挥重要作用：它通过将最终哈希值与非平凡常量进行异或来帮助减少冲突。事实上，如果arr.size = 1，那么result = (init * k) ^ arr[0]哪里init * k是一个常数。众所周知，身份哈希函数等于arr[0]是不好的，因为它往往会导致许多冲突（这是一个复杂的主题，但简单地说，arr[0]可以除以哈希表中的桶数）。因此，init应该是一个相对较大的数字，init * k也应该是一个很大的非平凡值（素数是一个很好的目标值）。
为什么 mul_xor_hash 在没有 nb.njit 的情况下无法工作
这取决于输入。输入必须是一维数组，并且具有可被 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
达美航空是做什么的
它set包含行的哈希值arr2，因此比不使用哈希值进行比较更快地找到匹配行。
如何为一维数组准备这个解决方案？
AFAIK，哈希仅用于避免行比较，但这是因为输入是二维数组。在一维中，不存在这样的问题。
这种方法有一个很大的问题：它只有在没有哈希冲突的情况下才有效。否则，即使值不相等，实现也会错误地假设它们相等！@norok 在评论中明确提到了这一点：
请注意，还应该实现散列的冲突处理
使用 @norok2 的 2D 解决方案来处理 1D 并不是一个好主意，因为散列不会使其使用速度更快。事实上，aset已经在内部使用了哈希函数。更不用说需要正确实现碰撞（这是由 a 完成的set）。
使用 aset是一个相对好的主意，因为它会导致复杂度为O(n + m)wheren = len(arr1)和m = len(arr2)。话虽这么说，如果arr1转换为 a set，那么它将太大而无法放入 L1 缓存（由于arr1您的情况的大小），导致缓慢的缓存未命中。此外，不断增长的大小set将导致值被重新散列，这是低效的。如果arr2转换为 a set，那么许多哈希表提取将不会非常有效，因为arr2在您的情况下非常小。这就是为什么这个解决方案不是最优的。
一种解决方案是分割arr1成块set，然后基于目标块构建。然后，您可以有效地检查某个值是否在集合中。由于规模不断增大，构建该集合的效率仍然不是很高。这个问题是由于Python本身没有提供像其他语言（例如C++）那样为数据结构保留一些空间的方法。避免此问题的一种解决方案是简单地重新实现哈希表，这并不简单且麻烦。实际上，布隆过滤器可以用来加速这个过程，因为它们可以快速发现两个集合之间是否没有冲突arr1以及arr2平均值（尽管它们实现起来并不简单）。
另一种优化是使用多个线程并行计算块，因为它们是独立的。话虽如此，并行有效地追加到最终数组并不容易，特别是因为您不希望修改顺序。一种解决方案是从并行循环中移走副本并串行执行，但这很慢，并且据我所知，目前在 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)
请注意，如果小数组中的所有值都包含在每个循环的主数组中，那么我们可以通过在找到它们时arr1[i] in arr2删除值来加快该部分的速度。arr2话虽这么说，碰撞和发现应该非常罕见，所以我不认为这会明显更快（更不用说它增加了一些开销和复杂性）。如果项目以块的形式计算，则可以直接复制最后的块而不进行任何检查，但好处仍然相对较小。请注意，此策略对于前面提到的简单 (C/C++) SIMD 实现非常有效（速度大约快 2 倍）。
本节重点介绍有关输入大小的算法。它特别详细介绍了基于 SIMD 的实现，并讨论了多线程的使用。
首先，关于值r，使用的最佳算法可能不同。进一步来说：
r为 0 时，最好的做法是返回arr1未修改的输入数组（可能是一个副本，以避免就地算法出现问题）；r为 1 时，我们可以使用一个基本循环迭代数组，但最好的实现可能是使用np.where对此进行了高度优化的 Numpyr很小，比如 <10 时，那么使用基于 SIMD 的实现应该特别高效，特别是如果arr2基于 的循环的迭代范围在编译时已知并且展开r仍然相对较小的较大值（例如r < 1000和r << n），提供的基于哈希的解决方案应该是最好的解决方案之一；r值r << n，可以通过将布尔值打包为位bloomFilter并使用多个哈希函数而不是一个哈希函数来优化基于哈希的解决方案，以便更好地处理冲突，同时更加缓存友好（事实上，这就是实际的绽放）过滤器确实如此）；r请注意，可以使用多线程，因此当很大并且时可以加快查找速度r << n；r很大并且不小于 太多时n，那么问题很难有效解决，最好的解决方案当然是对两个数组进行排序（通常使用基数排序）并使用基于合并的方法来删除重复项，可能使用多个当 和r都很大时线程n（难以实现）。让我们从基于 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 实现。np.where这可以防止在 Numba 代码中使用任何加速。
实际上，基于哈希的实现速度非常快，并且复制在我的机器上已经花费了大量时间。计算部分可以使用多线程来加速。这并不容易，因为 Numba 的并行模型非常有限。使用 Numba 无法轻松优化副本（可以使用非临时存储，但 Numba 尚不支持），除非计算可能就地完成。
要使用多个线程，一种策略是首先将范围分割为块，然后：
arr1是否在 中找到该项目arr2（完全并行）这是一个高效的基于并行哈希的实现：
@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。话虽这么说，块的数量需要足够大，以便多个线程能够对许多块进行操作。因此，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 倍。
| 归档时间： |
|
| 查看次数： |
8485 次 |
| 最近记录： |