Nontrivial algorithm for computing a sliding window median

Question

I need to calculate the running median:

• Input: $n$, $k$, vector $(x_1, x_2, \dotsc, x_n)$.

• Output: vector $(y_1, y_2, \dotsc, y_{n-k+1})$, where $y_i$ is the median of $(x_i, x_{i+1}, \dotsc, x_{i+k-1})$.

(No cheating with approximations; I would like to have exact solutions. Elements $x_i$ are large integers.)

There is a trivial algorithm that maintains a search tree of size $k$; the total running time is $O(n \log k)$. (Here a "search tree" refers to some efficient data structure that supports insertions, deletions, and median queries in logarithmic time.)

However, this seems a bit stupid to me. We will effectively learn all order statistics within all windows of size $k$, not just the medians. Moreover, this is not too attractive in practice, especially if $k$ is large (large search trees tend to be slow, overhead in memory consumption is non-trivial, cache-efficiency is often poor, etc.).

Can we do anything substantially better?

Are there any lower bounds (e.g., is the trivial algorithm asymptotically optimal for the comparison model)?

---

Edit: David Eppstein gave a nice lower bound for the comparison model! I wonder if it is nevertheless possible to do something slightly more clever than the trivial algorithm?

For example, could we do something along these lines: divide the input vector to parts of size $k$; sort each part (keeping track of the original positions of each element); and then use the piecewise sorted vector to find the running medians efficiently without any auxiliary data structures? Of course this would still be $O(n \log k)$, but in practice sorting arrays tends to be much faster than maintaining search trees.

---

Edit 2: Saeed wanted to see some reasons why I think sorting is faster than search tree operations. Here are very quick benchmarks, for $k = 10^7$, $n = 10^8$:

• ≈ 8s: sorting $n/k$ vectors with $k$ elements each
• ≈ 10s: sorting a vector with $n$ elements
• ≈ 80s: $n$ insertions & deletions in a hash table of size $k$
• ≈ 390s: $n$ insertions & deletions in a balanced search tree of size $k$

The hash table is there just for comparison; it is of no direct use in this application.

In summary, we have almost a factor 50 difference in the performance of sorting vs. balanced search tree operations. And things get much worse if we increase $k$.

(Technical details: Data = random 32-bit integers. Computer = a typical modern laptop. The test code was written in C++, using the standard library routines (std::sort) and data structures (std::multiset, std::unsorted_multiset). I used two different C++ compilers (GCC and Clang), and two different implementations of the standard library (libstdc++ and libc++). Traditionally, std::multiset has been implemented as a highly optimised red-black tree.)

Answer 1

Here's a lower bound from sorting. Given an input set $S$ of length $n$ to be sorted, create an input to your running median problem consisting of $n-1$ copies of a number smaller than the minimum of $S$, then $S$ itself, then $n-1$ copies of a number larger than the maximum of $S$, and set $k=2n-1$. The running medians of this input are the same as the sorted order of $S$.

So in a comparison model of computation, $\Omega(n\log n)$ time is required. Possibly if your inputs are integers and you use integer sorting algorithms you can do better.

Answer 2

Edit: This algorithm is now presented here: http://arxiv.org/abs/1406.1717

---

Yes, to solve this problem it is sufficient to perform the following operations:

• Sort $n/k$ vectors, each with $k$ elements.
• Do linear-time post-processing.

Very roughly, the idea is this:

• Consider two adjacent blocks of input, $a$ and $b$, both with $k$ elements; let the elements be $a_1, a_2, ..., a_k$ and $b_1, b_2, ..., b_k$ in the order of appearance in input vector $x$.
• Sort these blocks and learn the rank of each element within the block.
• Augment the vectors $a$ and $b$ with predecessor/successor pointers so that by following the pointer chains we can traverse the elements in an increasing order. This way we have constructed doubly linked lists $a'$ and $b'$.
• One by one, delete all elements from the linked list $b'$, in the reverse order of appearance $b_k, b_{k-1}, ..., b_1$. Whenever we delete an element, remember what was its successor & predecessor at the time of deletion.
• Now maintain "median pointers" $p$ and $q$ that point to lists $a'$ and $b'$, respectively. Initialise $p$ to the midpoint of $a'$, and initialise $q$ to the tail of the empty list $b'$.

• For each $i$:

  • Delete $a_i$ from list $a'$ (this is $O(1)$ time, just delete from the linked list). Compare $a_i$ with the element pointed by $p$ to see if we deleted before or after $p$.
  • Put $b_i$ back to list $b'$ in its original position (this is $O(1)$ time, we memorised the predecessor and successor of $b_i$). Compare $b_i$ with the element pointed by $q$ to see if we added the element before or after $q$.
  • Update pointers $p$ and $q$ so that the median of the joined list $a' \cup b'$ is either at $p$ or at $q$. (This is $O(1)$ time, just follow the linked lists one or two steps to fix everything. We will keep track of how many items are before/after $p$ and $q$ in each list, and we will maintain the invariant that both $p$ and $q$ point to elements that are as close to the median as possible.)

The linked lists are just $k$-element arrays of indexes, so they are lightweight (except that the locality of memory access is poor).

---

Here is a sample implementation and benchmarks:

• https://github.com/suomela/median-filter

Here is a plot of running times (for $n \approx 2\cdot 10^6$):

• Blue = sorting + post-processing, $O(n \log k)$.
• Green = maintain two heaps, $O(n \log k)$, implementation from https://github.com/craffel/median-filter
• Red = maintain two search trees, $O(n \log k)$.
• Black = maintain a sorted vector, $O(n k)$.
• X axis = window size ($\approx k/2$).
• Y axis = running time in seconds.
• Data = 32-bit integers and random 64-bit integers, from various distributions.

Answer 3

Given David's bound it's unlikely you can do better worst case, but there are better output sensitive algorithms. Specifically, if $m$ in the number of medians in the result, we can solve the problem in time $O(n \log m + m \log n)$.

To do this, replace the balanced binary tree with a balanced binary tree consisting of only those elements that were medians in the past, plus two Fibonacci heaps in between each pair of previous medians (one for each direction), plus counts so that we can locate which Fibonacci heap contains a particular element in the order. Don't bother ever deleting elements. When we insert a new element, we can update our data structure in $O(\log m)$ time. If the new counts indicate that the median is in one of the Fibonacci heaps, it takes an additional $O(\log n)$ to pull the new median out. This $O(\log n)$ charge occurs only once per median.

If there was a clean way to delete elements without damaging the nice Fibonacci heap complexity, we'd get down to $O(n \log m + m \log k)$, but I'm not sure if this is possible.