中位数（Median）简介
本文转载自苏剑林《中位数（Median）简介》。文中批注为我补充撰写。
最近重新学习了一下中位数的概念，趁新鲜记录一下要点。
做异常值剔除或者裁剪时，我们经常需要一个“基准”，比如对于一堆非负数据，我们可能认为大于基准的50倍就是异常值。那个基准如何选取呢？一个常用的指标是平均值，然而平均值容易被异常值“带偏”，因此以它为基准可能会偏向异常值，从而漏掉一些结果，这时我们可以考虑选取中位数为基准。
基本性质
对于一维数据点 \(x_1, x_2, \cdots, x_n\)，它们的平均值（Mean）定义为：
$$ \text{mean}(x_1, x_2, \cdots, x_n) = \frac{1}{n} \sum_{i=1}^n x_i \quad (1) $$
由于全体数据都直接参与平均计算，所以一旦有几个点特别大，那么平均值也会随之变大，从而干扰异常值的判断。
中位数（Median）的想法是找出数据中的“中立者”作为基准，具体来说，它希望找到一个分隔点，使得大于等于和小于等于该分隔点的数据量一致，所以它也叫做“50%分位数（Quantile-50%）”。如果 \(n\) 是奇数，那么它等于这批数中第 \((n+1)/2\) 大的数，如果 \(n\) 是偶数，那么第 \(n/2\) 大和第 \(n/2+1\) 大之间的任意数字，都可以视为中位数，一般情况下我们取它们俩的平均。
从中位数的定义就可以看出它的抗干扰能力：只要不改变每个数跟中位数的大小关系，那么中位数就不会改变。进一步地，我们可以引入“崩溃点（Breakdown Point）”来描述这种能力，它是“使得结果趋于无穷的最少异常值比例”。对于平均值来说答案是约等于0，因为只需1个点趋于无穷它就趋于无穷；但对于中位数来说是50%，因为要大于一半的数趋于无穷，中位数才趋于无穷。
中位数的缺点是计算起来比均值更为复杂，因为它至少需要对数据做某种程度上的排序，这尤其是对于分布式计算场景并不友好。庆幸的是，大多数场景也不需要非常精确的中位数，所以也有一些挪腾空间，比如局部算中位数后，再全局算一次平均值或者中位数，这样可以明显减少通信。
求解中位数的确比单纯的累加求均值常数项更大，但从计算机算法复杂度的严格意义上讲，寻找中位数绝对不需要进行时间复杂度为 \(O(n \log n)\) 的全局排序。我们可以通过快速选择算法（Quickselect），在期望时间复杂度 \(O(n)\) 内精准找到中位数。Quickselect也可以推广到分布式场景。
其核心思想与快速排序同源：选取一个“基准”元素，将数组划分为比它小和大两部分，此时基准已处在绝对正确的排序位置。求中位数本质是找第\(n/2\)大的元素，只需比对基准索引与目标索引，直接舍弃肯定不存在中位数的那一半数组，仅在另一半中单侧递归。它的计算量级是：\(n + n/2 + n/4 + n/8 + \cdots\) 。根据等比数列求和公式，这个极限严格收敛于 \(2n\)。因此，其期望时间复杂度为线性时间 \(O(n)\)
优化视角
有趣的是，我们可以在优化视角下，将平均值和中位数统一起来：
$$ \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) $$
\(\text{mean}\) 的证明比较简单，就留给读者了，
设 \(f(\mu) = \sum_ {i=1}^ {n} (x_ {i} - \mu)^ {2}\)：
$$ f'(\mu) = -2\sum_{i=1}^n (x_i - \mu) = 0 $$
展开并移项得：$$ \sum_{i=1}^n x_i = n\mu \quad \Rightarrow \quad \mu = \frac{1}{n}\sum_{i=1}^n x_i $$
因二阶导恒为正，故此处取得全局最小值。
这里简单介绍一下 \(\text{median}\) 的证明。不失一般性，假设 \(x_i\) 已经排好序，即 \(x_1 \le x_2 \le \cdots \le x_n\)，设 \(f(\mu) = \sum_{i=1}^n |x_i - \mu|\)，直接求导得：
$$ f'(\mu) = \sum_{i=1}^n \text{sign}(\mu - x_i) = \#\{x_i < \mu\} - \#\{x_i > \mu\} \quad (4) $$
记号 # 表示满足条件的数目。为了找最小值点，我们希望导数尽可能接近于0，即希望大于 \(\mu\) 和小于 \(\mu\) 的 \(x_i\) 数目尽可能相等，这已经符合中位数的思想。更具体的讨论要分情况：
$$ 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) $$
即当 \(n\) 是奇数 \(2k+1\) 时，\(\mu < x_ {k+1}\) 和 \(\mu > x_ {k+1}\) 都会让 \(f(\mu)\) 增大，所以 \(\mu^ {\ast} = x_ {k+1}\)；当 \(n\) 是偶数 \(2k\) 时，\(\mu < x_ {k}\) 和 \(\mu > x_ {k+1}\) 都会让 \(f(\mu)\) 增大，所以 \(\mu^ {\ast} \in [x_ {k}, x_ {k+1}]\)，可以检验当 \(\mu^ {\ast}\) 位于这个区间时，\(f(\mu^ {\ast})\) 始终不变，所以这个区间的数都是 \(\mu^ {\ast}\)。综上，\(\mu^ {\ast}\) 跟中位数的定义完全吻合。
高维空间
从优化视角出发，我们也能理解为什么中位数比均值更能抵御异常值的干扰：假如有一个特别大的 \(x_i\)，那么均值带来的损失是 \((x_i - \mu)^2\)，而中位数则是 \(|x_i - \mu|\)，一般情况下 \((x_i - \mu)^2 \gg |x_i - \mu|\)，所以均值会更偏向于异常值，这样才能尽可能降低损失。
此外，优化视角还有一个好处，就是它容易推广到高维空间。我们知道，平均值的概念很容易推广到高维，对于一批向量 \(\boldsymbol{x}_ {1}, \boldsymbol{x}_ {2}, \cdots, \boldsymbol{x}_ {n}\)，平均向量就是 \(\frac{1}{n} \sum_ {i=1}^ {n} \boldsymbol{x}_ {i}\)，然而中位数的概念却不易直接推广，因为它依赖于排序，但对于向量数据来说很难定义出一个良好的序。
然而，在优化视角下，这种推广是很自然的：
$$ \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) $$
其中 \(\lVert \cdot \rVert_ {2}\) 是欧氏距离，容易证明这样定义出来的平均向量正好是 \(\frac{1}{n} \sum_ {i=1}^ {n} \boldsymbol{x}_ {i}\)，与经验定义相符。至于中位向量，我们也称“几何中位数（Geometric Median）”，从目标函数可以看出，如果 \(n=3\)，那么其实就是经典的“费马点”，因此很多时候我们也直接称中位向量为费马点。
很遗憾，中位向量没有解析解，我们通常是基于如下Weiszfeld迭代来计算：
$$ \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) $$
继续推广
很明显，我们还可以考虑更一般的推广：
$$ \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) $$
在统计学习中，“最小化某种距离（误差）”与“假设数据服从某种概率分布”、“训练仅有一个常数项的机器学习模型来预测数据”是完全等价的。这一桥梁就是最大似然估计（Maximum Likelihood Estimation, MLE）。
假设我们有一组观测数据 \(X = {x_1, x_2, \cdots, x_n}\)，我们要寻找一个中心点 \(\mu\)。我们可以建立如下观测模型：
$$x_i = \mu + \epsilon_i$$
其中，\(\epsilon_i\) 是第 \(i\) 个观测点的随机噪声（误差）。
假设所有噪声 \(\epsilon_i\) 是独立同分布（i.i.d.）的，且其概率密度函数（PDF）为 \(f(\epsilon)\)。那么，观测到单个数据点 \(x_i\) 的概率密度可以写为 \(f(x_i - \mu)\)。对于整个数据集 \(X\)，其联合概率（即似然函数 Likelihood）为所有数据点概率的乘积：
$$L(\mu | X) = \prod_{i=1}^n f(x_i - \mu)$$
我们需要寻找一个 \(\mu\)，使得观测到当前数据集 \(X\) 的概率（似然）最大化。
为了方便计算，我们通常对似然函数取自然对数，并加上负号，将问题转化为最小化负对数似然（Negative Log-Likelihood, NLL）
$$\text{NLL}(\mu) = -\ln L(\mu | X) = -\sum_{i=1}^n \ln f(x_i - \mu)$$
我们回顾文章中给出的广义距离目标函数（以一维为例）：
$$J(\mu) = \sum_{i=1}^n |x_i - \mu|^\alpha$$
要让 MLE 与这个目标函数等价，就必须满足：
$$-\sum_{i=1}^n \ln f(x_i - \mu) \propto \sum_{i=1}^n |x_i - \mu|^\alpha$$
这意味着，对于单个误差 \(\epsilon = x_i - \mu\)，必须满足：
$$f(\epsilon) \propto \exp(-c |\epsilon|^\alpha)$$
这种形式的分布在统计学中已有现成结论，称为广义正态分布，其归一化后的完整形式为：
$$f(\epsilon; \alpha, \beta) = \frac{\alpha}{2\beta \Gamma(1/\alpha)} \exp\left( -\left( \frac{|\epsilon|}{\beta} \right)^\alpha \right)$$
如果我们将这个 \(f(\epsilon)\) 代回负对数似然公式，就会得到：
$$\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)$$
因为 \(\alpha\) 和 \(\beta\) 在求解 \(\mu\) 时是常数，求导时会被消去，所以最小化上述等式，完全等价于：
$$\operatorname*{argmin}_{\mu} \sum_{i=1}^n |x_i - \mu|^\alpha$$
在优化问题中调整的指数 \(\alpha\)，就是在对数据的真实世界物理模型做假设：假设高斯分布就设定 \(\alpha = 2\)；假设拉普拉斯分布就设定 \(\alpha = 1\)。
记 \(f(\boldsymbol{\mu}) = \sum_{i=1}^n \lVert \boldsymbol{x}_ {i} - \boldsymbol{\mu} \rVert_ {2}^ {\alpha}\)，那么：
$$ \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) $$
令 \(\nabla_{\boldsymbol{\mu}} f(\boldsymbol{\mu}) = \boldsymbol{0}\)，那么需要求解的方程可以写成：
$$ \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) $$
将 \(\boldsymbol{\mu}_ {t}\) 代入右端记为 \(\boldsymbol{\mu}_ {t+1}\)，即得不动点迭代法。代入 \(\alpha=2\) 即为平均向量，代入 \(\alpha=1\) 即得Weiszfeld迭代。当然严格来讲还要证明收敛性和唯一性，个中细节比较复杂，这里就省略了。
此外，高维空间的中位向量还有一种比较少用的形式，就是将欧氏范数换成L1范数：
$$ \operatorname*{argmin}_{\boldsymbol{\mu}} \sum_{i=1}^n \|\boldsymbol{x}_i - \boldsymbol{\mu}\|_1 \quad (12) $$
这叫做“坐标中位数（Coordinate-wise Median）”，因为它本质上就是逐分量取一维的中位数，计算起来比较简单，但由于没有明显的几何意义，使用场景相对少一些。
文章小结
本文简单总结了中位数的概念、性质以及它在高维空间中的推广。
中位数（Median）简介