极其频繁地调用 std::nth_element() 函数

Question

我在任何地方都没有找到这个特定的主题...

我在 23 个整数的 std::vector 中调用 nth_element() 算法大约 400,000 次，更精确的“无符号短”值。

我想提高计算速度，而这个特定的调用需要大量的 CPU 时间。 现在我注意到，与 std::sort() 一样，即使在最高优化级别和 NDEBUG 模式（Linux Clang 编译器）下，nth_element 函数在分析器中也是可见的，因此比较是内联的，而不是函数调用本身。好吧，更确切地说：不是 nth_element() 而是 std::__introselect() 是可见的。

由于数据量小，我尝试使用二次排序函数PIKSORT，当数据量小于20个元素时，往往比调用std::sort要快，可能是因为函数会内联.

template <class CONTAINER>
inline void piksort(CONTAINER& arr)  // indeed this is "insertion sort"
{
    typename CONTAINER::value_type a;

    const int n = (int)arr.size();
    for (int j = 1; j<n; ++j) {
        a = arr[j];
        int i = j;
        while (i > 0 && a < arr[i - 1]) {
            arr[i] = arr[i - 1];
            i--;
        }
        arr[i] = a;
    }
}

但是在这种情况下，这比使用 nth_element 慢。

另外，使用统计方法也不合适，Something faster than std::nth_element

最后，由于值在 0 到大约 20000 的范围内，因此直方图方法看起来不合适。

我的问题：有人知道一个简单的解决方案吗？我想我可能不是唯一一个必须经常调用 std::sort 或 nth_element 的人。

Answer 1

您提到数组的大小始终为23。此外，使用的类型是unsigned short。在这种情况下，您可以尝试使用大小为23 的排序网络；由于您的类型是unsigned short，因此使用排序网络对整个数组进行排序可能比使用std::nth_element 部分排序更快。这是大小为23 和118 比较交换单元的排序网络的一个非常简单的C++14 实现，如Using Symmetry and Evolutionary Search to Minimize Sorting Networks 所述：

