原地基数排序

212

这是一段很长的文本,请耐心阅读。简而言之,问题是:是否有可行的原地基数排序算法


初步

我有很多只使用字母“A”、“C”、“G”和“T”(是的,你猜对了:DNA)的小固定长度字符串需要排序。

目前,我使用std::sort,它在所有常见的STL实现中使用introsort。这个方法效果还不错。然而,我相信基数排序完全适合我的问题集,并且在实践中应该表现得更好。

细节

我用一个非常天真的实现来测试了这个假设,对于相对较小的输入(大约为10,000),这是正确的(至少快两倍)。然而,当问题规模变得更大(N > 5,000,000)时,运行时间会极度下降。

原因很明显:基数排序需要复制整个数据(实际上在我的天真实现中不止一次)。这意味着我已经将~4 GiB放入了我的主内存中,这显然会影响性能。即使没有这个问题,我也无法承受使用这么多内存,因为问题的规模实际上会更大。

用例

Ideally,这个算法应该适用于2到100个字符长度的DNA和DNA5(允许使用额外的通配符“N”),甚至包括带有IUPAC ambiguity codes的DNA(导致16个不同的值)。然而,我意识到并不是所有情况都能覆盖,所以我很高兴能获得任何速度上的提升。代码可以动态地决定派发哪种算法。

研究

不幸的是,Wikipedia article on radix sort对原地变体的部分完全没有用处。NIST-DADS section on radix sort几乎不存在。有一篇听起来很有前途的论文,叫做Efficient Adaptive In-Place Radix Sorting,描述了算法“MSL”。不幸的是,这篇论文也令人失望。
特别是以下几点。

首先,该算法存在几个错误并留下了很多未解释的问题。特别是,它没有详细说明递归调用(我只是假设它增加或减少某些指针来计算当前的移位和掩码值)。此外,它在不给出定义的情况下使用函数dest_groupdest_address。我无法看出如何高效地实现这些(也就是说,在O(1)中;至少dest_address不是微不足道的)。

最后,该算法通过交换数组索引与输入数组内部的元素来实现原地性。这显然仅适用于数字数组。我需要在字符串上使用它。当然,我可以放弃强类型并继续假设内存将容忍我在不属于它的位置存储索引。但是,只要我能将我的字符串压缩到32位内存中(假设32位整数),这种方法才有效。那只有16个字符(暂且忽略16>log(5,000,000))。

该作者的另一篇论文根本没有提供准确的描述,但它将MSL的运行时间描述为次线性,这是完全错误的。

总之:是否有希望找到一个适用于DNA字符串的工作参考实现或至少是工作在原地的基数排序的好伪代码/描述?


73
这是一个写得非常好的问题。 - JustinT
1
小固定长度字符串有多小? - EvilTeach
2
@Stephan:这一切都很好。但是在复制/缓存未命中的情况下,我只会遇到延迟。而在内存方面,我会遇到物理限制。这是绝对不可谈判的。所有那些将数据的部分存储在磁盘上的花哨技术肯定比当前的快速排序解决方案慢。 - Konrad Rudolph
2
dsimcha的解决方案则在某些输入情况下绝对比quicksort更快。移动次数可能会很高,缓存局部性也较小,但在现实世界中仍然表现良好。我还稍微调整了一下解决方案,以减少需要执行的交换次数。 - Konrad Rudolph
1
@PeterMortensen 对于未来,我欣赏您提供更正和链接以增加上下文。然而,我并不特别欣赏那些仅仅是风格问题的编辑(如“在...次序中”与“按...顺序”,“即”与“也就是”)。 - Konrad Rudolph
显示剩余7条评论
15个回答

2
看起来您已经解决了问题,但是记录一下,似乎可行的原地基数排序版本之一是“美国国旗排序”。它在这里描述:工程基数排序。总体思路是对每个字符进行两次处理-首先计算每个字符的数量,以便将输入数组分成段。然后再次遍历,将每个元素交换到正确的段中。现在在下一个字符位置上递归地对每个段进行排序。

