3

So I know how to do a fast O(log2(N)) sliding max or sliding min algorithm.

Brookes: "Algorithms for Max and Min Filters with Improved Worst-Case Performance" IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—II: ANALOG AND DIGITAL SIGNAL PROCESSING, VOL. 47, NO. 9, SEPTEMBER 2000

But I just cannot figure out how to do a sliding median without first sorting the buffer (or at least sorting it half way) and then picking the value halfway into the buffer (or midway between the bottom and the top of the sorted buffer).

Given a finite-length buffer (length N) that is chronologically sorted (like a FIFO) and another with the same N values sorted by value, then when the buffer window "slides" by one and the oldest value falls offa the edge and a new one is input, to insert that new value into the sorted buffer, that is an O(N) operation (probabilistically).

What's a faster way to do this?

6
  • i've stumbled upon "quickselect" just now but haven't yet grok'd it. if this is the current prevailing wisdom, i would appreciate if someone confirms that.
    – robert bristow-johnson
    Commented Apr 28, 2015 at 0:28
  • and, i don't think quickselect takes any advantage of the fact that this is a sliding median, so i already have the sorted buffer from the previous time, one new value is input, one old value is kicked out. we can insert the new one and remove the old one from the sorted buffer with an O(N) operation. quickselect appears to be O(N log2(N)) which is slower than O(N).
    – robert bristow-johnson
    Commented Apr 28, 2015 at 0:34
  • 2
    QuickSelect seems to do what it says on the tin. How many values are in this buffer? Because if it's not an enormous number of values, O(n) will work anyway. It's hard to imagine any algorithm that calculates median that doesn't touch every element.
    – Robert Harvey
    Commented Apr 28, 2015 at 1:36
  • @RobertHarvey I think you hit the nail on the head with that statement. Median says something about all values, but how can you find it without knowing all values? You could sample a percentage and get an approximate median, but that's still technically O(n) time. I think the algorithm doesn't take sorting into consideration.
    – Neil
    Commented Apr 28, 2015 at 7:30
  • 1
    @RobertHarvey, N can be as large as 2K or maybe 8K. this is single dimensional, so i will have to grok the 2D median algs that the answer below refers to. someone on comp.dsp referred me to Huang et. al.: A Fast Two-Dimensional Median Filtering Algorithm and another paper i haven't looked at yet. i'm starting to grok the possibility of an approximate median based on some running histogram.
    – robert bristow-johnson
    Commented Apr 28, 2015 at 18:41

3 Answers 3

Reset to default
2

See https://stackoverflow.com/q/5527437/10396 for an implementation of a rolling median that uses a min-max heap to process each new sample in O(lgN). The heap keeps the data loosely sorted into two groups, one bigger than the median, one smaller. For each new sample, it swaps the the oldest item in the heap with the newest one. Rebalancing the heap takes up to 2 O(lgN) steps: a sift-up to make sure the new item is closer to the median than any parents, and maybe a sift-down to push it down the other side.

5
  • thanks AShelly, but i wonder what the definition (and implementation) of maxSortDown() and minSortDown() are. and maxSortUp() and minSortUp(). what are they and how much do they cost?
    – robert bristow-johnson
    Commented Apr 28, 2015 at 20:02
  • okay, i found the definitions here. they still look like O(N), not O(log2(N)). can you show me how these component tools are less costly than O(N)?
    – robert bristow-johnson
    Commented Apr 28, 2015 at 20:11
  • Every new sample is inserted into a heap, which has lgN insert and delete, because every time you go down a level, you are eliminating half the candidates from consideration.
    – AShelly
    Commented Apr 28, 2015 at 20:22
  • hay @AShelly, thank you for your code. i was wondering if it would be more efficient to maintain two buffers for the two sets of samples above the median (on which we apply a sliding max) and below the median (one which we apply the sliding min algorithm). both halves are just as long as the delay and samples in the upper half that are less than the median are marked as +3.4e38 and samples in the lower half are marked as -3.4e38 so the don't get in the min or max search. then there needs to be a mechanism to transfer samples from one half to the other when needed.
    – robert bristow-johnson
    Commented Dec 19, 2016 at 1:47
  • similar to your binary search, there is this code that i wrote that does sliding min and sliding max very efficiently. there isn't the need for indirection, or that Mediator struct where each sample index in the queue is saved. when a new sample goes into the same half that the old sample that drops offa the edge, then there is no need to transfer a sample from one buffer to the other. it's when the two samples go into/come outa the different halves that a transfer must be made.
    – robert bristow-johnson
    Commented Dec 19, 2016 at 1:53
2

Lets assume N=100. You have the raw values in chronological order raw[1],...,raw[100] and you have the ordered list of the same values as ordered[1],...,ordered[100]. What is the median? The median will be (ordered[50]+ordered[51])/2.

