Download: STREAMING MAXIMUM-MINIMUM FILTER USING NO MORE THAN THREE COMPARISONS PER ELEMENT

STREAMING MAXIMUM-MINIMUM FILTER USING NO MORE THAN THREE COMPARISONS PER ELEMENT

Daniel Lemire University of Quebec at Montreal (UQAM), UER ST 100 Sherbrooke West, Montreal (Quebec), H2X 3P2 Canada email is hidden Abstract. The running maximum-minimum (-) filter computes the maxima and minima over running windows of size w. This filter has numerous applications in signal processing and time series analysis. We present an easy-to-implement online algorithm requiring no more than 3 comparisons per element, in the worst case. Comparatively, no algorithm is known to compute the running maximum (or minimum) filter in 1.5 com- parisons per element, in the worst case. Our algorithm has reduced latency and memory usage. ACM CCS Categories and Subject Descriptors: F.2.1 Numerical Algorithms and Prob- lems Key words: Design of Algorithms, Data Streams, Time Series, Latency, Monotonicity 1. Introduction The maximum and the minimum are the simplest form of order statistics. Com- puting either the global maximum or the global minimum of an array of n ele- ments requires n − 1 comparisons, or slightly less than one comparison per ele- ment. However, to compute simultaneously the maximum and the minimum, only 3⌈n/2⌉ − 2 comparisons are required in the worst case [Cormen et al. 2001], or slightly less than 1.5 comparisons per element. A related problem is the computation of the running maximum-minimum (- ) filter: given an array a1, ., an, find the maximum and the minimum over all windows of size w, that is max /mini∈[ j, j+w) ai for all j (see Fig. 1.1). The running maximum () and minimum () filters are defined similarly. The - filter problem is harder than the - problem, but a tight bound on the number of comparisons required in the worst case remains an open problem. Running maximum-minimum (-) filters are used in signal processing and pattern recognition. As an example, Keogh and Ratanamahatana [2005] use a pre- computed - filter to approximate the time warping distance between two time series. Time series applications range from music retrieval [Zhu and Shasha 2003] to network security [Sun et al. 2004]. The unidimensional - filter can be applied to images and other bidimensional data by first applying the uni- dimensional on rows and then on columns. Image processing applications include cancer diagnosis [He et al. 2005], character [Ye et al. 2001] and handwriting [Ye Received February 1, 2008. arXiv:cs/0610046v5 [cs.DS] 22 Mar 2007, 2 DANIEL LEMIRE running maximum running minimum Fig. 1.1: Example of a running - filter. et al. 2001] recognition, and boundary feature comparison [Taycher and Garakani 2004]. We define the stream latency of a filter as the maximum number of data points required after the window has passed. For example, an algorithm requiring that the whole data set be available before the running filter can be computed has a high stream latency. In effect, the stream latency is a measure of an algorithm on the batch/online scale. We quantify the speed of an algorithm by the number of comparisons between values, either a < b or b < a, where values are typically floating-point numbers. We present the first algorithm to compute the combined - filter in no more than 3 comparisons per element, in the worst case. Indeed, we are able to save some comparisons by not treating the - filter as the aggregate of the and filters: if x is strictly larger than k other numbers, then there is no need to check whether x is smaller than any of these numbers. Additionally, it is the first algorithm to require a constant number of comparisons per element without any stream latency and it uses less memory than competitive alternatives. Further, our algorithm requires no more than 2 comparisons per element when the input data is monotonic (either non-increasing or non-decreasing). We provide experimental evidence that our algorithm is competitive and can be substantially faster (by a factor of 2) when the input data is piecewise monotonic. A maybe surprising result is that our algorithm is arguably simpler to implement than the recently proposed algorithms such as Gil and Kimmel [2002] or Droogenbroeck and Buckley [2005]. Finally, we prove that at least 2 comparisons per element are required to compute the - filter when no stream latency is allowed., STREAMING MAXIMUM-MINIMUM FILTER3TI: Worst-case number of comparisons and stream latency for competitive - filter algo- rithms. Stream latency and memory usage (buffer) are given in number of elements. algorithm comparisons per ele- stream latency buffer ment (worst case) naive 2w − 2 0 O(1) van Herk [1992], Gil 6 − 8/w w 4w + O(1) and Werman [1993] Gil and Kimmel 3 + 2 log w/w w 6w + O(1) [2002] +O(1/w) New algorithm302w + O(1) 2. Related Work Pitas [1989] presented the filter algorithm requiring O(log w) com- parisons per element in the worst case and an average-case performance over inde- pendent and identically distributed (i.i.d.) noise data of slightly more than 3 com- parisons per element. Douglas [1996] presented a better alternative: the filter algorithm was shown to average 3 comparisons per element for i.i.d. input signals and Myers and Zheng [1997] presented an asynchronous implementation. More recently, van Herk [1992] and Gil and Werman [1993] presented an algo- rithm requiring 6− 8/w comparisons per element, in the worst case. The algorithm is based on the batch computation of cumulative maxima and minima over over- lapping blocks of 2w elements. For each filter ( and ), it uses a memory buffer of 2w + O(1) elements. We will refer to this algorithm as the H- G-W algorithm. Gil and Kimmel [2002] proposed an improved version (G-K) which lowered the number of comparisons per element to slightly more than 3 comparisons per element, but at the cost of some added memory us- age and implementation complexity (see Table I and Fig. 2.2 for summary). For i.i.d. noise data, Gil and Kimmel presented a variant of the algorithm requiring ≈ 2+ (2+ ln 2/2) log w/w comparisons per element (amortized), but with the same worst case complexity. Monotonic data is a worst case input for the G-K variant. Droogenbroeck and Buckley [2005] proposed a fast algorithm based on anchors. They do not improve on the number of comparisons per element. For window sizes ranging from 10 to 30 and data values ranging from 0 to 255, their implementation has a running time lower than their H-G-W implementation by as much as 30%. Their G-K implementation outperforms their H-G- W implementation by as much as 15% for window sizes larger than 15, but is outperformed similarly for smaller window sizes, and both are comparable for a window size equals to 15. The Droogenbroeck-Buckley filter pseudocode alone requires a full page compared to a few lines for H-G-W al- gorithm. Their experiments did not consider window sizes beyond w = 30 nor arbitrary floating point data values., 4 DANIEL LEMIRE van Herk 5.5 Gil-Kimmel streaming 4.5 3.5 2.5 10 20 30 40 50 60 70 80 90 100 window size Fig. 2.2: Worst-case number of comparisons per element with the H-G-W (van Herk) algorithm, the G-K algorithm, and our new streaming algorithm (less is better). 3. Lower Bounds on the Number of Comparisons Gil and Kimmel [2002] showed that the - (max /mini≤ j ai for all j) requires at least log 3 ≈ 1.58 comparisons per element, while they conjectured that at least 2 comparisons are required. We prove that their result applies directly to the - filter problem and show that 2 comparisons per element are required when no latency is allowed. T 1. In the limit where the size of the array becomes infinite, the - filter problem requires at least 2 comparisons per element when no stream latency is allowed, and log 3 comparisons per element otherwise. P. Let array values be distinct real numbers. When no stream latency is allowed, we must return the maximum and minimum of window (i − w, i] using only the data values and comparisons in [1, i]. An adversary can choose the array value ai so that ai must be compared at least twice with preceding values: it takes two comparisons with ai to determine that it is neither a maximum nor a minimum (ai ∈ (min j∈(i−w,i] a j,max j∈(i−w,i] a j)). Hence, at least 2(n − w) comparisons are required, but because 2(n −w)/n→ 2 as n→ ∞, two comparisons per element are required in the worst case. Next we assume stream latency is allowed. Browsing the array from left to right, each new data point ai for i ∈ [w, n] can be either a new maximum (ai = max j∈(i−w,i] a j), a new minimum (ai = min j∈(i−w,i] a j), or neither a new max- imum or a new minimum (ai ∈ (min j∈(i−w,i] a j,max j∈(i−w,i] a j)). For any ternary number of comparisons (worst case), STREAMING MAXIMUM-MINIMUM FILTER 5 sequence such as MAX-MAX-MIN-NOMAXMIN-MIN-MAX-., we can gener- ate a corresponding array. This means that a - filter needs to distinguish between more than 3n−w different partial orders over the values in the array a. In other words, the binary decision tree must have more than 3n−w leaves. Any binary tree having l leaves has height at least ⌈log l⌉. Hence, our binary tree must have height at least ⌈log 3n−w⌉ ≥ (n−w) log 3, proving that (1−w/n) log 3→ log 3 com- parisons per element are required when n is large. By the next proposition, we show that the general lower bound of 2 comparisons per element is tight. P 1. There exists an algorithm to compute the - filter in no more than 2 comparisons per element when the window size is 3 (w = 3), with no stream latency. P. Suppose we know the location of the maximum and minimum of the window [i − 3, i − 1]. Then we know the maximum and minimum of {ai−2, ai−1}. Hence, to compute the maximum and minimum of {ai−2, ai−1, ai}, it suffices to determine whether ai−1 > ai and whether ai−2 > ai. 4. The Novel Streaming Algorithm To compute a running - filter, it is sufficient to maintain a monotonic wedge (see Fig. 4.3). Given an array a = a1, ., an, a monotonic wedge is made of two lists U, L where U1 and L1 are the locations of global maximum and minimum, U2 and L2 are the locations of the global maximum and minimum in (U1,∞) and (L1,∞), and so on. Formally, U and L satisfy maxi>Uj1ai = aU j and min− i>L j−1 ai = aL j for j = 1, 2, .where, by convention, U0 = L0 = −∞. If all values of a are distinct, then the monotonic wedge U, L is unique. The location of the last data point n in a, is the last value stored in both U and L (see U5 and L4 in Fig. 4.3). A monotonic wedge has the property that it keeps the location of the current (global) maximum (U1) and minimum (L1) while it can be easily updated as we remove data points from the left or append them from the right: ◦ to compute a monotonic wedge of a2, a3, ., an given a monotonic wedge U, L for a1, a2, ., an, it suffices to remove (pop) U1 from U if U1 = 1 or L1 from L if L1 = 1; ◦ similarly, to compute the monotonic wedge of a1, a2, ., an, an+1, if an+1 > an, it suffices to remove the last locations stored in U until alast(U) ≥ an+1 or else, to remove the last locations stored in L until alast(L) ≤ an+1, and then to append the location n + 1 to both U and L. Fig. 4.4 provides an example of how the monotonic wedge for window [i−w, i−1] is updated into a wedge for [i − w + 1, i]. In Step A, we begin with a monotonic wedge for [i − w, i − 1]. In Step B, we add value ai to the interval. This new value is compared against the last value ai−1 and since ai > aU5 , we remove the index U5 from U. Similarly, because ai > aU4 , we also remove U4. In Step C, the index i is, 6 DANIEL LEMIRE U1 U2 U3 U4 U5 L4 L L32 L1 Fig. 4.3: Example of a monotonic wedge: data points run from left to right. appended to both U and L and we have a new (extended) monotonic wedge. Then, we would further remove L1, consider the next value forward, and so on. Algorithm 1 and Proposition 2 show that a monotonic wedge can be used to compute the - filter efficiently and with few lines of code. Algorithm 1 Streaming algorithm to compute the - filter using no more than 3 comparisons per element. 1: INPUT: an array a indexed from 1 to n 2: INPUT: window width w > 2 3: U, L← empty double-ended queues, we append to “back” 4: append 1 to U and L 5: for i in {2, ., n} do 6: if i ≥ w + 1 then 7: OUTPUT: afront(U) as maximum of range [i − w, i) 8: OUTPUT: afront(L) as minimum of range [i − w, i) 9: if ai > ai−1 then 10: pop U from back 11: while ai > aback(U) do 12: pop U from back 13: else 14: pop L from back 15: while ai < aback(L) do 16: pop L from back 17: append i to U and L 18: if i = w + front(U) then 19: pop U from front 20: else if i = w + front(L) then 21: pop L from front P 2. Algorithm 1 computes the - filter over n values using no more than 3n comparisons, or 3 comparisons per element. P. We prove by induction that in Algorithm 1, U and L form a monotonic wedge of a over the interval [max{i − w, 1}, i) at the beginning of the main loop (line 5). Initially, when i = 2, U, L = {1}, U, L is trivially a monotonic wedge. We, STREAMING MAXIMUM-MINIMUM FILTER 7 U1 U2 U STEP A3 L3 L1 L2 monotonic wedge i − w i − 1 i U1 ai > aU2 ai > ai−1 U2 U STEP B3 adding ai L3 removing U2 and U3 from U L1 L2 i − w i − 1 i U1 U2 STEP C L4 L3 L1 L2 i − w i − 1 i monotonic wedge Fig. 4.4: Algorithm 1 from line 5 to line 17: updating the monotonic wedge is done by either removing the last elements of U or the last elements of L until U, L form a monotonic wedge for [max{i − w, 1}, i]. have that the last component of both U and L is i − 1. If ai > ai−1 (line 11), then we remove the last elements of U until alast(U) ≥ an+1 (line 11) or if ai ≤ ai−1, we remove the last elements of L until alast(L) ≤ an+1 (line 15). Then we append i to both U and L (line 17). The lists U, L form a monotonic wedge of [max{i−w, 1}, i] at this point (see Fig. 4.4). After appending the latest location i (line 17), any location j < i will appear in either U or L, but not in both. Indeed, i − 1 is necessarily removed from either U or L. To compute the monotonic wedge over [max{i−w+1, 1}, i+1) from the monotonic wedge over [max{i−w, 1}, i], we check whether the location i − w is in U or L at line 18 and if so, we remove it. Hence, the algorithm produces the correct result. We still have to prove that the algorithm will not use more than 3n comparisons, no matter what the input data is. Firstly, the total number of elements that Algo- rithm 1 appends to queues U and L is 2n, as each i is appended both to U and L, 8 DANIEL LEMIRE (line 17). The comparison on line 9 is executed n − 1 time and each execution re- moves an element from either U or L (lines 10 and 14), leaving 2n−(n−1) = n+1 el- ements to be removed elsewhere. Because each time the comparisons on lines 11 and 15 gives true, an element is removed from U or L, there can only be n + 1 true comparisons. Morever, the comparisons on lines 11 and 15 can only be false once for a fixed ai since it is the exit condition of the loop. The number of false comparisons is therefore n. Hence, the total number of comparisons is at most (n − 1) + (n + 1) + n = 3n, as we claimed. While some signals such as electroencephalograms (EEG) resemble i.i.d noise, many more real-world signals are piecewise quasi-monotonic [Lemire et al. 2005]. While one G-K variant [Gil and Kimmel 2002] has a comparison complex- ity of nearly 2 comparisons per element over i.i.d noise, but a worst case complexity of slightly more than 3 comparisons for monotonic data, the opposite is true of our algorithm as demonstrated by the following proposition. P 3. When the data is monotonic, Algorithm 1 computes the - filter using no more than 2 comparisons per element. P. If the input data is non-decreasing or non-increasing, then the conditions at line 11 and line 15 will never be true. Thus, in the worse case, for each new element, there is one comparison at line 9 and one at either line 11 or line 15. The next proposition shows that the memory usage of the monotonic wedge is at most w+ 1 elements. Because U and L only store the indexes, we say that the total memory buffer size of the algorithm is 2w + O(1) elements (see Table I). P 4. In Algorithm 1, the number of elements in the monotonic wedge (size(U) + size(L)) is no more than w + 1. P. Each new element is added to both U and L at line 17, but in the next iteration of the main loop, this new element is removed from either U or L (line 10 or 14). Hence, after line 14 no element in the w possible elements can appear both in U and L. Therefore size(U) + size(L) ≤ w + 1. 5. Implementation and Experimental Results While interesting theoretically, the number of comparison per element is not neces- sarily a good indication of real-world performance. We implemented our algorithm in C++ using the STL deque template. A more efficient data structure might be possible since the size of our double-ended queues are bounded by w. We used 64 bits floating point numbers (“double” type). In the pseudocode of Algorithm 1, we append i to the two double-ended queues, and then we systematically pop one of them (see proof of proposition 2). We found it slightly faster to rewrite the code to avoid one pop and one append (see appendix). The implementation of our algo- rithm stores only the location of the extrema whereas our implementation of the H-G-W algorithm stores values. Storing locations means that we can, STREAMING MAXIMUM-MINIMUM FILTER 9 compute the arg max /min filter with no overhead, but each comparison is slightly more expensive. While our implementation uses 32 bits integers to store locations, 64 bits integers should be used when processing streams. For small window sizes, Gil and Kimmel [2002] suggests unrolling the loops, essentially compiling w in the code: in this manner we could probably do away with a dynamic data structure and the corresponding overhead. We ran our tests on an AMD Athlon 64 3200+ using a 64 bit Linux platform with 1 Gigabyte of RAM (no thrashing observed). The source code was compiled using the GNU GCC 3.4 compiler with the optimizer option “-O2”. We process synthetic data sets made of 1 million data points and report wall clock timings versus the window width (see Fig. 5.5). The linear time complexity of the naive algorithm is quite apparent for w > 10, but for small window sizes (w < 10), it remains a viable alternative. Over i.i.d. noise generated with the Unix rand function, the H-G-W and our algorithm are comparable (see Fig. 5(b)): both can process 1 million data points in about 0.15 s irrespective of the window width. For piecewise monotonic data such as a sine wave (see Fig. 5(a)) our algorithm is roughly twice as fast and can process 1 million data points in about 0.075 s. Our C++ implementation of the G-K algorithm [Gil and Kimmel 2002] performed slightly worse than the H-G-W algorithm. To insure reproducibility, the source code is available freely from the author. 6. Conclusion and Future Work We presented an algorithm to compute the - filter using no more than 3 com- parisons per element in the worst case whereas the previous best result was slightly above 3 + 2 log w/w + O(1/w) comparisons per element. Our algorithm has lower latency, is easy to implement, and has reduced memory usage. For monotonic input, our algorithm incurs a cost of no more than 2 comparisons per element. Ex- perimentally, our algorithm is especially competitive when the input is piecewise monotonic: it is twice as fast on a sine wave. We have shown that at least 2 comparisons per element are required to solve the - filter problem when no stream latency is allowed, and we showed that this bound is tight when the window is small (w = 3). Acknowledgements This work is supported by NSERC grant 261437. The author wishes to thank Owen Kaser of the University of New Brunswick for his insightful comments. References C, T.H., L, C.E., R, R.L., S., C. 2001. Introduction to Algorithms, second edition. MIT Press, Cambridge, MA. D, S.C. 1996. Running max/min calculation using a pruned ordered list. IEEE Transactions on Signal Processing 44 , 11, 2872–2877., 10 DANIEL LEMIRE 1.4 naive van Herk 1.2 Gil-Kimmel streaming 0.8 0.6 0.4 0.2 0 10 20 30 40 50 60 70 80 90 100 window size (a) Input data is a sine wave with a period of 10 000 data points. 1.2 naive van Herk 1 Gil-Kimmel streaming 0.8 0.6 0.4 0.2 0 10 20 30 40 50 60 70 80 90 100 window size (b) Input data is i.i.d. noise with a uniform distribution. Fig. 5.5: Running time to compute the - filter over a million data points using the naive algorithm, our H-G-W (van Herk) implementation, our G-K implementation, and our streaming implementation (less is better). time (s) time (s), STREAMING MAXIMUM-MINIMUM FILTER 11 D, M. B, M. J. 2005. Morphological Erosions and Openings: Fast Algo- rithms Based on Anchors. J. Math. Imaging Vis. 22 , 2-3, 121–142. G, J. K, R. 2002. Efficient Dilation, Erosion, Opening, and Closing Algorithms. IEEE Trans. Pattern Anal. Mach. Intell. 24 , 12, 1606–1617. G, J. W, M. 1993. Computing 2-D Min, Median, and Max Filters. IEEE Trans. Pattern Anal. Mach. Intell. 15, 5, 504–507. G, J (Y) K, R. 2002. Data filtering apparatus and method. US Patent Number 6,952,502. H, Y.-L., T, L.-F., Z, C.-M., C, P., L, B., M, Z.-Y. 2005. Development of intelligent diagnosis and report system based on whole body bone SPECT image. In Machine Learning and Cybernetics 2005, 5437–5441. K, E. R, C.A. 2005. Exact indexing of dynamic time warping. Knowledge and Information Systems 7, 3, 358–386. L, D., B, M., Y, Y. 2005. An Optimal Linear Time Algorithm for Quasi-Monotonic Segmentation. In ICDM’05, 709–712. M, C. Z, H. 1997. An Asynchronous Implementation of the MAXLIST Algorithm. In ICASSP’97. P, I. 1989. Fast algorithms for running ordering and max/min calculation. IEEE Transactions on Circuits and Systems 36 , 6, 795–804. S, H., L, JCS, Y, DKY. 2004. Defending against low-rate TCP attacks: dynamic detection and protection. In ICNP 2004 , 196–205. T, L. G, A. 2004. Machine vision methods and systems for boundary feature comparison of patterns and images. US Patent Number 6,687,402. H, M. 1992. A fast algorithm for local minimum and maximum filters on rectangular and octagonal kernels. Pattern Recogn. Lett. 13, 7, 517–521. Y, X., C, M., S, C. Y. 2001. A generic method of cleaning and enhancing handwritten data from business forms. International Journal on Document Analysis and Recognition 4 , 2, 84–96. Y, X., C, M., S, CY. 2001. Stroke-model-based character extraction from gray-level documentimages. IEEE Transactions on Image Processing 10 , 8, 1152–1161. Z, Y. S, D. 2003. Warping indexes with envelope transforms for query by humming. In SIGMOD’03, 181–192.