template<typename RandomIt, typename Compare = std::less<>>
void network_sort23(RandomIt first, Compare compare={})
{
    swap_if(first[1u], first[20u], compare);
    swap_if(first[2u], first[21u], compare);
    swap_if(first[5u], first[13u], compare);
    swap_if(first[9u], first[17u], compare);
    swap_if(first[0u], first[7u], compare);
    swap_if(first[15u], first[22u], compare);
    swap_if(first[4u], first[11u], compare);
    swap_if(first[6u], first[12u], compare);
    swap_if(first[10u], first[16u], compare);
    swap_if(first[8u], first[18u], compare);
    swap_if(first[14u], first[19u], compare);
    swap_if(first[3u], first[8u], compare);
    swap_if(first[4u], first[14u], compare);
    swap_if(first[11u], first[18u], compare);
    swap_if(first[2u], first[6u], compare);
    swap_if(first[16u], first[20u], compare);
    swap_if(first[0u], first[9u], compare);
    swap_if(first[13u], first[22u], compare);
    swap_if(first[5u], first[15u], compare);
    swap_if(first[7u], first[17u], compare);
    swap_if(first[1u], first[10u], compare);
    swap_if(first[12u], first[21u], compare);
    swap_if(first[8u], first[19u], compare);
    swap_if(first[17u], first[22u], compare);
    swap_if(first[0u], first[5u], compare);
    swap_if(first[20u], first[21u], compare);
    swap_if(first[1u], first[2u], compare);
    swap_if(first[18u], first[19u], compare);
    swap_if(first[3u], first[4u], compare);
    swap_if(first[21u], first[22u], compare);
    swap_if(first[0u], first[1u], compare);
    swap_if(first[19u], first[22u], compare);
    swap_if(first[0u], first[3u], compare);
    swap_if(first[12u], first[13u], compare);
    swap_if(first[9u], first[10u], compare);
    swap_if(first[6u], first[15u], compare);
    swap_if(first[7u], first[16u], compare);
    swap_if(first[8u], first[11u], compare);
    swap_if(first[11u], first[14u], compare);
    swap_if(first[4u], first[11u], compare);
    swap_if(first[6u], first[8u], compare);
    swap_if(first[14u], first[16u], compare);
    swap_if(first[17u], first[20u], compare);
    swap_if(first[2u], first[5u], compare);
    swap_if(first[9u], first[12u], compare);
    swap_if(first[10u], first[13u], compare);
    swap_if(first[15u], first[18u], compare);
    swap_if(first[10u], first[11u], compare);
    swap_if(first[4u], first[7u], compare);
    swap_if(first[20u], first[21u], compare);
    swap_if(first[1u], first[2u], compare);
    swap_if(first[7u], first[15u], compare);
    swap_if(first[3u], first[9u], compare);
    swap_if(first[13u], first[19u], compare);
    swap_if(first[16u], first[18u], compare);
    swap_if(first[8u], first[14u], compare);
    swap_if(first[4u], first[6u], compare);
    swap_if(first[18u], first[21u], compare);
    swap_if(first[1u], first[4u], compare);
    swap_if(first[19u], first[21u], compare);
    swap_if(first[1u], first[3u], compare);
    swap_if(first[9u], first[10u], compare);
    swap_if(first[11u], first[13u], compare);
    swap_if(first[2u], first[6u], compare);
    swap_if(first[16u], first[20u], compare);
    swap_if(first[4u], first[9u], compare);
    swap_if(first[13u], first[18u], compare);
    swap_if(first[19u], first[20u], compare);
    swap_if(first[2u], first[3u], compare);
    swap_if(first[18u], first[20u], compare);
    swap_if(first[2u], first[4u], compare);
    swap_if(first[5u], first[17u], compare);
    swap_if(first[12u], first[14u], compare);
    swap_if(first[8u], first[12u], compare);
    swap_if(first[5u], first[7u], compare);
    swap_if(first[15u], first[17u], compare);
    swap_if(first[5u], first[8u], compare);
    swap_if(first[14u], first[17u], compare);
    swap_if(first[3u], first[5u], compare);
    swap_if(first[17u], first[19u], compare);
    swap_if(first[3u], first[4u], compare);
    swap_if(first[18u], first[19u], compare);
    swap_if(first[6u], first[10u], compare);
    swap_if(first[11u], first[16u], compare);
    swap_if(first[13u], first[16u], compare);
    swap_if(first[6u], first[9u], compare);
    swap_if(first[16u], first[17u], compare);
    swap_if(first[5u], first[6u], compare);
    swap_if(first[4u], first[5u], compare);
    swap_if(first[7u], first[9u], compare);
    swap_if(first[17u], first[18u], compare);
    swap_if(first[12u], first[15u], compare);
    swap_if(first[14u], first[15u], compare);
    swap_if(first[8u], first[12u], compare);
    swap_if(first[7u], first[8u], compare);
    swap_if(first[13u], first[15u], compare);
    swap_if(first[15u], first[17u], compare);
    swap_if(first[5u], first[7u], compare);
    swap_if(first[9u], first[10u], compare);
    swap_if(first[10u], first[14u], compare);
    swap_if(first[6u], first[11u], compare);
    swap_if(first[14u], first[16u], compare);
    swap_if(first[15u], first[16u], compare);
    swap_if(first[6u], first[7u], compare);
    swap_if(first[10u], first[11u], compare);
    swap_if(first[9u], first[12u], compare);
    swap_if(first[11u], first[13u], compare);
    swap_if(first[13u], first[14u], compare);
    swap_if(first[8u], first[9u], compare);
    swap_if(first[7u], first[8u], compare);
    swap_if(first[14u], first[15u], compare);
    swap_if(first[9u], first[10u], compare);
    swap_if(first[8u], first[9u], compare);
    swap_if(first[12u], first[14u], compare);
    swap_if(first[11u], first[12u], compare);
    swap_if(first[12u], first[13u], compare);
    swap_if(first[10u], first[11u], compare);
    swap_if(first[11u], first[12u], compare);
}

swap_if 实用函数将x 和y 两个参数与谓词compare 进行比较，如果compare(y, x) 则交换它们。我的示例使用了一个通用的swap_if 函数，但是如果您知道无论如何您都会将unsigned short 值与operator< 进行比较，则可以使用优化版本（如果您的编译器识别并优化比较，您可能不需要这样的函数-exchange，但不幸的是，并非所有编译器都这样做 - 我正在使用带有 -O3 的 g++5.2，我仍然需要以下函数来提高性能）：

void swap_if(unsigned short& x, unsigned short& y)
{
    unsigned short dx = x;
    unsigned short dy = y;
    unsigned short tmp = x = std::min(dx, dy);
    y ^= dx ^ tmp;
}

