Selection and Medians, in depth: L1 optimality, the breakdown point, MAD outlier rules, median filters, median-of-means and the geometric median

The median minimises absolute error

The mean minimises squared error: it is the c that makes the sum of (x_i - c)^2 smallest. The median minimises absolute error, the sum of |x_i - c|. The argument takes one line. Move c to the right by a small amount d. Every point to the left of c gets d further away and every point to the right gets d closer, so the total changes by d times (points left minus points right). The total stops decreasing exactly when as many points sit on each side, which is the median.

For an even count the argument says more: any c between the two middle values minimises the error, because moving inside that interval keeps the two sides balanced. With 12 values whose middle pair is 13 and 14, every c in [13, 14] is optimal, and the conventional median 13.5 is just one choice. That is why libraries disagree. NumPy's median and Python's statistics.median average the middle pair, statistics.median_low and median_high return one of them, and a selection routine returns whichever rank you asked for. Pick one definition per system and write it down.

The L1 view matters in practice. Fitting a model with absolute loss estimates a conditional median, which is why it resists outliers. Quantile (pinball) loss generalises it: weighting errors above by q and below by 1 - q estimates the q-th quantile, the same balancing argument with unequal weights.

Breakdown point: why one bad value cannot move it

The breakdown point of an estimator is the smallest fraction of the data an adversary must corrupt to drive the estimate arbitrarily far. For the mean it is 1/n: one value set to a million moves the mean as far as you like. For the median it is about one half: until more than half the points are corrupted, the median stays between two honest values. No location estimator that treats points symmetrically can do better than 50 percent, because past that point the corrupted data could itself be the honest majority.

The same holds for spread. The standard deviation breaks with one bad point. The median absolute deviation, MAD = median(|x_i - median(x)|), also has a 50 percent breakdown point. For normal data, MAD multiplied by 1.4826 (that is, divided by 0.6745) estimates the standard deviation, so the two can be swapped where the data is clean and the MAD keeps working where it is not.

The price is efficiency. On perfectly normal data the sample median has higher variance than the mean, roughly pi / 2, about 1.57 times, for large n. You pay that insurance premium on every clean dataset to be safe on the dirty ones.

Worked example: two slow requests and a blind alert

Thirteen request latencies in milliseconds arrive at an alerting job: eleven between 11 and 16, and two from a stalled dependency at 900 and 950. A common rule flags points more than three standard deviations from the mean.

import numpy as np

lat = [12, 14, 13, 15, 11, 13, 14, 12, 950, 13, 16, 14, 900]
m, s = np.mean(lat), np.std(lat)
print(m, s)                                    # 153.6, 329.1
print([v for v in lat if abs(v - m) > 3 * s])  # []  -- nothing flagged

def mad_outliers(x, z=3.5):
    # Modified z-score of Iglewicz and Hoaglin: 0.6745 * (x - median) / MAD
    x = np.asarray(x, dtype=float)
    med = np.median(x)
    mad = np.median(np.abs(x - med))
    if mad == 0:                     # over half the values identical: no spread
        return np.zeros(len(x), dtype=bool), med, mad
    return np.abs(0.6745 * (x - med) / mad) > z, med, mad

flags, med, mad = mad_outliers(lat)
print(med, mad, np.array(lat)[flags])           # 14.0 1.0 [950 900]

The three-sigma rule flags nothing. The two outliers inflated the standard deviation to 329 ms, so each sits only 2.27 and 2.42 standard deviations from a mean of 153.6 that describes no actual request. This is masking: outliers hide themselves by corrupting the yardstick. The median of 14 and MAD of 1 are untouched, and both outliers score above 500 on the modified z-score against the usual 3.5 cut-off. With only one of the two outliers the three-sigma rule happens to fire (z = 3.32), which is exactly why the masking failure survives testing: it appears only when trouble arrives in groups.

Note the guard for MAD = 0. When more than half the values are identical, such as a counter that is usually zero, the MAD is zero and every other value would score as infinite. Fall back to a percentile-based spread or a fixed floor in that case rather than paging on every non-zero sample.

Median filters: remove spikes, keep edges

A median filter replaces each sample with the median of a window around it. Unlike a moving average, it removes isolated spikes completely instead of smearing them, and it preserves step edges, because the median of a window that straddles a step is one of the two levels, not something in between. That makes it the standard first stage for impulse (salt and pepper) noise in images and glitch removal in sensor data.

def median_filter_1d(x, w):
    # Odd window w, edges padded by reflection
    h = w // 2
    xp = np.pad(np.asarray(x, dtype=float), h, mode="reflect")
    win = np.lib.stride_tricks.sliding_window_view(xp, w)
    return np.median(win, axis=1)

sig = [1, 1, 1, 9, 1, 1, 5, 5, 5, 5, 0, 5, 5]
print(median_filter_1d(sig, 3))   # [1 1 1 1 1 1 5 5 5 5 5 5 5]

The spike to 9 and the dropout to 0 vanish and the step from 1 to 5 stays exactly where it was. A window of width w removes runs of up to (w - 1) / 2 bad samples; a run longer than that is treated as signal. Choose w from the longest glitch you need to remove, not from smoothness.

Edges need a policy. scipy.signal.medfilt zero-pads, which pulls the first and last samples toward 0 on a signal that does not live near 0; scipy.ndimage.median_filter reflects by default. On this example both give the same output, but on a signal sitting at 100 the zero-padded version distorts the ends. For images, the naive cost is a selection per pixel over r squared values. Huang's 1979 algorithm keeps a histogram of the window and updates it by one column as the window slides, which is fast for 8-bit images, and Perreault and Hebert's 2007 method makes the per-pixel cost constant in the radius. For streaming one-dimensional windows, the moving percentile article compares sorted windows, two heaps with lazy deletion and Fenwick trees.