实际上,我使用的解决方案与旗帜排序算法非常密切相关。我不知道是否有任何相关区别。 - Konrad Rudolph
2
从来没有听说过美国国旗排序(American Flag Sort),但显然这就是我编写的算法:http://coliru.stacked-crooked.com/a/94eb75fbecc39066。目前,它的性能优于`std::sort`,而且我相信一个多位数位化器还能更快,但我的测试套件遇到了内存问题(不是算法本身,而是测试套件本身)。 - Mooing Duck
@KonradRudolph:Flag排序和其他基数排序之间的重大区别在于计数传递。你说得对,所有基数排序都非常相似,但我不认为你的是Flag排序。 - Mooing Duck
@MooingDuck:刚从你的示例中获得了一些灵感——我在自己的独立实现中卡住了,而你的帮助我重新回到了正轨。谢谢!一个可能的优化——我还没有进展到足够的程度来确定它是否值得:如果您要交换到的位置上的元素恰好已经在它应该在的位置上,您可能希望跳过它并前进到另一个位置。当然,检测这将需要额外的逻辑和可能的额外存储,但由于相对于比较而言,交换是昂贵的,因此这样做可能是值得的。 - 500 - Internal Server Error

2
Radix-Sort不是缓存感知的,也不是用于大型数据集的最快排序算法。 你可以看一下: 您还可以在存储到排序数组之前对DNA的每个字母进行压缩,并将其编码为2位。

Bill:你能解释一下qsort函数相对于C++提供的std::sort函数有哪些优势吗?特别是后者在现代库中实现了高度复杂的introsort,并内联比较操作。我不相信它在大多数情况下都能以O(n)的速度运行,因为这需要一定程度的内省,而这在一般情况下是不可用的(至少没有很多开销)。 - Konrad Rudolph
我并不使用C++,但在我的测试中,内联QSORT比stdlib中的qsort快3倍。对于整数来说,ti7qsort是最快的排序算法(比内联QSORT还要快)。你也可以用它来排序小型固定大小的数据。你必须使用你自己的数据进行测试。 - bill

1

首先,考虑你的问题的编码。摆脱字符串,用二进制表示替换它们。使用第一个字节来指示长度+编码。或者,在四字节边界上使用固定长度表示。然后基数排序就变得容易得多了。对于基数排序,最重要的是不要在内部循环的热点处进行异常处理。

好的,我再想一想4元问题。你需要像Judy tree这样的解决方案。下一个解决方案可以处理可变长度的字符串;对于固定长度,只需删除长度位,实际上这使得它更容易。

分配16个指针的块。指针的最低有效位可以被重用,因为你的块总是对齐的。你可能需要一个特殊的存储分配器(将大的存储分成小块)。有许多不同类型的块:

  • 使用7个长度位的可变长度字符串进行编码。当它们填满时,你可以用以下方式替换它们:
  • 位置编码下两个字符,你有16个指向下一个块的指针,以此类推:
  • 位图编码字符串的最后三个字符。

对于每种块,您需要在LSB中存储不同的信息。由于您有可变长度的字符串,因此需要存储字符串结束符,并且最后一种块仅适用于最长的字符串。随着深入结构,7个长度位应该被替换为更少的位。

这为您提供了一种相当快速和非常内存高效的排序字符串存储方式。它的行为有点像trie。要使其正常工作,请确保构建足够的单元测试。您需要覆盖所有块转换。您需要从第二种块开始。

为了获得更高的性能,您可能需要添加不同的块类型和更大的块大小。如果块始终具有相同的大小并且足够大,则可以使用更少的指针位。使用16个指针的块大小,在32位地址空间中已经有一个字节可用。查看Judy树文档以获取有趣的块类型。基本上,您需要增加代码和工程时间来进行空间(和运行时)折衷。

