n0s*_*dow 5 python bit-manipulation numpy numpy-einsum
假设我有项目和一个代表这些项目包含在结果中的N多热值向量:{0, 1}
N = 4
# items 1 and 3 will be included in the result
vector = [0, 1, 0, 1]
# item 2 will be included in the result
vector = [0, 0, 1, 0]
Run Code Online (Sandbox Code Playgroud)
我还提供了一个冲突矩阵,指示哪些项目不能同时包含在结果中:
conflicts = [
[0, 1, 1, 0], # any result that contains items 1 AND 2 is invalid
[0, 1, 1, 1], # any result that contains AT LEAST 2 items from {1, 2, 3} is invalid
]
Run Code Online (Sandbox Code Playgroud)
给定这个冲突矩阵,我们可以确定早期 s 的有效性vector:
# invalid as it triggers conflict 1: [0, 1, 1, 1]
vector = [0, 1, 0, 1]
# valid as it triggers no conflicts
vector = [0, 0, 1, 0]
Run Code Online (Sandbox Code Playgroud)
检测给定是否vector“有效”(即不触发任何冲突)的简单解决方案可以通过 numpy 中的点积和求和运算来完成:
violation = np.dot(conflicts, vector)
is_valid = np.max(violation) <= 1
Run Code Online (Sandbox Code Playgroud)
是否有更有效的方法来执行此操作,也许通过np.einsum或完全绕过 numpy 数组以支持位操作?
我们假设被检查的向量的数量可能非常大(例如,2^N如果我们评估所有可能性),但一次可能只检查一个向量(以避免生成形状高达作为(2^N, N)输入的矩阵)。
TL;DR:您可以使用Numba进行优化,只对二进制值np.dot进行操作。更具体地说,您可以使用64 位视图一次对 8 个字节执行类似 SIMD 的操作。
首先,使用这种方法可以将列表有效地转换为相对紧凑的数组:
vector = np.fromiter(vector, np.uint8)
conflicts = np.array([np.fromiter(conflicts[i], np.uint8) for i in range(len(conflicts))])
Run Code Online (Sandbox Code Playgroud)
这比使用自动 Numpy 转换更快,或者np.array(在 Numpy 代码内部执行的检查更少,Numpy 知道要构建什么类型的数组,并且生成的数组在内存中更小,因此填充速度更快)。此步骤可用于加速np.dot基于 的解决方案。
如果输入已经是 Numpy 数组,则检查它们的类型为np.uint8或np.int8。否则,请使用conflits = conflits.astype(np.uint8)例如将它们转换为此类类型。
然后,一种解决方案可能是np.packbits将输入二进制值尽可能地打包在内存中的位数组中,然后执行逻辑与。但事实证明这np.packbits相当慢。因此,这个解决方案最终并不是一个好主意。事实上,任何创建形状类似的临时数组的解决方案conflicts都会很慢,因为在内存中写入这样的数组通常比np.dot(从内存中读取conflicts一次)慢。
由于np.dot经过了很好的优化,击败它的唯一解决方案是使用优化的本机代码。借助即时编译器,Numba 可用于在运行时从基于 Numpy 的 Python 代码生成本机可执行代码。vector这个想法是在每个块的 和 行之间执行逻辑与conflicts。检查每个块的冲突,以便尽早停止计算。通过比较两个数组的 uint64 视图(以 SIMD 友好的方式),可以按 8 个八位字节组有效地比较块。
import numba as nb
@nb.njit('bool_(uint8[::1], uint8[:,::1])')
def check_valid(vector, conflicts):
n, m = conflicts.shape
assert vector.size == m
for i in range(n):
block_size = 128 # In the range: 8,16,...,248
conflicts_row = conflicts[i,:]
gsum = 0 # Global sum of conflicts
m_limit = m // block_size * block_size
for j in range(0, m_limit, block_size):
vector_block = vector[j:j+block_size].view(np.uint64)
conflicts_block = conflicts_row[j:j+block_size].view(np.uint64)
# Matching
lsum = np.uint64(0) # 8 local sums of conflicts
for k in range(block_size//8):
lsum += vector_block[k] & conflicts_block[k]
# Trick to perform the reduction of all the bytes in lsum
lsum += lsum >> 32
lsum += lsum >> 16
lsum += lsum >> 8
gsum += lsum & 0xFF
# Check if there is a conflict
if gsum >= 2:
return False
# Remaining part
for j in range(m_limit, m):
gsum += vector[j] & conflicts_row[j]
if gsum >= 2:
return False
return True
Run Code Online (Sandbox Code Playgroud)
对于大量形状(没有冲突),这比我的机器上快大约 9 倍。这两种情况都不包括转换列表的时间。当存在冲突时,所提供的解决方案要快得多,因为它可以提前停止计算。np.dotconflicts(16, 65536)
从理论上讲,计算应该更快,但 Numba JIT 无法成功使用 SIMD 指令对循环进行矢量化。话虽这么说,似乎同样的问题也出现在np.dot. 如果数组更大,您可以并行计算块(如果函数返回 False,则计算速度会变慢)。
| 归档时间: |
|
| 查看次数: |
169 次 |
| 最近记录: |