An Introduction to the Median

Chinese

Reprinted from Jianlin Su’s article “An Introduction to the Median“. The annotations in this post were added by me.

I recently revisited the concept of the median, so I decided to write down the key points while they are still fresh.

When removing or clipping outliers, we often need a “baseline.” For example, for a batch of nonnegative data, we may regard any value greater than 50 times the baseline as an outlier. So how should that baseline be chosen? A common choice is the mean, but the mean is easily pulled away by outliers. As a result, using it as the baseline may bias us toward those outliers and cause some problematic results to be missed. In such cases, the median can be a better baseline.

Basic Properties

For one-dimensional data points \(x_1, x_2, \cdots, x_n\), their mean is defined as:

$$ \text{mean}(x_1, x_2, \cdots, x_n) = \frac{1}{n} \sum_{i=1}^n x_i \quad (1) $$

Because every data point directly participates in the averaging process, a few especially large values can drive the mean upward and interfere with outlier detection.

The idea of the median is to find a “neutral” point in the data and use it as the baseline. More specifically, it seeks a dividing point such that the number of data points greater than or equal to it is the same as the number less than or equal to it. For this reason, it is also called the “50% quantile.” If \(n\) is odd, it is the \((n+1)/2\)-th largest value in the dataset. If \(n\) is even, then any number between the \(n/2\)-th and the \((n/2+1)\)-th largest values can be regarded as a median, though in practice we usually take their average.

Its robustness is already clear from this definition: as long as the relative order of each value with respect to the median does not change, the median itself will not change. More generally, we can describe this robustness using the breakdown point, defined as the smallest proportion of outliers required to drive the estimator to infinity. For the mean, the answer is essentially 0, because one point going to infinity is enough to make the mean go to infinity. For the median, however, the answer is 50%, because more than half of the values must go to infinity before the median does.

The downside of the median is that it is more complicated to compute than the mean, because it requires some degree of sorting. This is especially inconvenient in distributed-computing scenarios. Fortunately, most applications do not require an extremely precise median, so there is room for compromise, for example by computing local medians first and then computing a global mean or median from them, which can significantly reduce communication overhead.

Solving for the median does indeed have a larger constant factor than simply summing values for a mean. But in the strict sense of algorithmic complexity, finding the median absolutely does not require a full sort with time complexity \(O(n \log n)\). We can use Quickselect to find the median exactly in expected time \(O(n)\). Quickselect can also be extended to distributed settings.

Its core idea comes from the same source as Quicksort: choose a “pivot” element and partition the array into elements smaller than it and larger than it. At that point, the pivot is already in its absolutely correct sorted position. Finding the median is essentially finding the element ranked \(n/2\), so we only need to compare the pivot index with the target index, discard the half that definitely cannot contain the median, and recurse on the remaining half. The total amount of computation is: \(n + n/2 + n/4 + n/8 + \cdots\). By the geometric-series sum formula, this limit converges strictly to \(2n\). Therefore, its expected time complexity is linear, namely \(O(n)\).

Optimization Perspective

Interestingly, from an optimization perspective, the mean and the median can be unified:

$$ \text{mean}(x_1, x_2, \cdots, x_n) = \operatorname*{argmin}_{\mu} \sum_{i=1}^n (x_i - \mu)^2 \quad (2) $$

$$ \text{median}(x_1, x_2, \cdots, x_n) = \operatorname*{argmin}_{\mu} \sum_{i=1}^n |x_i - \mu| \quad (3) $$

The proof for \(\text{mean}\) is straightforward, so I will leave it to the reader:

Let \(f(\mu) = \sum_ {i=1}^ {n} (x_ {i} - \mu)^ {2}\):

$$ f'(\mu) = -2\sum_{i=1}^n (x_i - \mu) = 0 $$

Expanding and rearranging gives:

$$ \sum_{i=1}^n x_i = n\mu \quad \Rightarrow \quad \mu = \frac{1}{n}\sum_{i=1}^n x_i $$

Since the second derivative is always positive, this point is the global minimum.

Here I will briefly introduce the proof for \(\text{median}\). Without loss of generality, assume the \(x_i\) have already been sorted, that is, \(x_1 \le x_2 \le \cdots \le x_n\). Let \(f(\mu) = \sum_{i=1}^n |x_i - \mu|\). Taking the derivative directly gives:

$$ f'(\mu) = \sum_{i=1}^n \text{sign}(\mu - x_i) = \#\{x_i < \mu\} - \#\{x_i > \mu\} \quad (4) $$

Here # denotes the number of values satisfying the condition. To find the minimizer, we want the derivative to be as close to 0 as possible, that is, we want the number of \(x_i\) greater than \(\mu\) and less than \(\mu\) to be as balanced as possible, which already matches the intuition behind the median. A more concrete discussion splits into cases:

$$ f'(\mu) = \begin{cases} \left. \begin{aligned} \underbrace{\#\{x_i < \mu\}}_{\le k} - \underbrace{\#\{x_i > \mu\}}_{\ge k+1} < 0, & \quad \mu < x_{k+1} \\ \underbrace{\#\{x_i < \mu\}}_{\ge k+1} - \underbrace{\#\{x_i > \mu\}}_{\le k} > 0, & \quad \mu > x_{k+1} \end{aligned} \right\} & n = 2k + 1 \\ \left. \begin{aligned} \underbrace{\#\{x_i < \mu\}}_{\le k-1} - \underbrace{\#\{x_i > \mu\}}_{\ge k+1} < 0, & \quad \mu < x_k \\ \underbrace{\#\{x_i < \mu\}}_{\ge k+1} - \underbrace{\#\{x_i > \mu\}}_{\le k-1} > 0, & \quad \mu > x_{k+1} \end{aligned} \right\} & n = 2k \end{cases} \quad (5) $$

Therefore, when \(n\) is odd, \(n = 2k + 1\), both \(\mu < x_ {k+1}\) and \(\mu > x_ {k+1}\) increase \(f(\mu)\), so \(\mu^ {\ast} = x_ {k+1}\). When \(n\) is even, \(n = 2k\), both \(\mu < x_ {k}\) and \(\mu > x_ {k+1}\) increase \(f(\mu)\), so \(\mu^ {\ast} \in [x_ {k}, x_ {k+1}]\). One can verify that when \(\mu^ {\ast}\) lies in this interval, \(f(\mu^ {\ast})\) remains unchanged, so every value in the interval is a valid \(\mu^ {\ast}\). In summary, \(\mu^ {\ast}\) matches the definition of the median exactly.

Higher-Dimensional Space

From the optimization perspective, we can also understand why the median is more resistant to outliers than the mean: if one \(x_i\) is extremely large, then the mean incurs a loss of \((x_i - \mu)^2\), whereas the median incurs only \(|x_i - \mu|\). In general, \((x_i - \mu)^2 \gg |x_i - \mu|\), so the mean is more strongly pulled toward the outlier in order to reduce the loss.

Another advantage of the optimization perspective is that it extends naturally to higher-dimensional spaces. We know that the notion of the mean generalizes easily: for a set of vectors \(\boldsymbol{x}_ {1}, \boldsymbol{x}_ {2}, \cdots, \boldsymbol{x}_ {n}\), the mean vector is simply \(\frac{1}{n} \sum_ {i=1}^ {n} \boldsymbol{x}_ {i}\). The notion of the median, however, does not generalize directly, because it relies on sorting, and there is no natural total order for vector-valued data.

Under the optimization perspective, however, the extension is completely natural:

$$ \text{mean}(\boldsymbol{x}_1, \boldsymbol{x}_2, \cdots, \boldsymbol{x}_n) = \operatorname*{argmin}_{\boldsymbol{\mu}} \sum_{i=1}^n \|\boldsymbol{x}_i - \boldsymbol{\mu}\|_2^2 \quad (6) $$

$$ \text{median}(\boldsymbol{x}_1, \boldsymbol{x}_2, \cdots, \boldsymbol{x}_n) = \operatorname*{argmin}_{\boldsymbol{\mu}} \sum_{i=1}^n \|\boldsymbol{x}_i - \boldsymbol{\mu}\|_2 \quad (7) $$

Here \(\lVert \cdot \rVert_ {2}\) denotes the Euclidean distance. It is easy to prove that the mean vector defined this way is exactly \(\frac{1}{n} \sum_ {i=1}^ {n} \boldsymbol{x}_ {i}\), which agrees with the empirical definition. As for the median vector, it is also called the geometric median. From the objective function, we can see that if \(n=3\), it is actually the classical Fermat point, so in many cases people also refer to the median vector directly as the Fermat point.

Unfortunately, the median vector has no closed-form solution. It is usually computed by the following Weiszfeld iteration:

$$ \boldsymbol{\mu}_{t+1} = \frac{\sum_{i=1}^n \boldsymbol{x}_i / \|\boldsymbol{x}_i - \boldsymbol{\mu}_t\|_2}{\sum_{i=1}^n 1 / \|\boldsymbol{x}_i - \boldsymbol{\mu}_t\|_2} \quad (8) $$

Further Generalization

Clearly, we can also consider a more general extension:

$$ \text{average}(\boldsymbol{x}_1, \boldsymbol{x}_2, \cdots, \boldsymbol{x}_n; \alpha) = \operatorname*{argmin}_{\boldsymbol{\mu}} \sum_{i=1}^n \|\boldsymbol{x}_i - \boldsymbol{\mu}\|_2^\alpha \quad (9) $$

In statistical learning, “minimizing some notion of distance (error),” “assuming that the data follow some probability distribution,” and “training a machine-learning model with only a constant term to predict the data” are completely equivalent. The bridge connecting them is maximum likelihood estimation (MLE).

