# Nontrivial algorithm for computing a sliding window median

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.)

-
I don't think you'll be able to improve $nlogk$. The reason is, if you look at a window $x_t,...,x_{t+k-1}$, you can never rule out any of the numbers $x_{t+\frac{k}{2}},...,x_{t+k-1}$ from being medians of future window. This means that in any time you have to keep at least $\frac{k}{2}$ integers in a data structure, and it doesn't seem to update in less than log time. – R B Mar 24 '14 at 17:55
Your trivial algorithm to me seems to be $O((n-k)\cdot k \cdot \log k)$ not $O(n \log k)$, am I misunderstood something? And I think because of this you have problem with big $k$, otherwise logarithmic factor is nothing in practical applications, also there is no big hidden constant in this algorithm. – Saeed Mar 24 '14 at 22:53
@Saeed: In the trivial algorithm, you process elements one by one; in step $i$ you add $x_i$ to the search tree and (if $i > k$) you also remove $x_{i-k}$ from the search tree. This is $n$ steps, each of which takes $O(\log k)$ time. – Jukka Suomela Mar 24 '14 at 22:58
So you mean you have a balanced search tree not casual search tree? – Saeed Mar 24 '14 at 23:02
@Saeed: Please note that in my benchmarks I did not even try to find medians. I just did $n$ insertions and $n$ deletions in a search tree of size $k$, and these operations are guaranteed to take $O(\log k)$ time. You just need to accept that search tree operations are very slow in practice, in comparison with sorting. You will see this easily if you try to write a sorting algorithm that works by adding elements to a balanced search tree — it certainly works in $O(n \log n)$ time, but it will be ridiculously slow in practice, and also waste a lot of memory. – Jukka Suomela Mar 30 '14 at 10:48

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.

-
This answer really makes me wonder if the converse holds as well: given an efficient sorting algorithm, do we get an efficient running median algorithm? (For example, does on efficient integer sorting algorithm imply an efficient running median algorithm for integers? Or does an IO-efficient sorting algorithm provide an IO-efficient running median algorithm?) – Jukka Suomela Mar 27 '14 at 14:59
Once again, many thanks for your answer, it really set me on the right track and gave inspiration for the sorting-based median filter algorithm! In the end I was able to find a paper from 1991 that presented basically the same argument as what you give here, and Pat Morin gave a pointer to another relevant paper from 2005; see refs. [6] and [9] here. – Jukka Suomela Jun 9 '14 at 8:42

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:

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.

-

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.

-
Oops, this doesn't work as written, since if you don't delete elements the counts won't reflect the new window. I'm not sure if it can be fixed, but I will leave the answer in case there is a way. – Geoffrey Irving Mar 25 '14 at 6:15
So I think this algorithm may in fact take $O(n \log m)$ if you delete nodes from the Fibonacci heaps, since the Fibonacci heap depth increases only when delete-min is called. Does anyone know nice bounds on Fibonacci heap complexity taking the number of delete-min calls into account? – Geoffrey Irving Mar 25 '14 at 6:58
side note: Question is not clear, underling data structure is not defined, we just know something very vague. how do you want improve something that you don't know what it is? how do you want compare your approach? – Saeed Mar 29 '14 at 0:18
I apologize for the incomplete work. I've asked the concrete question needed to fix this answer here: cstheory.stackexchange.com/questions/21778/…. If you think it's appropriate I can remove this answer until the secondary question is resolved. – Geoffrey Irving Mar 29 '14 at 1:14