Numpy没有yet的基数排序,所以我想知道是否可以使用预先存在的numpy函数编写一个。到目前为止,我已经有了下面的代码,它确实可以工作,但是比numpy的quicksort慢了大约10倍。
测试和基准:
a = np.random.randint(0, 1e8, 1e6)
assert(np.all(radix_sort(a) == np.sort(a)))
%timeit np.sort(a)
%timeit radix_sort(a)
mask_b
循环至少可以部分向量化,从&
跨掩码广播,并使用cumsum
和axis
arg,但这最终是一个pessimization,可能是由于内存占用增加。
如果有人能看到一种方法来改进我所拥有的东西,我会很有兴趣听到,即使它仍然比np.sort
慢.
请注意,您可以很容易地can implement一个快速计数排序,尽管这只与小整数数据有关。
**编辑1:**将np.arange(n)
从循环中取出会有一点帮助,但这并不是很令人兴奋。
编辑2:cumsum
实际上是冗余的(哦!),但这个更简单的版本只对性能有轻微的帮助。
def radix_sort(a):
bit_len = np.max(a).bit_length()
n = len(a)
cached_arange = arange(n)
idx = np.empty(n, dtype=int) # fully overwritten each iteration
for mask_b in xrange(bit_len):
is_one = (a & 2**mask_b).astype(bool)
n_ones = np.sum(is_one)
n_zeros = n-n_ones
idx[~is_one] = cached_arange[:n_zeros]
idx[is_one] = cached_arange[:n_ones] + n_zeros
# next three lines just do: a[idx] = a, but correctly
new_a = np.empty(n, dtype=a.dtype)
new_a[idx] = a
a = new_a
return a
**编辑3:**如果你在多个步骤中构造idx,你可以一次循环两个或更多个比特,而不是循环单个比特。使用2位有点帮助,我没有尝试更多:
idx[is_zero] = np.arange(n_zeros)
idx[is_one] = np.arange(n_ones)
idx[is_two] = np.arange(n_twos)
idx[is_three] = np.arange(n_threes)
**编辑4和5:**对于我正在测试的输入,转到4位似乎是最好的。此外,您可以完全摆脱idx
步骤。现在只比np.sort
(source available as gist)慢5倍,而不是10倍:
**编辑6:**这是上面的整理版本,但它也有点慢 *。80%的时间花在repeat
和extract
上-如果有一种方法可以广播extract
:(.
def radix_sort(a, batch_m_bits=3):
bit_len = np.max(a).bit_length()
batch_m = 2**batch_m_bits
mask = 2**batch_m_bits - 1
val_set = np.arange(batch_m, dtype=a.dtype)[:, nax] # nax = np.newaxis
for _ in range((bit_len-1)//batch_m_bits + 1): # ceil-division
a = np.extract((a & mask)[nax, :] == val_set,
np.repeat(a[nax, :], batch_m, axis=0))
val_set <<= batch_m_bits
mask <<= batch_m_bits
return a
**编辑7和8:**实际上,您可以使用numpy.lib.stride_tricks
从as_strided
广播提取,但它似乎对性能没有多大帮助:
最初,这对我来说是有意义的,因为extract
将在整个数组batch_m
上迭代,所以CPU请求的缓存行总数将与之前相同(只是在过程结束时,它已经请求了每个缓存行batch_m
次)。然而,实际情况是extract
不够聪明,无法覆盖任意阶跃数组,必须在开始之前扩展数组,即无论如何重复结束。事实上,在查看了extract
的源代码后,我现在看到我们可以用这种方法做的最好的事情是:
a = a[np.flatnonzero((a & mask)[nax, :] == val_set) % len(a)]
比extract
稍慢。然而,如果len(a)
是2的幂,我们可以用& (len(a) - 1)
代替昂贵的mod操作,这确实比extract
版本快一点(现在a=randint(0, 1e8, 2**20
大约是4.9x np.sort
)。我想我们可以通过零填充来实现非二次幂长度,然后在排序结束时裁剪额外的零.然而,这将是一种悲观,除非长度已经接近二次幂。
3条答案
按热度按时间yuvru6vn1#
我试过用Numba来看看基数排序有多快。Numba良好性能的关键(通常)是写出所有循环,这是非常有指导意义的。最后,我得到了以下结果:
从4个内部循环中,很容易看出这是第4个,这使得Numpy很难进行矢量化。
绕过这个问题的一种方法是从Scipy中引入一个特定的C++函数:
scipy.sparse.coo.coo_tocsr
。它的内部循环与上面的Python函数几乎相同,因此可以滥用它来在Python中编写更快的“向量化”基数排序。也许是这样的:像上面这样的LSB基数排序的一个低效之处是数组在RAM中被完全打乱了很多次,这意味着CPU缓存没有得到很好的利用。为了尝试减轻这种影响,可以选择首先使用MSB基数排序进行一次传递,将项目大致放在RAM的正确块中,然后使用LSB基数排序对每个结果组进行排序。下面是一个实现:
计时(在我的系统上每个都有最佳性能设置):
Numba表现得很好,正如预期的那样。而且,通过对现有C扩展的一些巧妙应用,有可能击败
numpy.sort
。IMO在优化的层面上,你已经得到了它的好处,它也考虑Numpy的附加组件,但我不会真正考虑我的答案中的实现“矢量化”:大部分工作是在外部专用功能中完成的。另一件让我印象深刻的事情是对基数选择的敏感性。对于我尝试的大多数设置,我的实现仍然比
numpy.sort
慢,因此在实践中需要某种启发式方法来提供全面的良好性能。daolsyd02#
你能把它改成一次8位的计数/基数排序吗?对于32位无符号整数,创建一个字节字段出现计数的矩阵[4][257],对要排序的数组进行一次读取。矩阵[][0]= 0,矩阵[][1]= 0的占用数,.然后将计数转换成索引,其中矩阵[][0]= 0,矩阵[][1]=字节数== 0,矩阵[][2]=字节数== 0+字节数= 1,.不使用最后一个计数,因为这将索引数组的末尾。然后进行4次基数排序,在原始数组和输出数组之间来回移动数据。一次处理16位需要一个矩阵[2][65537],但只需要2遍。示例C代码:
yhxst69z3#
我实际上使用Cython创建了一个基数排序。从我的测试来看,它比Rust或C/C++中的基数实现快了大约5%。而且它比np. sort快了很多(我想大概是10%)。地址:https://github.com/Ohmagar/Radix_cython/blob/main/parallel_radix_5.pyx
我做了一些漂亮的事情来减少处理,通过number_of_digits对元素进行预排序,确保数字只在有数字要排序时才被排序到桶中。因此,“10”将只被处理两次而不是8次(如果9_999_999 < max_element < 10_000_000)。我从头开始在python中构建它作为POC,并对它做了越来越多的工作。一旦我无法获得更多的速度,我就在Cython中重写了它,并开始进行更多的修补。最后一步是对每个“digit_chunk”进行并行处理,这最终导致我的实现比任何可比的东西都快,特别是numpy. sort。
我刚刚看到我可以让它更快一点,可能在处理器功能中,通过让预排序也并行完成。不知道我怎么没发现。
你可以去看看。