你可能想要从第一个四个字符开始使用256宽的直接基数排序。这提供了一个不错的空间/时间权衡。在这个实现中,与简单的trie相比,你会得到更少的内存开销;它大约小了三倍(我没有测量)。如果常数足够低,O(n)不是问题,就像你在与O(n log n)快速排序进行比较时注意到的那样。

你有兴趣处理双精度吗?对于短序列,将会有一些。调整块以处理计数是棘手的,但可以非常节省空间。


如果我使用位压缩表示,我不明白基数排序在我的情况下如何变得更容易。顺便说一下,我使用的框架实际上提供了使用位压缩表示的可能性,但对于我作为接口用户来说,这是完全透明的。 - Konrad Rudolph
不过当你看着秒表的时候就不一样了 :) - Stephan Eggermont
我一定会研究 Judy 树。然而,普通的香草树并没有什么特别之处,因为它们基本上像普通的 MSD 基数排序一样运行,只是在元素的遍历过程中需要额外的存储空间。 - Konrad Rudolph

1

dsimcha的MSB基数排序看起来很不错,但Nils通过观察到缓存局部性是在大问题规模下导致效率低下的根本原因而更接近问题的核心。

我建议采用非常简单的方法:

  1. 经验估计出最大的大小m,使得基数排序高效。
  2. 一次读取m个元素的块,对它们进行基数排序,并将它们写出(如果有足够的内存,则写入内存缓冲区,否则写入文件),直到耗尽输入。
  3. 归并排序结果已排序的块。

归并排序是我所知道的最友好缓存的排序算法:“从数组A或B中读取下一个项,然后将一个项写入输出缓冲区。” 它在磁带驱动器上运行效率高。 它确实需要2n的空间来排序n个项目,但我打赌你会看到大大改善的缓存局部性,这将使其变得不重要 - 如果您使用非就地基数排序,则无论如何都需要额外的空间。

最后请注意,归并排序可以在没有递归的情况下实现,事实上这样做可以清楚地显示真正的线性内存访问模式。


0

虽然接受的答案完美地解答了这个问题的描述,但是我在这里徒劳地寻找一个将数组内联分割为N部分的算法。我自己写了一个,所以在这里分享一下。

警告:这不是一个稳定的分区算法,因此对于多级分区,必须重新对每个结果分区进行分区,而不是整个数组。优点是它是内联的。

它对所提出的问题有帮助的方式是,您可以根据字符串的某个字母重复地进行内联分区,然后在它们足够小的时候使用您选择的算法对分区进行排序。

  function partitionInPlace(input, partitionFunction, numPartitions, startIndex=0, endIndex=-1) {
    if (endIndex===-1) endIndex=input.length;
    const starts = Array.from({ length: numPartitions + 1 }, () => 0);
    for (let i = startIndex; i < endIndex; i++) {
      const val = input[i];
      const partByte = partitionFunction(val);
      starts[partByte]++;
    }
    let prev = startIndex;
    for (let i = 0; i < numPartitions; i++) {
      const p = prev;
      prev += starts[i];
      starts[i] = p;
    }
    const indexes = [...starts];
    starts[numPartitions] = prev;
  
    let bucket = 0;
    while (bucket < numPartitions) {
      const start = starts[bucket];
      const end = starts[bucket + 1];
      if (end - start < 1) {
        bucket++;
        continue;
      }
      let index = indexes[bucket];
      if (index === end) {
        bucket++;
        continue;
      }
  
      let val = input[index];
      let destBucket = partitionFunction(val);
      if (destBucket === bucket) {
        indexes[bucket] = index + 1;
        continue;
      }
  
      let dest;
      do {
        dest = indexes[destBucket] - 1;
        let destVal;
        let destValBucket = destBucket;
        while (destValBucket === destBucket) {
          dest++;
          destVal = input[dest];
          destValBucket = partitionFunction(destVal);
        }
  
        input[dest] = val;
        indexes[destBucket] = dest + 1;
  
        val = destVal;
        destBucket = destValBucket;
      } while (dest !== index)
    }
    return starts;
  }

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接