现在，为了确保它确实更快，我决定在需要时对 std::nth_element 计时，仅对前 10 个元素进行部分排序，而不是使用排序网络对整个 23 个元素进行排序（1000000 次不同的洗牌数组）。这是我得到的：

std::nth_element    1158ms
network_sort23      487ms

也就是说，我的电脑已经运行了一段时间，速度有点慢，但性能差异很明显。我相信当我重新启动计算机时，这种差异将保持不变。稍后我可能会尝试并通知您。

关于这些时间是如何产生的，我使用了 this benchmark 的修改版本，来自我的 cpp-sort library。原始排序网络和swap_if 函数也来自那里，因此您可以确定它们已经过多次测试:)

编辑：这是我重新启动计算机后的结果。 network_sort23 版本仍然比 std::nth_element 快两倍：

std::nth_element    369ms
network_sort23      154ms

EDIT²：如果您只需要中位数，您可以简单地删除计算将位于第 11 位的最终值不需要的比较交换单元。生成的大小为 23 的中值查找网络使用与前一个不同的大小为 23 的排序网络，它产生的结果略好：

swap_if(first[0u], first[1u], compare);
swap_if(first[2u], first[3u], compare);
swap_if(first[4u], first[5u], compare);
swap_if(first[6u], first[7u], compare);
swap_if(first[8u], first[9u], compare);
swap_if(first[10u], first[11u], compare);
swap_if(first[1u], first[3u], compare);
swap_if(first[5u], first[7u], compare);
swap_if(first[9u], first[11u], compare);
swap_if(first[0u], first[2u], compare);
swap_if(first[4u], first[6u], compare);
swap_if(first[8u], first[10u], compare);
swap_if(first[1u], first[2u], compare);
swap_if(first[5u], first[6u], compare);
swap_if(first[9u], first[10u], compare);
swap_if(first[1u], first[5u], compare);
swap_if(first[6u], first[10u], compare);
swap_if(first[5u], first[9u], compare);
swap_if(first[2u], first[6u], compare);
swap_if(first[1u], first[5u], compare);
swap_if(first[6u], first[10u], compare);
swap_if(first[0u], first[4u], compare);
swap_if(first[7u], first[11u], compare);
swap_if(first[3u], first[7u], compare);
swap_if(first[4u], first[8u], compare);
swap_if(first[0u], first[4u], compare);
swap_if(first[7u], first[11u], compare);
swap_if(first[1u], first[4u], compare);
swap_if(first[7u], first[10u], compare);
swap_if(first[3u], first[8u], compare);
swap_if(first[2u], first[3u], compare);
swap_if(first[8u], first[9u], compare);
swap_if(first[2u], first[4u], compare);
swap_if(first[7u], first[9u], compare);
swap_if(first[3u], first[5u], compare);
swap_if(first[6u], first[8u], compare);
swap_if(first[3u], first[4u], compare);
swap_if(first[5u], first[6u], compare);
swap_if(first[7u], first[8u], compare);
swap_if(first[12u], first[13u], compare);
swap_if(first[14u], first[15u], compare);
swap_if(first[16u], first[17u], compare);
swap_if(first[18u], first[19u], compare);
swap_if(first[20u], first[21u], compare);
swap_if(first[13u], first[15u], compare);
swap_if(first[17u], first[19u], compare);
swap_if(first[12u], first[14u], compare);
swap_if(first[16u], first[18u], compare);
swap_if(first[20u], first[22u], compare);
swap_if(first[13u], first[14u], compare);
swap_if(first[17u], first[18u], compare);
swap_if(first[21u], first[22u], compare);
swap_if(first[13u], first[17u], compare);
swap_if(first[18u], first[22u], compare);
swap_if(first[17u], first[21u], compare);
swap_if(first[14u], first[18u], compare);
swap_if(first[13u], first[17u], compare);
swap_if(first[18u], first[22u], compare);
swap_if(first[12u], first[16u], compare);
swap_if(first[15u], first[19u], compare);
swap_if(first[16u], first[20u], compare);
swap_if(first[12u], first[16u], compare);
swap_if(first[13u], first[16u], compare);
swap_if(first[19u], first[22u], compare);
swap_if(first[15u], first[20u], compare);
swap_if(first[14u], first[15u], compare);
swap_if(first[20u], first[21u], compare);
swap_if(first[14u], first[16u], compare);
swap_if(first[19u], first[21u], compare);
swap_if(first[15u], first[17u], compare);
swap_if(first[18u], first[20u], compare);
swap_if(first[15u], first[16u], compare);
swap_if(first[17u], first[18u], compare);
swap_if(first[19u], first[20u], compare);
swap_if(first[0u], first[12u], compare);
swap_if(first[2u], first[14u], compare);
swap_if(first[4u], first[16u], compare);
swap_if(first[6u], first[18u], compare);
swap_if(first[8u], first[20u], compare);
swap_if(first[10u], first[22u], compare);
swap_if(first[2u], first[12u], compare);
swap_if(first[10u], first[20u], compare);
swap_if(first[4u], first[12u], compare);
swap_if(first[6u], first[14u], compare);
swap_if(first[8u], first[16u], compare);
swap_if(first[10u], first[18u], compare);
swap_if(first[8u], first[12u], compare);
swap_if(first[10u], first[14u], compare);
swap_if(first[10u], first[12u], compare);
swap_if(first[1u], first[13u], compare);
swap_if(first[3u], first[15u], compare);
swap_if(first[5u], first[17u], compare);
swap_if(first[7u], first[19u], compare);
swap_if(first[9u], first[21u], compare);
swap_if(first[3u], first[13u], compare);
swap_if(first[11u], first[21u], compare);
swap_if(first[5u], first[13u], compare);
swap_if(first[7u], first[15u], compare);
swap_if(first[9u], first[17u], compare);
swap_if(first[11u], first[19u], compare);
swap_if(first[9u], first[13u], compare);
swap_if(first[11u], first[15u], compare);
swap_if(first[11u], first[13u], compare);
swap_if(first[11u], first[12u], compare);