Appendix: C++ source code for the streaming algorithm

/ / input: arraya, integerwindowwidthw/ / output: arraysmaxval and minva l / / buffer: listsUand L / / requires: STLfordequesupportdeque< in t > U , L ; for( uinti= 1 ; i < a . size( ) ; ++ i ) { i f ( i>=w) { maxval [ i −w] = a [U . size( ) >0 ? U . front( ) : i −1 ] ; minval [ i −w] = a [L . size( ) >0 ? L . front( ) : i −1 ] ; } / / endifif( a [ i ] > a [ i −1] ) { L . push back ( i −1) ; i f ( i == w+L . front( ) ) L . popfront( ) ; while(U . size( ) >0) { i f ( a [ i ]<=a [U . back ( ) ] ) { i f ( i == w+U . front( ) ) U . popfront( ) ; break ; } / / endifU. pop back ( ) ; } / / endwhile, 12 DANIEL LEMIRE } else{ U . push back ( i −1) ; i f ( i == w+U . front( ) ) U . popfront( ) ; while(L . size( ) >0) { i f ( a [ i ]>=a [L . back ( ) ] ) { i f ( i == w+L . front( ) ) L . popfront( ) ; break ; } / / endifL. pop back ( ) ; } / / endwhile} / / endifelse} / / endformaxval [ a . size( )−w] = a [U . size( ) >0 ? U . front( ) : a . size( ) −1] ; minval [ a . size( )−w] = a [L . size( ) >0 ? L . front( ) : a . size( ) −1] ;]

15