# How can I generate a $2n$-point moving maximum of a list as efficiently as MaxFilter?

This is based on the prior question, Moving maximum function?, where @rasher provided two winning solutions (i.e. clearly the most efficient):

MaxFilter[list, 1]


and

Max /@ Transpose[{Rest[Append[#, 0]], #, Most[Prepend[#, 0]]}] &[list]


Taking the same example list,

{5, 6, 9, 3, 2, 6, 7, 8, 1, 1, 4, 7}


I'd like to be able to specify a window $w$ such that a rolling maximum starts at index $w$, generating, for example for $w=5$ (lined up on purpose):

            {9, 9, 9, 8, 8, 8, 8, 8}


That's not too hard; as @TylerDurden hinted, we can simply use MaxFilter with range $(w-1)/2=2$ and drop each end ("shifting," as he imagined):

w=5;
Drop[Drop[MaxFilter[list, (w - 1)/2], (w - 1)/2], -(w - 1)/2]
(* {9, 9, 9, 8, 8, 8, 8, 8} *)


But of course, this can't be done with an even number $w$.

So, giving up, I decided to try to generalize @rasher's second solution—i.e. staggering $w$ copies of the list (in his case, essentially hardcoded to $w=3$) and getting maximums across the transpose:

AnotherMaxFilter[list_, w_] :=
Max /@ (
Array[
ConstantArray[0, w - #]~Join~list &,
w
]~Flatten~{2}
)[[w~Range~Length@list, Range@w]]


Note: I'm using Flatten to do a jagged transpose.

However, my attempt wasn't efficient at all:

list = RandomInteger[10, 10^7];

AbsoluteTiming[AnotherMaxFilter[list, 3]][[1]]
(* 5.257301 *)

AbsoluteTiming[AnotherMaxFilter[list, 5]][[1]]
(* 7.153409 *)

AbsoluteTiming[AnotherMaxFilter[list, 7]][[1]]
(* 14.786846 *)


Is there a way to improve my generalization to reach @rasher's original efficiency? Or, is there a clever way to get MaxFilter to max over even-numbered windows?

I'd also very much appreciate if someone could explain what about my code could have introduced such inefficiency.

• The example output of  {9, 9, 9, 8, 8, 8, 8, 8} does not match the output of AnotherMaxFilter for window 5 for your example list {5, 6, 9, 3, 2, 6, 7, 8, 1, 1, 4, 7}... which one is the desired result? – ciao Mar 25 '15 at 5:45
• It appears I'd written a forward-looking moving maximum. Now fixed. Thanks, @rasher. – Andrew Cheong Mar 25 '15 at 5:57
• Based on what you've done, mfw = (Max /@ Partition[#1, #2, 1]) &;, mfw[<list>,<window>] should do the same more efficiently. I'll ponder a silver bullet... – ciao Mar 25 '15 at 6:01
• @rasher - Thanks, indeed your solution gets much, much closer to MaxFilter's efficiency (without any compilation, which is impressive to me). – Andrew Cheong Mar 26 '15 at 4:05

I found a clever way to make MaxFilter work with even windows, and used conditionals to combine the odd- and even- cases. This solution is, surprisingly, more efficient than even @MrWizard's compiled function, and increasingly so as the window size increases.

# tl;dr

MovingMax[list_, w_] := If[w > Length@list, {}, Module[{r, tmp},
r = Floor[w/2];
If[
OddQ@w
,
MaxFilter[list, r][[r + 1 ;; -(r + 1)]]
,
(* MaxFilter only supports odd windows; here's a hack for even windows. *)
tmp = MaxFilter[Riffle[list, Min@list, w], r][[r + 1 ;; -(r + 1)]];
If[
w > Length@tmp + 1
,
(* If window is greater than number of resulting elements, drop nothing. *)
tmp
,
Drop[tmp, {w, Length@tmp, w}]]
]
]
];


To create MovingMin, it will work to just swap all references to Max and Min with the opposites.

# How It Works

If the window size $w$ is odd, then we can just use MaxFilter with a radius of $(w-1)/2$ and drop the same amount of elements (i.e. as the radius) from each end of the resulting list.

If the window size $w$ is even, my idea was to insert a $-\infty$ for every $w$th element, e.g. turning

{1, 2, 3, 4, 5, 6, 7}


into

{1, 2, 3, -∞, 4, 5, 6, -∞, 7}


if $w$ were $4$. Because then, MaxFilter would essentially be finding the maximums of...

         1,  2,  3    = 3     // incomplete list; drop
1,  2,  3, -∞    = 3     // incomplete list; drop
1,  2,  3, -∞,  4    = 4
2,  3, -∞,  4,  5    = 5
3, -∞,  4,  5,  6    = 6
-∞,  4,  5,  6, -∞    = 6     // incomplete list; drop
4,  5,  6, -∞,  7    = 7
5,  6, -∞,  7        = 7     // incomplete list; drop
6, -∞,  7            = 7     // incomplete list; drop


...so after dropping the first and last two elements (same as we do for the odd case), we just have to additionally drop every $w$th element, and we're left with the maximums we're looking for.

Later, I replaced $-\infty$ with Min@list as the former method was adding an inefficiency (see @MrWizard's comment).

# Efficiency

Testing on a list of 20,000,000 random integers and varying window sizes, here were the results:

• @MrWizard's cf performs well, but increases in time linearly with the window size.

• @rasher's Partition-based method ran my machine out of memory for larger window sizes. (But @rasher meant to improve on my inefficient AnotherMaxFilter, not submit a serious competitor, so it's unfair to compare his solution; I only included it out of curiosity.)

• The MaxFilter functions (split between odd and even since they're essentially two unrelated functions), in contrast, decrease in time consumption as window size increases! I imagine this is because MaxFilter probably optimizes by caching the index of its latest maximum as it runs, e.g.

 6 8 12 14 6 9 11 7 13 17 3 9 20 20 12 18 18 1 3 16
|       |
+-------+


where the first maximum is 14 and the cached index is 4. This way, as it moves forward, it only needs to compare one number, e.g.

 6 8 12 14 6 9 11 7 13 17 3 9 20 20 12 18 18 1 3 16
|       |
+-------+


where since 6 is not greater than the current maximum of 14, there's no need to compare the new set of 4 elements. This can continue until the cached index "expires," e.g.

 6 8 12 14 6 9 11 7 13 17 3 9 20 20 12 18 18 1 3 16
|       |
+-------+


when a new maximum must be calculated (e.g. 11, now caching index 7.)

With such an algorithm, greater window sizes would mean greater savings.

# Edits and Follow-Up

### Original Code

MovingMax[list_, w_, lowerBound_: - Infinity] := Module[{r, tmp},
r = Floor[w/2];
If[
OddQ@w
,
Drop[Drop[MaxFilter[list, r], r], -r]
,
tmp = Drop[Drop[MaxFilter[Riffle[list, lowerBound, w], r], r], -r];
Drop[tmp, {w, Length@tmp, w}]
]
];


This code is much more efficient if given an actual lower bound, e.g. -10^-6, rather than using -∞. @MrWizard explains the reason in the comments.

### Original Test Code

Here's the (unedited, messy) test code I used.

FCompiled =
Compile[{{x, _Integer, 1}, {n, _Integer}},
Module[{i = n,
a = Take[x, n]}, (a[[Mod[i, n, 1]]] = #; i++; Max[a]) & /@
Drop[x, n - 1]]];
FPartition = (Max /@ Partition[#1, #2, 1]) &;
FMaxFilterOdd =
If[OddQ@#2,
Drop[Drop[MaxFilter[#, (#2 - 1)/2], (#2 - 1)/2], -(#2 - 1)/2]] &;
FMaxFilterEven =
If[EvenQ@#2,
Module[{tmp},
Drop[tmp =
Drop[Drop[
MaxFilter[Riffle[#, -Infinity, #2], #2/2], #2/2], -#2/
2], {#2, Length@tmp, #2}]]] &;

list = RandomInteger[{-1*^6, 1*^6}, 100000];
timings = results = ConstantArray[Null, {4, 8}];

test[fn_, row_] := (
timings[[row, 1]] = ((results[[row, 1]] = fn[list, 3]); // AbsoluteTiming // First);
timings[[row, 2]] = ((results[[row, 2]] = fn[list, 6]); // AbsoluteTiming // First);
timings[[row, 3]] = ((results[[row, 3]] = fn[list, 9]); // AbsoluteTiming // First);
timings[[row, 4]] = ((results[[row, 4]] = fn[list, 12]); // AbsoluteTiming // First);
If[
! SameQ[fn, FPartition]
,
timings[[row, 5]] = ((results[[row, 5]] = fn[list, 301]); // AbsoluteTiming // First);
timings[[row, 6]] = ((results[[row, 6]] = fn[list, 602]); // AbsoluteTiming // First);
timings[[row, 7]] = ((results[[row, 7]] = fn[list, 903]); // AbsoluteTiming // First);
timings[[row, 8]] = ((results[[row, 8]] = fn[list, 1204]); // AbsoluteTiming // First);
]
);

test[FCompiled, 1];
test[FPartition, 2];
test[FMaxFilterOdd, 3];
test[FMaxFilterEven, 4];

TableForm[timings // Transpose,
TableHeadings -> {{3, 6, 9, 12, 301, 602, 903, 1204}, {"@MrWizard's",
"@rasher's", "MF-based (odd)", "MF-based (even)"}}]


### Analysis of New Version

@MrWizard suggested two improvements: replacing Drop[Drop[..., r, -r]] with [[r+1;;-r-1]], and using Min@list as the lower bound. I ran some new tests:

The first improvement's impact was clear. The second improvement would of course slow it down (since computing a Min can't be faster than a pre-supplied lower bound), but it was worth eliminating a parameter. Therefore both improvements have been incorporated into the new "tl;dr."

• Very nice! I tried to implement the caching thing you hypothesize for MaxFilter but it ended up slower than the code I posted. Perhaps I should try again. By the way I think you can replace the Drop[Drop[ thing with MaxFilter[list, r][[r + 1 ;; -r - 1]] for a slight improvement. – Mr.Wizard Mar 26 '15 at 9:57
• Regarding lowerBound the use of -∞ will force Riffle to unpack. I believe you can replace lowerBound with Min[list]. – Mr.Wizard Mar 26 '15 at 10:01
• Incorporated your suggestions into a new tl;dr (after a new analysis, edited into the end of the post), and also fixed the function for corner cases like window size being greater than the list size. If you have any other suggestions / simplifications / improvements, including style, please let me know. – Andrew Cheong Mar 26 '15 at 12:49

Compilation is useful here (with a method from Efficient circular buffer?):

cf =
Compile[{{x, _Integer, 1}, {n, _Integer}},
Module[{i = n, a = Take[x, n]},
(a[[ Mod[i, n, 1] ]] = #; i++; Max[a]) & /@ Drop[x, n - 1]
]
];


Test:

x = RandomInteger[{-1*^6, 1*^6}, 200000];

(r1 = AnotherMaxFilter[x, 200]); // AbsoluteTiming // First
(r2 = cf[x, 200]);               // AbsoluteTiming // First

r1 === r2

6.364364

0.039002

True


We get more than two orders of magnitude improvement in this example. Compiling to C could likely be faster still, but note that the InputForm of cf indicates that there remain callbacks as it contains "Take", "Drop" and "MaxIT". These would probably need to be removed for optimal C generation. Perhaps:

cf2 =
Compile[{{x, _Integer, 1}, {n, _Integer}},
Module[{i = n, a = x[[1 ;; n]]},
(a[[ Mod[i, n, 1] ]] = #; i++; Fold[If[#2 > #, #2, #] &, a[[1]], a[[2 ;; n]]]) & /@
x[[ n - 1 ;; Length[x] ]]
],
CompilationTarget -> "C", RuntimeOptions -> "Speed"
];


However I do not have a C compiler installed to test this.

• surprisingly on your example it is slightly slower... – chris Mar 25 '15 at 9:20
• @chris What is? cf2? Something else? – Mr.Wizard Mar 25 '15 at 9:21
• cf = Compile[{{x, _Integer, 1}, {n, _Integer}}, Module[{i = n, a = Take[x, n]}, (a[[Mod[i++, n, 1]]] = #; Max[a]) & /@ Drop[x, n - 1]]]; cf2 = Compile[{{x, _Integer, 1}, {n, _Integer}}, Module[{i = n, a = x[[1 ;; n]]}, (a[[Mod[i++, n, 1]]] = #; Fold[If[#2 > #, #2, #] &, a[[1]], a[[2 ;; n]]]) & /@ x[[n - 1 ;; Length[x]]]], CompilationTarget -> "C", RuntimeOptions -> "Speed"]; x = RandomInteger[{-1*^6, 1*^6}, 200000]; (r = cf[x, 200]) // AbsoluteTiming // First (r2 = cf2[x, 200]) // AbsoluteTiming // First Out[11]= 0.096349 Out[12]= 0.105093 – chris Mar 25 '15 at 9:22
• Hmm., I get CompiledFunction::cflist: Nontensor object generated; proceeding with uncompiled evaluation trying this. Did something change in compilability between 9.x and 10.x? – ciao Mar 25 '15 at 23:17
• @Mr.Wizard: Ah, looks like a subtle bug in compiler perhaps - changing to (a[[Mod[i, n, 1]]] = #; i++; Max[a]) works fine on 9.x, but (a[[Mod[i++, n, 1]]] = #; Max[a]) poops the error. Strange and surprising... – ciao Mar 25 '15 at 23:52

Here's a version allowing MaxFilter to work with windows of even length k. It runs MaxFilter with window radius k/2-1, then corrects the output. MovingMaxEven is slower than Andrew's MovingMax, for example, 1.01 s versus his 0.78 s on 10 million points.

MovingMaxEven[s_List, k_?EvenQ] :=
Block[{r = k/2 - 1, f},
f = MaxFilter[s, r][[r + 1 ;; -r - 2]];
f + (UnitStep[#]*#)&[s[[k ;;]] - f]
]


The compiled version runs in 0.90 s.

MovingMaxEvenC = Compile[{{s, _Integer, 1}, {k, _Integer}},
Block[{r = 0, f},
r = Quotient[k, 2] - 1;
f = MaxFilter[s, r][[r + 1 ;; -r - 2]];
f + (UnitStep[#]*#)&[s[[k ;; Length[s]]] - f]
], CompilationTarget -> "C", RuntimeOptions -> "Speed"];