可能有更聪明的方法来生成中值查找网络，但我认为尚未对这个主题进行过广泛的研究。因此，这可能是您目前可以使用的最佳方法。结果并不好，但它仍然使用 104 个比较交换单元而不是 118 个。

Answer 2

总体思路

查看MSVC2013中std::nth_element的源代码，似乎N 的情况是通过插入排序解决的。这意味着 STL 实现者意识到，尽管对于该大小有更好的渐近性，但进行随机分区会更慢。

提高性能的方法之一是优化排序算法。 @Morwenn's answer 展示了如何使用排序网络对 23 个元素进行排序，这是已知的对小型恒定大小数组进行排序的最快方法之一。 我将研究另一种方法，即在没有排序算法的情况下计算中位数。事实上，我根本不会置换输入数组。

既然我们谈论的是小数组，我们需要以最简单的方式实现一些O(N^2)算法。理想情况下，它应该根本没有分支，或者只有可预测的分支。此外，算法的简单结构可以让我们对其进行矢量化，进一步提高其性能。

算法

我决定采用计数方法，使用here 来加速小型线性搜索。首先，假设所有元素都是不同的。选择数组的任何元素：元素的数量小于它在排序数组中定义的位置。我们可以遍历所有元素，并为每个元素计算小于它的元素数。如果排序后的索引有期望的值，我们可以停止算法。

