Faster way to do this smoothing?

Question

Looking for a pleasing way to smooth a particular dataset I found something I like but it is rather slow to compute: I take the mean of a series of MedianFilter results computed over a range of window sizes.

SeedRandom[0];
x = RandomReal[1, 5000];

smooth = Mean @ Table[MedianFilter[x, n], {n, 500, 900, 5}];  (* slow *)

ListLinePlot[smooth]

(The closer the spacing the better the result to my eye, i.e. keeping 5 a small value.)

Is there a faster way to get the same smoothing, or something very similar?

Note:

Separate from this question is the issue of whether or not this smoothing actually represents the data in a useful or accurate way. I would rather not focus on that but instead accept this result being desirable as an assumption.

Answer 1

All smooths in this post

Here are all smooths in this post:

The line "smooth" is the one derived in the question. The lines for "smooth2" and "smooth3" were the motivation to derive the Quantile regression solution, "smoothQR". The speed-ups are around 7 and 4.5 times respectively. (See the related discussion in Rahul's answer.)

Simple modifications of question's code

Switching Mean and Median makes the computations ~ 7 times faster on my laptop. (The rationale for doing that is that finding the median is slower than finding the mean.)

SeedRandom[0];
x = RandomReal[1, 5000];

AbsoluteTiming[
 smooth = Mean@Table[MedianFilter[x, n], {n, 500, 900, 5}];
](* original code in the question *)

(* {8.44073, Null} *)

ListLinePlot[smooth]

Mean and Median switched around:

AbsoluteTiming[
 smooth2 = Median@Table[MeanFilter[x, n], {n, 500, 900, 5}];
]

(* {1.31114, Null} *)

ListLinePlot[smooth2]

Or using just Mean:

AbsoluteTiming[
 smooth3 = Mean@Table[MeanFilter[x, n], {n, 500, 900, 5}];
]

(* {1.41859, Null} *)

ListLinePlot[smooth3]

Using Quantile regression

Here is a more complicated way to speed up things together with some error measurements.

1. Select parameters.

2. Apply MedianFilter as in the question but with larger step between the min and max window sizes.

3. Find the mean or median of the obtained array. (We get a vector as a result.)

4. Apply quantile regression with B-spline basis over a sample of the points from step 3.

5. Using the found regression quantile function extract points.

6. Plot and measure errors.

The process of steps 2-5 is around 3 to 5 times faster on my laptop than the original code in the question and with a reasonable max relative error, less than 0.5%.

{a, b, step} = {500, 900, 5};
{factor, sampleSize, nKnots} = {16, 800, 25};

AbsoluteTiming[
 smooth = Mean@Table[MedianFilter[x, n], {n, a, b, s}];
]

(* {8.62565, Null} *)

Import["https://raw.githubusercontent.com/antononcube/MathematicaForPrediction/master/QuantileRegression.m"]

AbsoluteTiming[
 smoothArr = Table[MedianFilter[x, n], {n, a, b, factor*step}]; 
 ts = N@Transpose[{Rescale@Range[Length[Mean[smoothArr]]], 
     Mean[smoothArr]}]; 
 qfuncs = QuantileRegression[
   Join[ts[[1 ;; 10]], 
    RandomSample[ts[[21 ;; -11]], UpTo[sampleSize]], ts[[-10 ;; -1]]],
    nKnots, {0.5}, 
   Method -> {LinearProgramming, Method -> "InteriorPoint", 
     Tolerance -> 1.}, InterpolationOrder -> 3];
 smoothQR = 
  Table[qfuncs[[1]][y], {y, ts[[1, 1]], ts[[-1, 1]], (
    ts[[-1, 1]] - ts[[1, 1]])/(Length[x] - 1)}];
 ]

(* {1.72669, Null} *)

grOpts = {PlotRange -> All, PlotTheme -> "Detailed", 
  ImageSize -> 450}; ListLinePlot[{smooth, smoothQR}, 
 PlotLegends -> SwatchLegend[{"smooth", "smoothQR"}], grOpts]
ListPlot[Abs[smooth - smoothQR]/Abs[smooth], 
 PlotLegends -> 
  RecordsSummary[
   Abs[smooth - smoothQR]/Abs[smooth], {"Relative error"}], 
 PlotLabel -> 
  Row[{"Relative error, ", 
    HoldForm[Abs[smooth - smoothQR]/Abs[smooth]]}], grOpts]

Answer 2

One way to look at your method is that for each sample, the samples within distance 500 of it get considered in all of the medians you are taking the mean of, the samples at distance 700 get considered in half the medians, the samples at distance 900 get considered in only one, and so on.

What if we take just one median, but with each sample weighted according to the number of times it would have been considered in your original approach?

wmin = 500; 
wmax = 900; 
newSmooth = 
 Table[With[{range = Range[Max[i - wmax, 1], Min[i + wmax, Length@x]]}, 
   Median@WeightedData[x[[range]], 
     Min[#, wmax - wmin + 1] & /@ (wmax + 1 - Abs[range - i])]], 
   {i, Length@x}];

ListLinePlot[{smooth, newSmooth}, PlotStyle -> {Gray, Red}]

Pretty close, and about 10 times faster on my machine; switching to ParallelTable adds another 3x speedup. Though the result does have those telltale MedianFilter-esque discontinuities that you might not like.

(Note: The code above is designed so that if you choose wmin = wmax the result should be identical to MedianFilter[x, wmax], while if you set wmin = 0 you should get an exactly triangular weight distribution.)

---

Since Mr.Wizard has been cagey about their actual data and motivation ^{(They should really know better by now. Maybe a moderator could knock some sense into them.)}, here is an example dataset that illustrates why one might prefer a MedianFilter-based approach.

x = ImageData[ExampleData[{"TestImage", "Flower"}]][[1500, All, 1]];

Here's Mr.Wizard's smoothing method, though I changed the smoothing range to 100-500 instead of 500-900:

Show[ListPlot[x, PlotStyle -> Lighter@Gray], 
 ListLinePlot[smooth, PlotStyle -> Black]]

As you can see, MedianFilter ignores narrow spikes of outliers without blurring across significant discontinuities. Here are the rest of the smoothing methods for comparison:

Show[ListPlot[x, PlotStyle -> LightGray], 
 ListLinePlot[{smooth, smooth2, smooth3, newSmooth}, 
  PlotStyle -> {Black, Darker@Cyan, Lighter@Blue, Lighter@Red}, 
  PlotLegends -> {"Mr.Wizard", "Anton2", "Anton3", "Rahul"}]]

Answer 3

Only considering speeding it up you may use ParallelMap to spread the filtering over multiple cores.

LaunchKernels[];
smooth = ParallelMap[MedianFilter[x, #] &, Range[500, 900, 5]];
CloseKernels[];

smooth = Mean@smooth;

This executes in 13 seconds on my laptop whereas the non-parallel code takes 32 seconds.

Hope this helps.

Answer 4

Rather than an answer, this is an extended comment. All of the answers so far speed things up considerably. However, the size of the observed gyrations of the smoothed data needs to be considered relative to the scatter of points:

Show[ListPlot[x],
 ListLinePlot[smooth, PlotStyle -> {Red, Thickness[0.02]}]]

The process that generated the example dataset (which I assume is much different than your real dataset of interest) should result in a very straight horizontal line with the vertical axis value of 0.5 with there being random but serially correlated deviations from that value.

The point is that some sort of regression approach (loess, quantile regression, generalized additive models, etc.) might be considered prior to smoothing to test that any of the observed "features" might be real (rather than expected random deviations).

Answer 5

Edit: after working this up I realized I've mimicked MeanFilter, not MedianFilter

SeedRandom[0];
x = RandomReal[1, 5000];
wt = Join[#, ConstantArray[81, 1000], Reverse@#] &@
   Flatten[Table[ ConstantArray[ j, 5] , {j, 80}]];
filter[x_, i_] := Module[{wl = Round[Length@wt/2]}, Which[ 
   i < wl , (x[[;; i + wl]].#/Total@#) &@wt[[-(wl + i) ;;]],
   i > Length@x - wl , (x[[i - wl ;;]].#/Total@#) &@
    wt[[;; Length@x - wl - i]],
   True, x[[i - Length@wt/2 + 1 ;; i + wl]].wt/(Total@wt)]]

ListLinePlot[filter[x, #] & /@ Range[Length@x]]

Answer 6

This is fast, and gives a similar smoothing, though probably not as similar as you would like. It relies on the data in x being uniformly distributed though I think it could be adapted for other cases.

Consider a run of $m$ elements from $x$. If there are $m/2$ points less than 0.5 and $m/2$ greater than 0.5 then the median will be somewhere very close to 0.5. Changing one of the <0.5 points to become >0.5 will shift the median by one place, on average an increase of $1/m$. So we can approximate the median by looking at the number of points either side of 0.5. Each point <0.5 shifts the median by $-1/(2m)$ and each point >0.5 shifts the median by $+1/(2m)$.

The result can be computed quickly using ListCorrelate. It runs in a few hundredths of a second.

y = Sign[x - 0.5];

ca[n_] := ConstantArray[0.5/(2 n + 1), 2 n + 1]

smooth2 = 
  Mean @ Table[0.5 + ListCorrelate[ca[n], y, {n + 1, n + 1}, 0], {n, 500, 900, 5}];

ListLinePlot[{smooth, smooth2}]