Benchmark a sorting method for fixed-size arrays

Question

During the past few days, I have written a small C++14 library which mainly provides the cppsort::sort function. It is designed to efficiently sort small arrays when their time is known at compile time, and to fall back to the standard library std::sort otherwise.

While I am rather confident that the new sorting method achieves its goal, I still wrote a benchmarking program to compare cppsort::sort and std::sort. Since I am not used to writing benchmarks, I would like this benchmarking method reviewed in order to know whether there are bias which might affect the results. Any kind of review is welcome :)

The program generates shuffled arrays for every small array size (\$0\$ to \$64\$ in our case) and sorts them one million times per algorithm and returns series of integers to the standard output which correspond to the execution times of the algorithms in milliseconds (one line corresponds to the time of one algorithm, where the \$n\$th integer in the line corresponds to the time it took to sort one million of shuffled arrays of size \$n\$).

#include <algorithm>
#include <array>
#include <chrono>
#include <ctime>
#include <iostream>
#include <iterator>
#include <numeric>
#include <random>
#include <utility>
#include <cpp-sort/sort.h>

template<
    typename RandomAccessIterable,
    typename Compare = std::less<>
>
auto std_sort(RandomAccessIterable& iterable, Compare&& compare={})
    -> void
{
    std::sort(std::begin(iterable), std::end(iterable), std::forward<Compare>(compare));
}

template<
    typename T,
    std::size_t N,
    typename SortFunction1,
    typename SortFunction2
>
auto time_compare(SortFunction1 sort1, SortFunction2 sort2, std::size_t times)
    -> std::array<std::chrono::milliseconds, 2u>
{
    // Random numbers generator
    thread_local std::mt19937_64 engine(std::time(nullptr));

    // Generate shuffled array, the same for both algorithms
    std::array<T, N> array;
    std::iota(std::begin(array), std::end(array), 0);
    std::shuffle(std::begin(array), std::end(array), engine);

    // Time first algorithm
    auto start = std::chrono::high_resolution_clock::now();
    for (std::size_t i = 0 ; i < times ; ++i)
    {
        auto unsorted = array;
        sort1(unsorted, std::less<>{});
    }
    auto end = std::chrono::high_resolution_clock::now();
    auto duration1 = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);

    // Time second algorithm
    start = std::chrono::high_resolution_clock::now();
    for (std::size_t i = 0 ; i < times ; ++i)
    {
        auto unsorted = array;
        sort2(unsorted, std::less<>{});
    }
    end = std::chrono::high_resolution_clock::now();
    auto duration2 = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);

    return { duration1, duration2 };
}


template<typename T, std::size_t... Ind>
auto time_them(std::size_t size, std::index_sequence<Ind...>)
    -> std::array<
            std::array<
                std::chrono::milliseconds,
                sizeof...(Ind)
            >,
            2u
        >
{
    // Benchmark the sorts std::sort
    std::array<std::array<std::chrono::milliseconds, 2u>, sizeof...(Ind)> results = {
        time_compare<T, Ind>(
            &cppsort::sort<T, Ind, std::less<>>,
            &std_sort<std::array<T, Ind>, std::less<>>,
            size
        )...
    };

    // Results for cppsort::sort
    std::array<std::chrono::milliseconds, sizeof...(Ind)> first = {
        std::get<Ind>(results)[0u]...
    };

    // Results for std::sort
    std::array<std::chrono::milliseconds, sizeof...(Ind)> second = {
        std::get<Ind>(results)[1u]...
    };

    return { first, second };
}

int main()
{
    using indices = std::make_index_sequence<64u>;
    auto sorts_times = time_them<int>(1000000u, indices{});

    for (auto&& sort_times: sorts_times)
    {
        for (auto&& time: sort_times)
        {
            std::cout << time.count() << ' ';
        }
        std::cout << '\n';
    }
}

To use it, I call benchmark.exe > results.txt then forward the results to a Python script which displays them on a graph with python plot.py results.txt. Here is the said Python script:

import sys
import matplotlib.pyplot as plt


def fetch_results(fresults):
    results = fresults.readline().split(' ')
    results.pop()
    return [int(elem) for elem in results]