不幸的是，在一般情况下可能存在相等的元素。我们必须使我们的算法明显更慢和更复杂来处理它们。我们可以计算它的可能排序索引的间隔，而不是计算元素的唯一排序索引。对于任何元素，计算小于它的元素数（L）和等于它的元素数（E）就足够了，然后排序索引适合范围[左，左+右）。如果此区间包含所需的排序索引（即 N/2），那么我们可以停止算法并返回考虑的元素。

for (size_t i = 0; i < n; i++) {
    auto x = arr[i];
    //count number of "less" and "equal" elements
    int cntLess = 0, cntEq = 0;
    for (size_t j = 0; j < n; j++) {
        cntLess += arr[j] < x;
        cntEq += arr[j] == x;
    }
    //fast range checking from here: https://stackoverflow.com/a/17095534/556899
    if ((unsigned int)(idx - cntLess) < cntEq)
        return x;
}

矢量化

构造的算法只有一个分支，这是相当可预测的：它在所有情况下都会失败，除了我们停止算法的唯一情况。该算法很容易使用每个 SSE 寄存器的 8 个元素进行矢量化。由于我们必须在最后一个元素之后访问一些元素，因此我假设输入数组填充了 max=2^15-1 值，最多 24 或 32 个元素。

第一种方法是通过j 向量化内循环。在这种情况下，内循环将只执行 3 次，但在完成后必须执行两次 8 宽的归约。他们比内循环本身吃更多的时间。结果，这样的向量化不是很有效。

第二种方法是通过i 向量化外循环。在这种情况下，我们一次处理 8 个元素 x = arr[i]。对于每个包，我们将其与内部循环中的每个元素arr[j] 进行比较。在内循环之后，我们对整个 8 个元素包执行矢量化范围检查。如果其中任何一个成功，我们使用简单的标量代码确定确切的数字（无论如何它消耗的时间很少）。

__m128i idxV = _mm_set1_epi16(idx);
for (size_t i = 0; i < n; i += 8) {
    //load pack of 8 elements
    auto xx = _mm_loadu_si128((__m128i*)&arr[i]);
    //count number of less/equal elements for each element in the pack
    __m128i cntLess = _mm_setzero_si128();
    __m128i cntEq = _mm_setzero_si128();
    for (size_t j = 0; j < n; j++) {
        __m128i vAll = _mm_set1_epi16(arr[j]);
        cntLess = _mm_sub_epi16(cntLess, _mm_cmplt_epi16(vAll, xx));
        cntEq = _mm_sub_epi16(cntEq, _mm_cmpeq_epi16(vAll, xx));
    }
    //perform range check for 8 elements at once
    __m128i mask = _mm_andnot_si128(_mm_cmplt_epi16(idxV, cntLess), _mm_cmplt_epi16(idxV, _mm_add_epi16(cntLess, cntEq)));
    if (int bm = _mm_movemask_epi8(mask)) {
        //range check succeeds for one of the elements, find and return it 
        for (int t = 0; t < 8; t++)
            if (bm & (1 << (2*t)))
                return arr[i + t];
    }
}

在这里，我们在最里面的循环中看到 _mm_set1_epi16 内在函数。 GCC 似乎有一些性能问题。无论如何，每次最内层迭代都在消耗时间，如果我们在最内层循环中一次处理 8 个元素也可以减少时间。在这种情况下，我们可以执行 1 个矢量化加载和 14 个解包指令来获得 8 个元素的vAll。此外，我们必须为循环体中的八个元素编写比较和计数代码，因此它也可以作为 8x 展开。生成的代码是最快的，可以在下面找到它的链接。

比较

我在 Ivy Bridge 3.4 Ghz 处理器上对各种解决方案进行了基准测试。您可以在下面看到 2^23 ~= 8M 次调用的总计算时间（以秒为单位）（第一个数字）。第二个数字是结果的校验和。

MSVC 2013 x64 (/O2) 上的结果：

memcpy only: 0.020
std::nth_element: 2.110 (1186136064)
network sort: 0.630 (1186136064)              //solution by @Morwenn (I had to change swap_if)
trivial count: 2.266 (1186136064)             //scalar algorithm (presented above)
vectorized count: 0.692 (1186136064)          //vectorization by j
vectorized count (T): 0.602 (1186136064)      //vectorization by i (presented above)
vectorized count (both): 0.450 (1186136064)   //vectorization by i and j

MinGW GCC 4.8.3 x64 (-O3 -msse4) 上的结果：

memcpy only: 0.016
std::nth_element: 1.981 (1095237632)
network sort: 0.531 (1095237632)              //original swap_if used
trivial count: 1.482 (1095237632)
vectorized count: 0.655 (1095237632)
vectorized count (T): 2.668 (1095237632)      //GCC generates some crap
vectorized count (both): 0.374 (1095237632)

如您所见，针对 23 个 16 位元素提出的矢量化算法比基于排序的方法要快一点（顺便说一句，在较旧的 CPU 上，我只看到 5% 的时间差异）。 如果你能保证所有元素都是不同的，你就可以简化算法，让它变得更快。

所有算法的完整代码可在here获取，包括所有测试代码。

Answer 3

我发现这个问题很有趣，所以我尝试了所有我能想到的算法。
结果如下：

testing 100000 repetitions
variant 0, no-op (for overhead measure)
5 ms
variant 1, vector + nth_element
205 ms
variant 2, multiset + advance
745 ms
variant 2b, set (not fully conformant)
787 ms
variant 3, list + lower_bound
589 ms
variant 3b, list + block-allocator
269 ms
variant 4, avl-tree + insert_sorted
645 ms
variant 4b, avl-tree + prune
682 ms
variant 5, histogram
1429 ms

我想我们可以得出结论，你已经在使用最快的了 算法。 男孩是我错了。但是，如果你能接受一个大概的答案， 可能有更快的方法，例如median of medians。
如有兴趣，来源here。