最快的代码C / C ++来select一组27个浮点值的中位数
这是众所周知的selectalgorithm。 请参阅http://en.wikipedia.org/wiki/Selection_algorithm 。
我需要它来find一组3×3×3体素值的中值。 由于体积是由十亿个体素组成，并且algorithm是recursion的，所以最好稍微快一点。 总体而言，可以预期的是价值相对接近。
迄今为止我已经试过的最快速的已知algorithm使用快速sorting分区function。 我想知道是否有更快的。
我已经“发明”了使用两堆的速度提高了20％，但是使用散列的速度更快。 在执行这个之前，我想知道是否已经有一个快速的解决scheme。
我使用浮点数的事实应该是无关紧要的，因为它们可以在反转符号位之后被认为是无符号整数。 订单将被保留。
编辑：基准和源代码转移到一个单独的答案，由戴维兰德曼build议。 请看下面的chmike的答案。
编辑 ：到目前为止，最有效的algorithm是由Boojum引用作为快速中值和双边过滤纸，现在是这个问题的答案的链接。 这种方法的第一个聪明的想法是使用基数sorting，其次是结合中间search相邻像素谁共享很多像素。
由于这听起来像是在大量卷数据上执行中值滤波，所以您可能需要查看SIGGRAPH 2006年的快速中值滤波和双边滤波 。本文讨论二维image processing，但您可能会能够适应3D体积的algorithm。 如果没有别的，它可能会给你一些想法，如何退后一步，从稍微不同的angular度来看问题。
selectalgorithm是线性时间（O（n））。 复杂性方面，你不可能比线性时间做得更好，因为需要线性时间来读取所有数据。 所以你不可能有更复杂的东西。 也许你在某些投入上有更快的速度？ 我怀疑这会有很大的不同。
C ++已经包含线性时间selectalgorithm。 为什么不使用它？
std::vector<YourType>::iterator first = yourContainer.begin(); std::vector<YourType>::iterator last = yourContainer.end(); std::vector<YourType>::iterator middle = first + (last - first) / 2; std::nth_element(first, middle, last); // can specify comparator as optional 4th arg YourType median = *middle; 编辑：从技术上讲，这只是一个奇数长度容器的中位数。 对于一个甚至长度，它将获得“上”中位数。 如果你想要中等长度的传统定义，你可能需要运行两次，每次在first + (last - first) / 2和first + (last - first) / 2 - 1然后平均他们什么的。
编辑：我必须道歉。 下面的代码是错误的。 我有固定的代码，但需要find一个icc编译器来重做测量。
到目前为止所考虑的algorithm的基准结果
有关algorithm的协议和简短说明，请参见下文。 第一个值是超过200个不同序列的平均时间（秒），第二个值是stdDev。
HeapSort : 2.287 0.2097 QuickSort : 2.297 0.2713 QuickMedian1 : 0.967 0.3487 HeapMedian1 : 0.858 0.0908 NthElement : 0.616 0.1866 QuickMedian2 : 1.178 0.4067 HeapMedian2 : 0.597 0.1050 HeapMedian3 : 0.015 0.0049 <-- best
协议：使用从rand（）获得的随机位生成27个随机数。 连续500万次应用每个algorithm（包括之前的数组拷贝），并计算200个随机序列的平均值和stdDev。 用icc -S -O3编译的C ++代码，并运行在带有8GB DDR3的Intel E8400上。
algorithm：
堆sorting：使用堆sorting和select中间值的完整sorting。 天真的实现使用下标访问。
快速sorting：使用快速sorting并选取中间值，完成适当的sorting。 天真的实现使用下标访问。
QuickMedian1：快速selectalgorithm与交换。 天真的实现使用下标访问。
HeapMedian1：采用先前交换的平衡堆方法。 天真的实现使用下标访问。
NthElement：使用nth_element STLalgorithm。 使用memcpy（vct.data（），rndVal，…）将数据复制到vector中;
QuickMedian2：使用指针的快速selectalgorithm并复制到两个缓冲区中以避免交换。 基于MSalters的build议。
HeapMedian2：我发明的algorithm使用双堆共享头的变体。 左堆具有最大的值作为头，右有最小的值作为头。 初始化第一个值作为公共头和第一个中值猜测。 如果小于head，则将后续值添加到左堆中，否则将其添加到右堆中，直到其中一个堆已满。 它包含14个值时已满。 然后只考虑整个堆。 如果它的权利堆，所有值大于头，popup头和插入值。 忽略所有其他值。 如果它的左堆，对于所有小于头的值，popup头并将其插入堆中。 忽略所有其他值。 当所有的数值已经进行时，共同的头是中间值。 它使用整数索引到数组中。 使用指针（64位）的版本似乎慢了近一倍（〜1s）。
HeapMedian3：与HeapMedian2相同的algorithm，但优化。 它使用无符号字符索引，避免了价值交换和其他各种小事情。 平均值和stdDev值是在1000个随机序列上计算的。 对于nth_element，我用相同的1000个随机序列测量了0.508s和一个0.159537的stdDev。 HeapMedian3因此比nth_element stl函数快33倍。 检查每个返回的中值是否与heapSort返回的中值相匹配。 我怀疑使用散列的方法可能会明显更快。
编辑1：这个algorithm可以进一步优化。 根据比较结果将元素分派到左侧或右侧堆的第一阶段不需要堆积。 简单地将元素附加到两个无序序列就足够了。 一旦一个序列满了，第一阶段就会停止，这意味着它包含了14个元素（包括中值）。 第二阶段首先对整个序列进行堆积处理，然后按照HeapMedian3algorithm进行处理。 我会尽快提供新的代码和基准。
编辑2：我实施和基准优化algorithm。 但是heapMedian3没有显着的性能差异。 它的平均速度甚至更慢。 显示的结果被确认。 可能会有更大的集合。 还要注意，我只是select第一个值作为初始中位数猜测。 如所暗示的那样，可以从我们在“重叠”值集合中search中值的事实中受益。 使用中值algorithm的中位数将有助于select一个更好的初始中值猜测。
HeapMedian3的源代码
// return the median value in a vector of 27 floats pointed to by a float heapMedian3( float *a ) { float left[14], right[14], median, *p; unsigned char nLeft, nRight; // pick first value as median candidate p = a; median = *p++; nLeft = nRight = 1; for(;;) { // get next value float val = *p++; // if value is smaller than median, append to left heap if( val < median ) { // move biggest value to the heap top unsigned char child = nLeft++, parent = (child - 1) / 2; while( parent && val > left[parent] ) { left[child] = left[parent]; child = parent; parent = (parent - 1) / 2; } left[child] = val; // if left heap is full if( nLeft == 14 ) { // for each remaining value for( unsigned char nVal = 27 - (p - a); nVal; --nVal ) { // get next value val = *p++; // if value is to be inserted in the left heap if( val < median ) { child = left[2] > left[1] ? 2 : 1; if( val >= left[child] ) median = val; else { median = left[child]; parent = child; child = parent*2 + 1; while( child < 14 ) { if( child < 13 && left[child+1] > left[child] ) ++child; if( val >= left[child] ) break; left[parent] = left[child]; parent = child; child = parent*2 + 1; } left[parent] = val; } } } return median; } } // else append to right heap else { // move smallest value to the heap top unsigned char child = nRight++, parent = (child - 1) / 2; while( parent && val < right[parent] ) { right[child] = right[parent]; child = parent; parent = (parent - 1) / 2; } right[child] = val; // if right heap is full if( nRight == 14 ) { // for each remaining value for( unsigned char nVal = 27 - (p - a); nVal; --nVal ) { // get next value val = *p++; // if value is to be inserted in the right heap if( val > median ) { child = right[2] < right[1] ? 2 : 1; if( val <= right[child] ) median = val; else { median = right[child]; parent = child; child = parent*2 + 1; while( child < 14 ) { if( child < 13 && right[child+1] < right[child] ) ++child; if( val <= right[child] ) break; right[parent] = right[child]; parent = child; child = parent*2 + 1; } right[parent] = val; } } } return median; } } } }
这个问题不能简单地回答，原因很简单，一个algorithm相对于另一个algorithm的性能取决于编译器/处理器/数据结构的组合，就像algorithm本身一样，你肯定知道
所以你的方法尝试一下就好了。 是的，快速sorting应该是相当快的。 如果你还没有这样做，你可能会想尝试insertionsort，它往往在小数据集上performance更好。 这就是说，只要做一个sortingalgorithm就足够快。 select“正确的”algorithm通常不会快10倍。
为了获得实质性的加速，更好的方法是使用更多的结构。 一些过去为我工作的想法带来了大规模的问题：
-
你可以有效地预先计算，而创build体素和存储28而不是27浮点数？
-
大致的解决scheme是否足够好？ 如果是这样的话，只要看一下9个数值的中位数就可以了，因为“总的来说，可以预期数值相对接近”。 或者只要值相对接近，就可以用平均值代替它。
-
你真的需要所有数十亿体素的中位数吗？ 也许你有一个简单的testing，你是否需要中位数，然后才能计算相关的子集。
-
如果没有其他帮助：查看编译器生成的asm代码。 你也许可以写更快的asm代码（例如通过使用寄存器来完成所有的计算）。
编辑：为了什么是值得的，我已经附上下面的评论（完全未经testing）中提到的（部分）insertionsort代码。 如果numbers[]是一个大小为N的数组，并且您希望在数组开头处sorting的最小P浮点数，请调用partial_insertionsort<N, P, float>(numbers); 。 因此，如果你调用partial_insertionsort<27, 13, float>(numbers); ， numbers[13]将包含中位数。 为了获得更多的速度，你也必须展开while循环。 正如上面所讨论的，为了获得真正的速度，你必须使用你对数据的知识（例如，数据已经部分sorting了吗？你知道数据分布的属性吗？我想，你会得到漂移）。
template <long i> class Tag{}; template<long i, long N, long P, typename T> inline void partial_insertionsort_for(T a[], Tag<N>, Tag<i>) { long j = i <= P+1 ? i : P+1; // partial sort T temp = a[i]; a[i] = a[j]; // compiler should optimize this away where possible while(temp < a[j - 1] && j > 0) { a[j] = a[j - 1]; j--;} a[j] = temp; partial_insertionsort_for<i+1,N,P,T>(a,Tag<N>(),Tag<i+1>());} template<long i, long N, long P, typename T> inline void partial_insertionsort_for(T a[], Tag<N>, Tag<N>){} template <long N, long P, typename T> inline void partial_insertionsort(T a[]) {partial_insertionsort_for<0,N,P,T>(a, Tag<N>(), Tag<0>());}
在你的第一次尝试中使用的最可能的algorithm只是nth_element; 它几乎直接给你你想要的东西。 只是要求第十四个要素。
第二次尝试时，目标是利用固定的数据大小。 根本不用分配任何内存。 因此，将您的体素值复制到27个元素的预分配数组中。 select一个枢纽，并将其复制到53元素数组的中间。 将其余值复制到数据透视的任一侧。 这里你保留两个指针（ float* left = base+25, *right=base+27 ）。 现在有三种可能性：左侧更大，右侧更大，或者两者都有12个元素。 最后一个例子是微不足道的。 你的支点是中位数。 否则，请在左侧或右侧调用nth_element。 N的确切值取决于有多less值大于或小于主轴。 例如，如果除法是12/14，那么你需要最小的元素大于主元素，所以Nth = 0，如果除法是14/12，则需要最小的元素小于主元素，所以Nth = 13。 最糟糕的情况是26/0和0/26，当你的枢轴是一个极端的，但那些只发生在所有情况的二十七分之二。
第三个改进（或者第一个，如果你必须使用C并且没有nth_element）完全replacenth_element。 你仍然有53个元素的数组，但这次你直接从体素值中填充它（把你的临时副本保存到一个float[27] ）。 在这个第一次迭代的主轴只是voxel [0] [0] [0]。 对于后续的迭代，你使用第二个预先分配的float[53] （如果两者的大小相同，则更容易），并在两者之间复制浮点数。 这里的基本迭代步骤仍然是：将数据透视图复制到中间，将剩余的数据向左和向右sorting。 在每个步骤结束时，您将知道中位数是否小于或大于当前枢轴，因此您可以丢弃大于或小于该枢轴的浮点数。 每次迭代，这消除了1到12个元素，平均剩下的25％。
最后一个迭代，如果你还需要更多的速度，是基于观察你的体素大部分重叠显着。 您预先计算每个3x3x1切片的中位数值。 然后，当你需要一个3x3x3体素立方体的初始支点时，你需要三者的中位数。 您先前知道有9个体素较小，9个体素比中位数（4 + 4 + 1）大9个体素。 所以，在第一个转折步骤之后，最坏的情况是9/17和17/9分割。 所以，你只需要在浮点[17]中find第4或第13个元素，而不是在浮点[12]中的第12或第14个元素。
背景：使用左指针和右指针首先复制一个枢轴，然后将float [N]的其余部分复制到float [2N-1]的想法是，您将围绕枢轴填充一个float [N]子arrays，所有元素小于左侧（较低的指数）和较高的指数（较高的指数）。 现在，如果你想要第M个元素，你可能会发现自己很幸运，并且有M-1元素小于枢轴，在这种情况下，枢轴是你需要的元素。 如果有多于（M-1）个元素小于枢轴，第M个元素就在它们之间，所以你可以丢弃枢轴和大于枢轴的任何东西，而对于第M个元素，可以舍弃所有较低值。 如果小于 （M-1）的元素小于主元，则您正在查找比主元更高的值。 所以，你会放弃这个枢纽和比它小的东西。 让元素的数量小于枢轴，即枢轴的左边为L.在下一次迭代中，要（NL-1）个浮点数大于枢轴的第（ML-1）个元素。
这种nth_elementalgorithm是相当有效的，因为大部分的工作都是在两个小数组之间复制浮点数，而这两个小数组都将在caching中，而且由于您的状态大部分时间都是由3个指针（源指针，左目的指针，正确的目标指针）。
显示基本代码：
float in[27], out[53]; float pivot = out[26] = in[0]; // pivot float* left = out+25, right = out+27 for(int i = 1; i != 27; ++1) if((in[i]<pivot)) *left-- = in[i] else *right++ = in[i]; // Post-condition: The range (left+1, right) is initialized. // There are 25-(left-out) floats <pivot and (right-out)-27 floats >pivot
我想你最好的select是采取一个现有的sortingalgorithm，并试图找出是否可以适应它，使该集不需要完全sorting。 为了确定中位数，您最多需要sorting的值的一半，或者更低或更高的一半就足够了：
original: | 5 | 1 | 9 | 3 | 3 | sorted: | 1 | 3 | 3 | 5 | 9 | lower half sorted: | 1 | 3 | 3 | 9 | 5 | higher half sorted: | 3 | 1 | 3 | 5 | 9 |
另一半是一堆未分类的价值，只是分享大/小或等于最大/最小sorting值的属性。
但是，我还没有准备好的algorithm，这只是一个想法，你可能会在你的sorting中采取捷径。
使用Bose-Nelsonalgorithm生成的sortingnetworking将使用173个比较直接find没有循环/recursion的中值。 如果您可以并行执行比较（比如使用向量算术指令），那么您可以将比较分组为less至28个并行操作。
如果你确定浮点数已经归一化，而不是（qs）NaN，那么你可以使用整数运算来比较IEEE-754浮点数，这些浮点数可以在某些CPU上更好地执行。
将这个sortingnetworking直接转换为C（gcc 4.2）会在Core i7上产生388个时钟周期的最坏情况。
sortingnetworking
亚历克斯·斯捷潘诺夫（Alex Stepanov）的新书“编程的元素 ”（ Elements of Programming）在一定程度上讲述了使用平均比较的最小数量来查找顺序统计量，同时最小化运行时开销 不幸的是，需要大量的代码来计算5个元素的中位数，即使这样，他作为一个项目寻找一个替代解决scheme，平均使用比较less的比较，所以我不会梦想延长框架find27个元素的中位数。 而这本书甚至在2009年6月15日之前都不可用。重点在于，因为这是一个固定大小的问题，所以有一个直接的比较方法是可行的。
而且，这个algorithm不是孤立地运行一次，而是多次运行，在大多数运行之间，只有9个值会改变。 这意味着理论上已经有一些工作已经完成了。 但是，我还没有听说过任何利用这个事实的image processing中的中值滤波algorithm。
对于每个提到了nth_element的人来说都是+1，但是这种代码是手写algorithm比STL更好的地方，因为你想要为在特定数据集的CPU上运行的那个编译器生成最高效的代码。 例如，对于一些CPU /编译器组合std :: swap（int，int），可能比使用XOR手写交换要慢（在回复之前，我知道这可能是20年前的事情，但现在不再）。 有时性能是通过手工编写特定于您的CPU的汇编代码而获得的。 如果您打算利用GPU的stream处理器，则可能需要相应地devisealgorithm。
你提到使用2堆，并跟踪插入的中位数。 这就是我刚才在一个项目中所做的。 我改变了arrays，只使用了一个堆。 我想不出任何更快的algorithm，但我想告诉你有关内存使用情况，特别是CPUcaching内存。 你想要小心访问内存。 CPU高速caching是按页面交换的，所以你希望你的algorithm能够触碰紧密靠近的内存，以尽量减lessCPUcaching未命中。
当我们说一百万个不同的值时，你需要中位数。 是否有可能把你的中位数作为百万分之一，比如10％。 因此，中位数接近第二个元素，它将两个相等（或几乎相等）的子集中的值相除。 因此，为了find中位数，你需要less于O（n）次（在这个例子中是O（1 / 10n）），因此在O（nlogn）中用快速sorting来接近最优sorting？
如果你想看Donald E. Knuth所着的algorithm，
PS。 如果您认为自己已经发明了更好的东西，那么您应该能够certificate复杂性与已知algorithm的复杂性相似或更好。 另一方面，基于桶和基数的变化是O（n），快速sorting只是O（n.log（n））。 一个快20％的方法仍然是O（n.log（n）），直到你可以显示algorithm:-)
我敢打赌，你可以计算它们的零成本 – 在从磁盘加载的单独的线程中（或者它们被生成）。
我真正说的是'速度'不是来自一点点的混乱，因为27个值不足以使大O符号成为一个真正的因素。
你可能想看Knuth的练习5.3.3.13。 它描述了一个由Floyd计算的algorithm，它使用（3/2）n + O（n ^（2/3）log n）比较来找出n个元素的中位数，并且隐藏在O（·）中的常量似乎不是在实践中太大了。
如果有3x3x3 = 27个可能的值（如果是这样的话为什么浮点数？），你能创build一个由27个元素组成的数组，并且可以一次遍历数据来计算每个可能性吗？
我的一个一维数据集的中位数计算的超级快速algorithm在三遍中完成了工作，并且不需要对数据集进行sorting（!!!）。
一个非常通用的描述如下：
- 通过1：扫描一维数据集并收集数据集的一些统计信息
- 通过2：使用数据集的统计信息并应用一些数据挖掘来创build中间（辅助）数组
- 通过3：扫描中间（帮助者）arrays，以find中位数
该algorithmdevise用于查找单精度浮点值（在具有32GB物理内存和128GB虚拟内存的桌面系统上）超过8GE（千兆元素）的超大1-D数据集的中值，或用于查找中位数在硬实时环境中的小数据集。
该algorithm是：
- 比传统的基于Heap或Mergesortingalgorithm的algorithm快60〜75倍
- 以纯C语言实现
- 不使用任何英特尔内部函数
- 不使用任何内联汇编程序指令
- 绝对可移植的C / C ++编译器之间，如MS，Intel，MinGW，Borland，Turbo和Watcom
- 平台之间绝对可移植
最好的问候，谢尔盖Kostrov