Lets now move the window one position so raw data will be raw[2],...,raw[101]. How do you generate an updated list for ordered[1],...,ordered[100]? Simple. Do a binary search for raw[1] value in the ordered array (ordered[1],...,ordered[100]) and remove it from the list. I am assuming a binary search tree or skip list implementation so you don't need to move all values higher than raw[1] one position down. Now do another binary search for raw[101] value in the ordered array and add it to the list. I am again assuming a binary search tree or skip list so you don't need to move all values higher than raw[101] one position up.

How you calculate the median now? The same way as before, it will be (ordered[50]+ordered[51])/2. But now ordered[50] and ordered[51] may be different values now.

So the sliding median is O(Log(N)) because that is what a binary search takes which is what you are looking for. Yes, you need to sort the first N values, but once you do that, you just keep the numbers ordered as you slide thought the chronological values (signal). In the long run, the time you took to order those first values will just be diluted.

7
  • OP is looking for a way to do this "without first sorting the buffer ... and then picking the value halfway into the buffer". So he's probably already aware of this method.
    – Ixrec
    Commented Apr 28, 2015 at 6:36
  • yeah, Mandrill, the way i know how to do it now has two binary searches, one to find the value of raw[1] and kick it out, which creates a hole. the other binary search is for the location to insert the value of raw[101]. both cost O(log2(N)). then between the locations of the insert and the hole in the sorted array, i must bump all values over by 1 place. that is a probabilistic number of places anywhere from 0 to 99. i think that cost is O(N).
    – robert bristow-johnson
    Commented Apr 28, 2015 at 18:50
  • @robertbristow-johnson Note that this answer included the sentence "I am assuming a binary search tree or skip list implementation so you don't need to move all values higher than raw[1] one position down." You don't have the array O(N) performance hit if you don't use an array.
    – Ordous
    Commented Apr 28, 2015 at 19:00
  • okay, @Ordous, being an EE and not CS, i am not aware of "skip list implementation". i can tell you about all sorts of DSP algs, but these sorting algs are not in my skill set. i dunno what a "skip list" is. (care to tell me about it or to refer to a good primer or tutorial?)
    – robert bristow-johnson
    Commented Apr 28, 2015 at 19:04
  • 3
    @robertbristow-johnson BST, Skip List. Both are quite common data-structures for sorted data that have O(Log(N)) lookup, insertion and deletion times.
    – Ordous
    Commented Apr 28, 2015 at 19:07
2

For certain types of data, you can achieve constant time median filtering (Perreault et al, 2007).

That paper describes 2D median filtering on images, assuming the pixels are 8-bit integers.

Note that "constant time" refers to constant time in the size of the window; it is not constant time in the size of data, or in the precision (bits) of data. As explained below, using this algorithm with high-precision data will dramatically increase the memory usage by the algorithm, due to the histogram.

Firstly it is necessary to understand multi-level histograms.

  • When the alphabet set (the set of allowed pixel values) has 256 symbols, the histogram has 256 bins.
  • A multi-level histogram for 256 bins will have two levels. The first level will have 16 bins, and the second level will have the full 256 bins.
  • Each of the 16 first-level bin corresponds to some 16 consecutive bins on the second level, in the most straightforward way.
  • Each increment (add-sample) operation will increment the bin on the second level, as well as one of the 16 bins on the first level.
  • Likewise, decrement (remove-sample) operation will decrement a bin on both the first and second levels.
  • When one needs to search for the median, one will first search among the first-level bins, as there are fewer of them. Once the correct first-level bin is found, one will then search among the 16 second-level bins, where it is guaranteed that one of those will contain the median.
  • From this, you can see why the algorithm time complexity is proportional to the number of bits in the data precision.
  • Then, for each dimension of data, you will maintain the histogram for a sliding window, each time adding one sample and removing one sample.
  • If the data has multiple dimensions, then an array of histograms (as many as the width of the input in the first dimension) will need to be maintained.

Obviously this quickly grows out of hand for higher dimensions, but for some combinations of input sizes and window sizes, this scheme turns out to have better performance than previously known approaches.

Whether this algorithm can be used in your application will depend on the type of data you need to sort.

For example, if your data contains 32-bit integers, and all 2^32 different values are equally likely to appear in the data, then you will need a histogram having 2^32 bins, which is a somewhat insane requirement, but might still be doable. If your data contains IEEE double-precision floating point values, and all 1023 * 2^52 values are likely to occur, then the histogram will apparently not fit in any kind of computers currently available.

You can reduce the number of histogram bins needed by lowering the resolution (precision) of your data.

3
  • I can imagine using a Bloom filter or a HyperLogLog counter to approximate the 2^32 histogram, and maybe keeping more levels of the histogram.
    – 9000
    Commented Apr 28, 2015 at 20:19
  • @9000: I'm not exactly sure how Bloom filter or HyperLogLog can be used for the purpose of calculating the median. Your comment on increasing the levels of multi-level histogram is correct; I will edit my answer.
    – rwong
    Commented Apr 28, 2015 at 20:45
  • Well, yes, on the second thought, neither of them are applicable here. Sorry.
    – 9000
    Commented Apr 28, 2015 at 20:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.