Suppose we have a set of observed data \(X = {x_1, x_2, \cdots, x_n}\), and we want to find a center point \(\mu\). We can build the following observation model:

$$x_i = \mu + \epsilon_i$$

where \(\epsilon_i\) is the random noise (error) of the \(i\)-th observation.

Assume that all noises \(\epsilon_i\) are independent and identically distributed (i.i.d.), and that their probability density function (PDF) is \(f(\epsilon)\). Then the probability density of observing a single data point \(x_i\) can be written as \(f(x_i - \mu)\). For the entire dataset \(X\), the joint probability, that is, the likelihood function, is the product of the probabilities of all observations:

$$L(\mu | X) = \prod_{i=1}^n f(x_i - \mu)$$

We need to find a \(\mu\) that maximizes the probability of observing the current dataset \(X\).

For convenience, we usually take the natural logarithm of the likelihood function and add a minus sign, thereby turning the problem into minimizing the negative log-likelihood (NLL):

$$\text{NLL}(\mu) = -\ln L(\mu | X) = -\sum_{i=1}^n \ln f(x_i - \mu)$$

Now recall the generalized distance objective given in the article, using the one-dimensional case as an example:

$$J(\mu) = \sum_{i=1}^n |x_i - \mu|^\alpha$$

For MLE to be equivalent to this objective, we must have:

$$-\sum_{i=1}^n \ln f(x_i - \mu) \propto \sum_{i=1}^n |x_i - \mu|^\alpha$$

This means that for a single error term \(\epsilon = x_i - \mu\), we must have:

$$f(\epsilon) \propto \exp(-c |\epsilon|^\alpha)$$

This form of distribution already has a standard name in statistics: the generalized normal distribution. Its fully normalized form is:

$$f(\epsilon; \alpha, \beta) = \frac{\alpha}{2\beta \Gamma(1/\alpha)} \exp\left( -\left( \frac{|\epsilon|}{\beta} \right)^\alpha \right)$$

If we substitute this \(f(\epsilon)\) back into the negative log-likelihood formula, we obtain:

$$\text{NLL}(\mu) = \sum_{i=1}^n \left( \left( \frac{|x_i - \mu|}{\beta} \right)^\alpha - \ln \left( \frac{\alpha}{2\beta \Gamma(1/\alpha)} \right) \right)$$

Since \(\alpha\) and \(\beta\) are constants when solving for \(\mu\), they disappear under differentiation. Therefore, minimizing the expression above is completely equivalent to:

$$\operatorname*{argmin}_{\mu} \sum_{i=1}^n |x_i - \mu|^\alpha$$

In this optimization problem, adjusting the exponent \(\alpha\) amounts to imposing a physical model assumption on the data-generating process in the real world: if we assume a Gaussian distribution, we set \(\alpha = 2\); if we assume a Laplace distribution, we set \(\alpha = 1\).

Let \(f(\boldsymbol{\mu}) = \sum_{i=1}^n \lVert \boldsymbol{x}_ {i} - \boldsymbol{\mu} \rVert_ {2}^ {\alpha}\). Then:

$$ \nabla_{\boldsymbol{\mu}} f(\boldsymbol{\mu}) = \alpha \sum_{i=1}^n \|\boldsymbol{x}_i - \boldsymbol{\mu}\|_2^{\alpha-2}(\boldsymbol{\mu} - \boldsymbol{x}_i) \quad (10) $$

Setting \(\nabla_{\boldsymbol{\mu}} f(\boldsymbol{\mu}) = \boldsymbol{0}\), the equation to solve can be written as:

$$ \boldsymbol{\mu} = \frac{\sum_{i=1}^n \|\boldsymbol{x}_i - \boldsymbol{\mu}\|_2^{\alpha-2}\boldsymbol{x}_i}{\sum_{i=1}^n \|\boldsymbol{x}_i - \boldsymbol{\mu}\|_2^{\alpha-2}} \quad (11) $$

Substituting \(\boldsymbol{\mu}_ {t}\) into the right-hand side and denoting the result by \(\boldsymbol{\mu}_ {t+1}\) gives a fixed-point iteration. Setting \(\alpha=2\) recovers the mean vector, while setting \(\alpha=1\) gives the Weiszfeld iteration. Strictly speaking, one should also prove convergence and uniqueness, but those details are somewhat involved, so I will omit them here.

In addition, there is a less commonly used form of the median vector in high-dimensional spaces, obtained by replacing the Euclidean norm with the L1 norm:

$$ \operatorname*{argmin}_{\boldsymbol{\mu}} \sum_{i=1}^n \|\boldsymbol{x}_i - \boldsymbol{\mu}\|_1 \quad (12) $$

This is called the coordinate-wise median, because it is essentially just taking the one-dimensional median coordinate by coordinate. It is relatively easy to compute, but since it has no particularly clear geometric interpretation, it is used less often.

Summary

This article briefly summarized the concept of the median, its properties, and its extension to higher-dimensional spaces.