http://dpaste.dzfl.pl/05a82699acc8
So over the past few days I've been in the zone working on a smooth
implementation of "Median of Medians"
(https://en.wikipedia.org/wiki/Median_of_medians). Its performance is
much better compared to the straightforward implementation. However, in
practice it's very hard to get it to beat simple heuristics (such as
median of five, random pivot etc). These heuristics run `n` comparisons
for an input of length `n`, whereas the MoM needs over twice as many
comparisons. Also it does a bunch of swapping around the data. Overall I
got it within 2.5x of the heuristic-based topN for most data sizes up to
tens of millions.
While thinking of MoM and the core reasons of its being slow (adds nice
structure to its input and then "forgets" most of it when recursing), I
stumbled upon a different algorithm. It's much simpler, also
deterministic and faster than MoM for many (most?) inputs. But it's not
guaranteed to be linear. After having pounded at this for many hours, it
is clear that I am in need of some serious due destruction. I call it a
"quick median of medians" or in short "quick mom".
Consider the algorithm defined as follows over a range `r` of length
`n`. It returns an index to an element likely to be in the second
tertile of r. (A tertile is a third of the range. So the expectation is
that the returned index x is such that e<=r[x] for at least one third of
the range, and r[x]<=e again for at least one third of the range.)
0. If n<=3, compute median by rote and return its index.
1. Divide r in three equal adjacent subranges r0 = r[0 .. $/3], r1 =
r[$/3 .. $*2/3], r2 = r[$*2/3 .. $].
2. Recurse to get the median of these three, call them m0, m1, m2.
3. Return the median of m0, m1, m2.
Note that no data has been written yet; we only have an estimate of a
pivot. On random inputs it's expected that the median thus gotten is
greater than 1/6 + 1/6 = 1/3 elements, and less than the same fraction.
The algorithm completes in linear time.
However, the pivot obtained is just an approximation. In the worst case
it's possible that e.g. all recursions return the leftmost of the
allowed index, so the bound deteriorates by 1/3 for each recursion
depth, or generally to (1/3)^^log3(n).
Not good! That's where the quick mom pays attention. After it partitions
the data using the pivot obtained by the quick mom method, the
partitioning stage checks whether the resulting pivot falls within
bounds. If not, it runs a precise selection method (such as proper
median of medians - "thorough mom") to bring the pivot where it needs to
be. The data patterns that cause the quick mom to fail systematically
are rather odd but worst case is what it is.
Overall the partitioning either succeeds with the quick mom pivot in
linear time, or does one more linear pass if "unlucky". So overall
partition is linear. (Micro-optimization: not all data needs to be
repartitioned, only that outside the not-so-good pivot.) Unlike the
quick mom, the partition guarantees a pivot in the mid tertile.
After partitioning the classic quickselect algorithm may be implemented.
There's one more _really_ juicy detail. In fact, the quick mom may
finish without looking at all data. Look at the implementation - there
are two overloads of quickMom. The second gets the bounds and works as
follows: if you've already computed two elements of a median, the third
may be chosen only if in between the two. Otherwise, you only care
whether it's larger or smaller than the other two. This allows the
algorithm to finish certain recursion branches early.
Destroy!
Andrei