if __name__ == '__main__':
    # Results of timing functions
    cppsort = []
    stdsort = []

    # Fetch the results
    with open(sys.argv[1]) as f:
        cppsort = fetch_results(f)
        stdsort = fetch_results(f)

    # Plot the results
    xaxis = list(range(len(cppsort)))
    line_cpp, = plt.plot(xaxis, cppsort)
    line_std, = plt.plot(xaxis, stdsort)
    plt.legend([line_cpp, line_std], ['cppsort::sort', 'std::sort'])
    plt.xlabel('Number of elements to sort')
    plt.ylabel('Execution time (ms)')
    plt.show()

It generates the following kind of graph:

As we can see, it does it job to "prove" that cppsort::sort is more performant than std::sort for values of \$n\$ smaller than \$32\$ and it also makes it pretty obvious that the algorithm falls back to std::sort for bigger values. That said I would like to know how I could stil improve the benchmarking method (is it possible to make the curve smoother?), and whether there is any bias in it.

---

Note: how the cppsort::sort algorithm is implemented is off-topic for this question; you can find the explanation in the linked library's README. That said, you will at least need its interface to test the benchmarking method. For the sake of simplicity, I will provide a function with the exact same interface, but it will only forward the job to std::sort instead. Trust me, you don't want 10k additional lines of code in the question.

namespace cppsort
{
    template<
        typename RandomAccessIterable,
        typename Compare = std::less<>
    >
    auto sort(RandomAccessIterable& iterable, Compare&& compare={})
        -> void
    {
        std::sort(
            std::begin(iterable),
            std::end(iterable),
            std::forward<Compare>(compare)
        );
    }

    template<
        typename T,
        std::size_t N,
        typename Compare = std::less<>
    >
    auto sort(std::array<T, N>& array, Compare&& compare={})
        -> void
    {
        std::sort(
            std::begin(array),
            std::end(array),
            std::forward<Compare>(compare)
        );
    }
}

Again, for the sake of simplicity, I didn't provide the cppsort::sort overload for fixed-size C arrays since it's not tested by the benchmarking method.

Answer 1

I don't have time to do a full review but I will note on one important thing:

// Time first algorithm
auto start = std::chrono::high_resolution_clock::now();
for (std::size_t i = 0 ; i < times ; ++i)
{
    auto unsorted = array;
    sort1(unsorted, std::less<>{});
}
auto end = std::chrono::high_resolution_clock::now();
auto duration1 = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);

Here, as unsorted is never used, the compiler is allowed under the "as-if" rule to completely omit generating any code for the for loop. And a good optimizing compiler will. This is indicated by the 0 time to sort for small values in your results. Probably the compiler can do this for your code but not for std::sort.

What's more the compiler can deduce that the entire loop body is loop invariant and just calculate the result once.

For your test to be accurate you need to:

1. Generate n shuffled arrays in a vector (std::vector<std::vector<T>>).
2. Start timing
3. For each of the n shuffled arrays, sort it in place.
4. Stop timing
5. Write all the n arrays somewhere (cout or file) to prevent removal by as-if.

Answer 2

• Your code never even looks at the sorted array. An all to smart compiler could figure out that the whole sorting is useless and optimize it out. You could avoid this, for example, by assigning a randomly picked element of the sorted array to a volatile variable. Don't make the whole array volatile, however, or you'll lose a lot of useful optimizations.

• You are sorting the same permutation over and over again. You'd get more representative results if you'd use different permutations. This is important especially since the run-time of some sorting algorithms depends on the inputs. Again, a very smart compiler might also figure out that sorting the same data again is pointless. Since the expected run-time of sorting such small arrays is tiny, it might seem reasonable to time a larger number of runs to become less dependent on the clock accuracy. You could, however, sort the same array again using a different comparator. For example, you could XOR a random value to the values before comparing them. Pick a different value each time you sort.

• The second algorithm might have an unfair advantage since the data will already be cached. If you are worried about this, randomize the order in which the algorithms run.

• It is implementation defined whether std::chrono::high_resolution_clock is a real-time clock or steady. If it is not steady and you have running an NTP daemon in the background, your timing results will be worthless. You could use std::chrono::steady_clock instead or statically select the best available clock at compile-time. In your case, it might also be worth consideration to use std::clock instead, if you are not interested in system effects.

• Instead of only plotting the average values, compute mean and standard deviation and plot them as data points with error bars. Then overlay a weighted linear regression of α + β n log(n) (or whatever curve you deem appropriate). This will give you the desired “smoother” and more honest curve.