Median-of-means for heavy tails

Sometimes you need the mean, not the median, for example the expected cost per request, but the data is heavy-tailed. The sample mean is unbiased, yet its deviations have heavy tails too: a rare huge value moves it a long way. Median-of-means splits the data into k groups, averages each, and returns the median of the k averages. With only a finite variance assumed, it achieves deviation bounds close to what the mean gets under Gaussian data, with the confidence level set by k.

def median_of_means(x, k, rng):
    x = rng.permutation(np.asarray(x, dtype=float))   # shuffle so groups are exchangeable
    return float(np.median([g.mean() for g in np.array_split(x, k)]))

rng = np.random.default_rng(42)
true = 2.1 / 1.1                     # mean of a Pareto with shape 2.1 and scale 1
for _ in range(2000):
    x = rng.pareto(2.1, 1000) + 1
    # record abs(x.mean() - true) and abs(median_of_means(x, 20, rng) - true)

The table shows the trade honestly. Median-of-means cuts the 99th percentile error by a third and the worst case by two thirds, but its typical error is worse. On skewed data each group mean is itself right-skewed, so the median of group means sits below the true mean: the estimator buys tail safety with bias. Use it when the occasional large miss is the expensive outcome, such as a budget or capacity estimate, and keep the plain mean when you average over many runs anyway.

Medians in more than one dimension

For vectors there is no single ordering, so there is no single median. The coordinate-wise median takes the median of each coordinate separately. It is cheap but depends on the coordinate system, and it can land outside the data altogether: for the three points (1,0,0), (0,1,0) and (0,0,1) it returns (0,0,0), which is not in their convex hull.

The geometric median minimises the sum of Euclidean distances, the direct analogue of the one-dimensional L1 property. It is rotation-equivariant, has a 50 percent breakdown point, and has no closed form. Weiszfeld's 1937 iteration computes it as a repeatedly reweighted average, each point weighted by the inverse of its distance to the current estimate:

def weiszfeld(points, iters=200, tol=1e-9):
    y = points.mean(axis=0)
    for _ in range(iters):
        d = np.linalg.norm(points - y, axis=1)
        if np.any(d < 1e-12):       # simplification: stop on a data point
            return y
        w = 1.0 / d
        y_new = (points * w[:, None]).sum(axis=0) / w.sum()
        if np.linalg.norm(y_new - y) < tol:
            return y_new
        y = y_new
    return y

pts = np.array([[0, 0], [1, 0], [0, 1], [1, 1], [50, 50]], float)
# centroid (10.4, 10.4); coordinate-wise (1, 1); geometric median (0.789, 0.789)

For the unit square plus one point at (50, 50) the centroid is dragged to (10.4, 10.4), the coordinate-wise median snaps to the corner (1, 1), and the geometric median stays inside the square at (0.789, 0.789). Stopping when the iterate hits a data point is a simplification, since that point may not be the optimum; the Vardi and Zhang modification handles the case correctly. Geometric medians are used to aggregate model updates robustly in federated learning, to place facilities, and as a robust centre in clustering (k-medians).

Operational guidance

Medians do not merge. You cannot combine shard medians into a global median. Values 1 to 3 on one shard and 4 to 10 on another have shard medians 2 and 7; their average is 4.5, but the true median of 1 to 10 is 5.5. Distributed systems either ship all the data (or a sample) to one place, run a distributed selection over several rounds as described in the selection algorithms article, or keep a mergeable quantile sketch such as a t-digest and accept an approximation. Exact one-pass medians need memory linear in n, a lower bound due to Munro and Paterson.

Compute with selection, not sorting. numpy.partition(a, k) and C++ std::nth_element run in average linear time and place the k-th value at index k; np.partition also accepts a list of indices to fetch both middle values at once. Sort only when you need many quantiles of the same array.

Decide NaN handling explicitly. numpy.median returns nan if any value is nan; numpy.nanmedian skips them; pandas skips them by default. A dashboard whose median silently ignores failed requests that were logged as NaN will look healthy during an outage, so count missing values beside the median.

Failure modes

• Averaging percentiles across hosts. The mean of per-host medians is not the fleet median, and the error grows with load imbalance. Merge sketches instead.
• Median of a multimodal distribution. With two clusters, the median can sit in the empty gap between them and describe no real case. Plot a histogram before trusting any single centre.
• Zero MAD. Constant or mostly constant data makes MAD zero and every deviation infinite. Guard it.
• Window too short for the glitch. A median filter of width w passes runs longer than (w - 1) / 2 unchanged, so a three-sample outage survives a width-5 filter.
• Medians hide the tail. A stable median can coexist with a degraded 99th percentile. Alert on tail quantiles too.

Trade-offs

What to do next

1. Find one alert that uses mean plus k standard deviations and replay last month's data through the MAD rule above; count the masked incidents.
2. Write down the median definition (low, high or average of the middle pair) and the NaN policy for every service metric you own.
3. Replace any average of per-host percentiles with a mergeable sketch.
4. Before smoothing a noisy signal, measure the longest glitch and size a median filter to it; compare with the moving average you have now.
5. If you estimate costs from heavy-tailed data, run the median-of-means experiment above on your own data and compare the tail errors.
6. Learn the selection engine underneath: read the k-th order statistics article for quickselect and multiselect.