14.3 聚类分析

聚类分析（cluster analysis），也称为数据分割（data segmentation），有多种不同的目的。所有的目的都可以理解为是将一个集合中的对象汇聚或分割到不同的子集或“簇”（cluster）中，从而使每个簇中的那些对象彼此更接近、不同簇中的对象的差异更大。可以用一组度量（特征变量）来描述一个对象，或者也可以用这对象与其他对象之间的关系来描述。此外，有时我们的目的是让簇可以形成一个自然的层级结构。这就需要依次地再对簇本身进行分组，从而使得在每一个层级上同一组内的簇彼此更相似，而不同组之间簇的差异更大。

聚类分析也可用于构建一些描述性的统计量，来明确数据是否是由一些性质存在显著差别的不同子集组成的。为实现这个目标，我们就需要一个可以描述不同簇内的对象差异程度的评估指标。

不论聚类分析的目的是什么，其核心的概念是被聚类的对象个体之间的相似（或不相似）程度的度量。某个聚类方法对对象的分组是依赖于它所用到的相似性定义的。这个定义只能取决于具体问题。这在某种程度上与预测问题（监督学习）中损失函数或成本函数的函数选取有相似之处。一个不准确的预测所带来的成本定义并不是从数据样本中得出的（而是依赖具体的问题和其他领域的知识）。

图 14.4 展示了通过常用的 K 均值算法分成了三个组的模拟数据。在这个例子中，有两个簇并没有很好地被分开；因此“分割”比“聚类”是对这个算法中这部分过程的一个更准确的描述。K 均值聚类从对三个簇中心的猜测值开始，然后交替进行以下两个步骤，直到收敛：

• 为每个数据点，确定与之（按欧式距离）最接近的簇中心，并将其纳入到这个簇。
• 用簇中所有的点的坐标平均值来更新每一个簇的中心。

我们会在后面介绍更多关于 K 均值聚类的细节，包括如何选择簇的数量的问题（这个例子中是三个）。K 均值聚类是一个自上而下的过程，而我们将介绍的其他方法是自下而上的。所有聚类方法的基础是对两个对象之间的距离或差异性测量的选择。在介绍不同的聚类算法之前，我们先探讨下距离的测度。

14.3.1 相似矩阵

有时候数据是直接以对象彼此之间的相似性（proximity / alikeness / affinity）的形式表达的。这些数据可以是 相似性，也可以是 差异性。例如在社会科学的实验中，会请被试评估一个对象与另一个对象不同的确定程度。然后可以对这些评估的样本集做平均，计算出“差异性”。这种类型的数据可以表达为一个 $N\times N$ 的矩阵 $\mathbf{D}$，其中的 $N$ 是对象的个数，矩阵中每个元素 $d_{i{i}^{\prime }}$ 代表的是第 $i$ 和第 $i^{\prime }$ 个对象之间的相似性。这个矩阵是作为聚类算法的输入参数。

大部分算法都假定输入的是一个 差异矩阵（dissimilarity matrix），其所有元素是非负的，对角元素为零 $d_{ii}=0,i=1,2,\dots ,N$。若原始数据中记录了的是相似性，则可用一个恰当的单调递减函数将相似性转换为差异性。此外，大部分算法假定差异矩阵是对称的，所以如果原始的矩阵 $\mathbf{D}$ 不对称，可以用 $(\mathbf{D}+{\mathbf{D}}^T)/2$ 代替。来自主观判断的差异性数据几乎不会符合严格意义上的“距离”定义，因为它不一定会满足距离的三角不等式关系 $d_{i{i}^{\prime }}\le d_{ik}+d_{i^{\prime }k},\mathrm{\forall }k\in \{1,\dots ,N\}$。因此对于这类数据，一些依赖距离假设的算法就无法使用了。

14.3.2 基于属性的差异度量

通常我们的数据是对变量（也称为 属性，attribute） $j=1,2,\dots ,p$ 的测量值 $x_{ij},i=1,2,\dots ,N$。既然大部分常用的聚类算法都以差异矩阵作为输入数据，我们需要首先构建观测样本两两之间的差异性。在大多常见的场景中，我们先定义第 $j$ 个属性取值之间的差异函数为 $d_j(x_{ij},x_{i^{\prime }j})$，然后再定义对象 $i$ 和 $i^{\prime }$ 之间的差异性为：

目前最常用的选择是平方距离：

不过也有其他可能的选择，不同的选择会产生可能不同的结果。而对非数量的属性（例如分类数据），平方距离可能不再适用。此外，有时会希望对属性给予不同的权重，而不是像在 14.20 中那样都是相同的权重。

这里先根据属性类型的不同介绍几个其他的距离定义：

• 数值变量（quantitative variable）
  这类变量或属性的测量值是在连续的实数域上。自然地会将它们之间的“误差”定义为它们之间绝对差值的一个单调递增函数： $$d(x_i,x_{i^{\prime }})=l(|x_i-x_{i^{\prime }}|)$$ 除了平方误差损失 $(x_i-x_{i^{\prime }}{)}^2$ 之外，另一个常见的选择是单位函数（即绝对误差）。平方误差对大的差异比对小的差异的重视程度更高。另外在聚类算法中也可使用相关性： $$\begin{array}{}\text{(14.22)}& \rho (x_i,x_{i^{\prime }})=\frac{\sum _j(x_{ij}-{\overline{x}}_i)(x_{i^{\prime }j}-{\overline{x}}_{i^{\prime }})}{\sqrt{\sum _j(x_{ij}-{\overline{x}}_i{)}^2\sum _j(x_{i^{\prime }j}-{\overline{x}}_{i^{\prime }}{)}^2}}\end{array}$$ 其中 ${\overline{x}}_i=\sum _j{x}_{ij}/p$。需要注意这里是对变量（维度）的平均，而不是对观测样本的平均。若已将观测样本进行了标准化处理，那么： $$\sum _j(x_{ij}-x_{i^{\prime }j}{)}^2\propto 2(1-\rho (x_i,x_{i^{\prime }}))$$ 因此，基于相关性（相似性度量）的聚类也就等价于基于平方距离（差异性度量）的聚类。

• 有序变量（ordinal variable）
  这个类型的变量通常是连续的整数，而且它的取值是在一个有序集合中。例如，学习成绩（A、B、C、D、F），偏好程度（厌恶、反感、一般、喜欢、极好）。排序数据也是一个特殊的有序变量。有序变量的误差测度定义，一般会按照其原始取值的排序，用下式代替其 $M$ 个原始取值： $$\begin{array}{}\text{(14.23)}& \frac{i-1/2}{M},i=1,\dots ,M\end{array}$$ 然后在这个取值尺度上将它们按照数值变量的方式处理。

• 分类变量（catgorical variable）
  对于无序的分类变量（unordered categorical，nominal），取值两两之间的差异度必须显式地直接定义出来。若变量有 $M$ 个不同的可能取值，这些差异度可以表达为一个对称的 $M\times M$ 矩阵，满足 $L_{r{r}^{\prime }}=L_{r^{\prime }r}$、$L_{rr}=0$、$L_{r{r}^{\prime }}\ge 0$。最常用的选择是 $L_{r{r}^{\prime }}=1,\mathrm{\forall }r\ne r^{\prime }$，不过也可以用不相等的损失值来区别对待不同的误差。

14.3.3 对象的差异性

接下来我们将 $p$ 个单独的属性差异性 $d_j(x_{ij},x_{i^{\prime }j})$ $j=1,2,\dots ,p$ 组合起来，构建包含这些属性的两个对象或观测样本 $(x_i$，$x_{i^{\prime }})$ 之间一个总体差异性的测量 $D(x_i,D_{i^{\prime }})$。这一般都是通过一个加权平均（凸性组合）的方法。

其中 $w_j$ 为分配给第 $j$ 个属性的权重，控制着这个变量在对象之间的总体差异性中的相对影响作用。权重的选择应该基于具体问题的分析。

必须要指出的是将每个变量的权重 $w_j$ 设置为相同的数值（例如，$w_j=1,\mathrm{\forall }j$ 并不一定意味着所有的属性（对差异性）有一样的影响作用。第 $j$ 个属性 $X_j$ 对 14.24 的对象差异性 $D(x_i,x_{i^{\prime }})$ 的影响，依赖于它对数据集所有样本两两之间对象差异性测度平均值的相对贡献：

其中第 $j$ 个属性的平均相似性为：

因此，第 $j$ 个变量的相对影响作用是 $w_j\cdot {\overline{d}}_j$；若令 $w_j\sim 1/{\overline{d}}_j$，所有的属性在对象之间的总体差异性中有相同的影响作用。例如，若有 $p$ 个数值变量以及在每个坐标上使用平方误差距离，则式 12.24 就变成了 ${\mathbb{R}}^p$ 空间上两个点之间的（加权）平方欧式距离：

其中的数值变量就作为了坐标轴。这里的式 14.25 就变成了：

其中的 ${\mathrm{var}}_j$ 为 $\mathrm{Var}(X_j)$ 的样本估计。因此，每个变量的相对重要程度与它在数据集上的方差呈正比。一般来说，无论变量的类型而令所有属性的 $w_j=1/{\overline{d}}_j$，会使每个变量在对象之间 $(x_i,x_{i^{\prime }})$ 的总体差异性的影响相同。尽管这看似合理而且有时是被推荐的做法，但是这可能非常地事与愿违。如果目的是将数据分割成相似对象组成的组中，那么所有的属性在对象之间差异性（依赖具体问题的）概念上的贡献程度可能有差异。在不同的问题背景中，某些属性的取值差异可能对实际的对象差异性影响更大。

如果目的是发掘出数据中的天然组群，某些属性可能比其他属性携带了更多的分组趋势。应该给对群组分割更有价值的变量在定义对象相似性中分配更高的影响因子。而在这里给所有的属性分配相同的影响因子会将群组变得模糊不清，以至于聚类算法无法发现这些群组。[图 14.5] 展示了一个例子。

虽然对个体属性差异性 $d_j(x_{ij},x_{i^{\prime }j})$ 和它们的权重 $w_j$ 一般性的描述比较容易理解，但在每个问题中对具体场景的深思熟虑是不可或缺的。在一个成功的聚类中，选择一个合理的差异性度量要比选择聚类算法重要得多。由于这取决于具体问题领域中的知识和常识，不便于进行一般性地研究，所以与算法本身相比，聚类相关的文献中较少强调这方面的问题。

最后一点，观测样本中通常在一个或多个属性上有缺失值。在式 14.24 的差异性计算中最常用的缺失值处理方法，是在计算样本 $x_i$ 和 $x_{i^{\prime }}$ 之间的差异性时忽略每个至少缺失了一个值的那对（属性）观测值 $(x_{ij}，x_{i^{\prime }j}$。当两个观测样本没有同时不缺失的属性值时，这个方法不再可用。这时可以将这两个样本都从聚类分析中排除。另一种方法是用每个属性非缺失数据的均值或中位数来插补（impute）缺失值。而对于分类变量而言，如果认为在相同的变量上有缺失值的两个对象也是一种相似性的话，可以将“缺失”作为一种分类取值。

14.3.4 聚类算法

聚类分析的目的是将观测样本划分到若干个组（group）或簇（cluster）中，使得同一个簇中的两两差异性一般要小于不同簇中的差异性。聚类算法可分成三个不同的类型：组合算法、混合模型、和众数搜索。

组合算法（combinatorial algorithm） 不需要直接参照一个底层的概率模型，而直接对观测样本进行操作。混合模型（mixture modeling） 则假设数据来自一个概率密度函数所描述的某个总体样本的独立同分布（i.i.d.）抽样。这个密度函数是由许多成分密度函数混合而成的一个参数化的模型；每一个成分密度函数代表了一个簇。然后使用最大似然或相应的贝叶斯方法在数据上拟合这个模型。众数搜索（mode seeking），或“凸块搜索”（bump hunting），从非参数方法的角度试图直接估计出概率密度函数不同的众数。按最近的众数点划分的观测点就组成了相应的众数点处的簇。

第 6.8 节介绍了混合模型。第 9.3 节和第 14.2.5 节中的 PRIM 算法，是众数搜索或凸块搜索中的一种。接下来将介绍组合算法。

14.3.5 组合算法

最常见的聚类算法是直接将每个观测样本分配到一个组或簇中，而不需要考虑这个数据背后的概率模型。将每个观测点用整数 $i\in \{1,\dots ,N\}$ 做唯一标识。预先假定存在 $k<N$ 个簇。则要将每个观测点分配到一个且仅限一个簇中。这些分配规则可以表达为一个多对一映射；或写为 $k=C(i)$，将第 i 个观测点分配到第 $k$ 个簇的 编码器（encoder）。聚类算法就是在基于观测点两两之间的差异性定义 $d(x_i,x_{i^{\prime }})$，寻找可实现分组目标（后面将详细说明）的特定编码器 $C^{\ast }(i)$。如上文所述，差异性的定义和分组的目标是由使用者指定的。一般来说，编码器是通过它为每个观测点 $i$ 分配的（簇标识）值来描述的。因此，这个编码的“参数”就是为 $N$ 个观测点的每一个分配的簇标识值。然后通过对一个描述了聚类目标未达成的程度的“损失”函数进行最小化，来对这些参数进行调节。

一个方法就是直接地指定一个量化的损失函数，然后通过某个组合最优化算法对其进行最小化。既然聚类的目标是将接近的点分配到同一个簇中，那么损失函数（或能量函数）自然的形式是：

这个准则函数描绘了同一个簇中的观测样本点彼此接近的程度。有时它被称作为“簇内分散度”（within cluster point scatter），因为分散度可分解为：

或者写为：

其中 $d_{i{i}^{\prime }}=d(x_i,d_x{i}^{\prime })$。这里的 $T$ 为“总分散度”（total cuslter），在给定的数据上它是一个不依赖于聚类结果的常数。而下式则为“簇间分散度”（between cluster point scatter）

当分配到不同簇中的观测样本点之间的距离较远时，这个分散度一般会比较大。因此可见：

对 $W(C)$ 的最小化等价于对 $B(C)$ 的最大化。

组合最优化的聚类分析从原理上很简单。只需要在 $N$ 个数据样本点和 $K$ 个簇的所有可能分配方式中选择使 $W$ 最小或使 $B$ 最大的即可。不过这样全量列举的最优化方式只对较小的数据集是可行的。不同的分配方式的个数为（Jain and Dubes, 1988）：

例如，$S(10,4)=34105$，尚且可行。但 $S(N,K)$ 会随着它的参数取值变大而迅速变大。比如 $S(19,4)\simeq {10}^{10}$，而且几乎所有分类问题的数据集都要大于 $N=19$。因此，实用的聚类算法可以只对所有可能分类器 $k=C(i)$ 中非常小的一部分进行考察。它们的目标是能识别出可能包含了最优解的一个小的集合，或至少是包含了次优的分配规则。

这些可行的最优策略是基于重复的贪心下降（greedy descent）。首先指定一个初始化的分配规则。然后在每一次迭代步骤中，对簇的分配规则进行改动，使得准则函数的取值较之前有所提升。这个类型的聚类算法在每次迭代中修改簇分配规则的方法上有所区别。当在某次迭代中对分配规则的改动无法得到提升时，算法终止并将当前的分配规则作为它最终的解。由于在每一次迭代中的观测点和簇的分配规则是基于上一次迭代结果的微小变动，所以在所有可能的分配规则（式 14.30）中只有很小的一部分会被考虑到。然而，这些算法收敛到的局部最优可能与全局最优的差距很大。

14.3.6 K 均值

K 均值（K-means） 算法是最常用的迭代下降聚类方法之一。它适用于所有变量都是数值变量并且使用欧式距离的平方作为差异性度量的场景中。

请注意也可以通过重新定义 $x_{ij}$ 的值的形式实现加权欧式距离（练习 14.1）。

则式 14.28 的簇内分散度可写为：

其中的 ${\overline{x}}_k=({\overline{x}}_{1k},\dots ,{\overline{x}}_{pk})$ 为第 $k$ 个簇内的均值向量，并且 $N_k=\sum _{i=1}^{N}I(C(i)=k)$。因此，可以使准则的最小化的 $N$ 个观测点与 $K$ 个簇的分配方式，也会使簇内观测点与簇内均值（簇内点的平均）平均距离最小化。

这个最小化问题为：

上式可通过一个迭代下降算法求解。对任意的观测点集合 $S$ 来说：

因此我们可以通过对一个扩大的最优化问题求解而得到 $C^{\ast }$：

这可以用算法 14.1 中给出的交替最优化过程进行最小化。

算法 14.1 K 均值聚类

1. 给定一个簇分配规则 $C$，对总的簇内方差 14.44 进行最小化，则 $\{m_1,\dots ,m_K\}$ 的解为当前簇内的均值（见式 14.32）。
2. 给定一组当前的均值 $\{m_1,\dots ,m_K\}$，式 14.33 的最小化解即为将每个观测样本点分配给距离（当前）簇均值最近的簇： $$\begin{array}{}\text{(14.34)}& C(i)=\underset{1\le k\le K}{\arg min}\Vert x_i-m_k{\Vert }^2\end{array}$$
3. 迭代进行步骤 1 和步骤 2，直到分配规则不再改变。

步骤 1 和步骤 2 都会降低准则 14.33 的值，所以可保证过程的收敛。不过这个结果有可能是一个次优的局部最小化解。Hartigan and Wong (1979) 的算法更进一步地保证了最终结果中不存在可使目标函数继续降低的变更单个观测点所属簇的操作。此外，使用者可以随机生成多个初始化均值来启动算法，然后在其中选择目标函数值最小的一个解。

图 14.6 展示了在图 14.4 的数据上 K 均值的一些迭代结果。其中的圆圈代表着簇中心点。直线划分出了观测点的区域，每个区域代表着距离每个中心点最近的观测点。这个分割图被称作 沃罗诺伊图（Voronoi tessellation）¹。在 20 次迭代后，计算过程达成了收敛。

14.3.7 高斯混合与 K 均值软聚类

K 均值聚类计算过程与估计某个特定高斯混合模型使用的 EM 算法有紧密的联系。（参考第 6.8 节和第 8.5.1 节）。EM 算法中的期望（E）步骤基于在每个混合成分分布下的相对密度为每个数据点分配“责任值”（responsibility），而后在最大化（M）步骤基于这些责任值重新计算成分密度函数的参数。假设我们指定了有 $K$ 个混合成分，每个成分高斯密度函数的协方差矩阵为常数矩阵 ${\sigma }^2\mathbf{I}$。那么在每个混合成分分布下的相对密度，为数据点与混合中心点欧式距离的一个单调函数。因此在这个背景下 EM 算法即为一个“软性”版本的 K 均值聚类，为数据点属于某个簇中心分配一个概率值（而不是确定的规则）。随着方差 ${\sigma }^2\to 0$，这些概率会变成 0 或 1，这两个方法是一样的。具体细节见练习 14.2。图 14.7 在实数线上的两个簇对此进行了演示。

14.3.8 示例：人类肿瘤微阵列数据

我们在第一章中介绍过的人类肿瘤微阵列（microarray）数据上应用 K 均值聚类。这是一个高维空间聚类问题。

数据为一个 $6830\times 64$ 的实数值矩阵，每个元素代表了一个基因（一行）在一个样本（一列）上的表达水平测量。我们要对样本进行聚类，每个样本是一个长度为 6830 的向量，对应着 6830 个表达水平值。每个样本都有一个标签，例如 breast（乳腺癌）、melanoma（黑色素癌）、等等；在聚类中我们不会使用这些标签，而是在过后检查各个标签落入到哪个类别中。

我们对 $K$ 的取值从 1 到 10 分别应用了 K 均值聚类，然后计算每次聚类的总簇内（离差）平方和，如图 14.8 所示。一般来说会在平方和曲线（或它的对数曲线）中寻找一个折点来作为聚类的最优簇数量（见第 14.3.11 节）。这里没有明确的显示，为了演示我们选择了 $K=3$，三个聚类的结果列在表 14.2 中。

表 14.2：人类肿瘤数据：K 均值聚类的三个簇中的每个类型癌症样本的个数。

可见聚类方法成功地将同一类癌症的样本归在同一组中。实际上第二个簇中的两个乳腺癌样本在之后被发现是误诊，是发生了转移的黑色素癌。不过在这个应用中 K 均值聚类也有不足之处。首先，它没有给出簇内对象的一个线性排序，我们只是简单地按字母顺序列举出来。其次，当改变簇的数量 $K$ 时，簇中的对象可能随意地发生变化。比如说更改为四个簇，则新的簇不需要嵌套在上面的三个簇中。由于这些问题，（稍后将介绍的）层次聚类可能更适用于这种应用场景。

14.3.9 向量量化

K 均值聚类算法在看似并不相关的图片和信号压缩领域也是一个关键的工具，特别是在向量量化（vector quantization，VQ）问题中（Gersho and Gray, 1992）。图 14.9² 的左图为著名统计学家罗纳德·费希尔爵士（Sir Ronald Fisher）的一张数字化照片。其中包含了 $1024\times 1024$ 个像素点，每个像素点为从 0 到 255 的灰度值，因此每个像素点需要 8 位（bit）的储存空间。整个图片占据了 1M（megabyte）的储存空间。中图为左图的一个 VQ 压缩的版本，只需要原储存空间的 0.239（存在一些图片质量损失）。右图为更进一步的压缩，只需要原储存空间的 0.0625（存在相当大的图片质量损失）。

这里的 VQ 版本首先将图片分解成很多小方块，在这个例子中是 $2\times 2$ 的的像素块。这些 $512\times 512$ 个包含了 4 个数字的块都可作为 ${\mathbb{R}}^4$ 空间上的向量。在这个空间应用 K 均值聚类算法（在这个场景中也被称为劳埃德（Lloyd）算法）。中图中使用了 $K=200$，而右图中使用了 $K=4$。这 $512\times 512$ 中的每个块（或点）可以用它距离最近的簇中心来近似，这被称作码字（codeword）。聚类的过程被称作编码（encoding），中心点的集合被成为码书（codebook）。

为了描述近似后的图片，我们需要为每个块记录它的近似值在码书中位置。每个块需要 ${\log }_2(K)$ 位空间。同时我们也需要记录码书本身，$K\times 4$ 个实数（一般可忽略不计）。总的来说，压缩后的图片占原图片空间的比例为 ${\log }_2(K)/(4\times 8)$（当 $K=200$ 时为 0.239，当 $k=4$ 时为 0.063）。这一般表达为以比特每像素（bits per pixel）衡量的码率（rate）：${\log }_2(K)/4$，这里分别为 1.91 和 0.50。从中心点（码书）构建近似后图片的过程被称作解码（decoding）。

为什么 VQ 图片压缩可以奏效？其原因是对于一般的日常图片，比如照片，很多块是一样的。在这个例子中，有很多几乎为纯白色的块，而且类似地也有很多不同灰度的纯灰色块。这些信息只需要每种保留一个块，以及指向这种块的像素指针。

由于被近似的图片是原图的质量降级版本，上面介绍的也被称作有损（lossy）压缩。通常可以用均方误差（mean squared error）测量这个降级（degradation）或失真（distortion）。在这个例子中，$K=200$ 对应着 $D=0.89$，$K=4$ 的 $D=16.95$。更一般地，可以用一个率压缩曲线（rate/distortion curve）来评估这个权衡。也可以利用块的聚类进行无损（lossless）压缩，也仍是在利用重复的像素块。如果对这个例子中的原始图片进行无损压缩，那么可达到的最大码率为 4.48 比特每像素。

上文中说在码书中识别 $K$ 个码字中的每一个都需要 ${\log }_2(K)$ 个比特。这里使用的是定长编码，如果图片中的某些码字比其他码字出现的次数多得多，那么这个方式是低效的。根据香农（Shannon）编码理论可知，一般来说变长编码更佳，那么码率就变成了 $-\sum _{\ell =1}^K{p}_{\ell }{\log }_2(p_{\ell })/4$。分子项为图片中码字分布 $p_{\ell }$ 的熵值（entropy）。使用了变长编码后，码率分别下降到 1.42 和 0.39。已有很多对 VQ 的推广被发展出来了：例如在第 14.3.12 节会稍微提及的树结构 VQ，通过一个自上而下、二均值类型的算法寻找中心点。这个推广可以进行逐次细化的压缩。更多细节可参考 Gersho and Gray (1992)。

14.3.10 K 中心点

如上所述，K 均值算法适用于以欧式距离平方 $D(x_i,x_{i^{\prime }})$（式 14.21）作为差异性度量的情况。这就要求所有的变量都是数值类型。此外，使用欧式距离平方令最大的距离有最高的影响力。这会导致计算过程在出现非常大距离的异常点时缺少稳健性。这些限制可通过额外的计算而去除。

K 均值算法中仅在以下用到了欧式距离平方假设：最小化步骤的式 14.32；式 14.33 中的簇中心向量 $\{m_1,\dots ,m_K\}$ 为当前簇中分配对象的均值。通过将这个步骤的式 14.33 中替换为一个对 $\{m_1,\dots ,m_K\}$ 的显式的最优化，可以将算法推广到任意的差异性定义 $D(x_i,x_{i^{\prime }})$ 上。在大多情况下，每个簇的中心点被限制为只能取自簇中的一个观测点，算法 14.2 汇总了计算过程。这个算法假设了有属性（变量）数据，不过这个方法也可以使用在只用相似矩阵描述的数据上（第 14.3.1 节）。这时就不需要再显式地计算出簇中心，而只需要记录簇中心的标识 $i_k^{\ast }$ 即可。

算法 14.2 K 中心点聚类

1. 给定一个簇分配规则 $C$，找出每个簇中距离同簇中其他点的总距离最小的观测点： $$\begin{array}{}\text{(14.35)}& i_k^{\ast }=\underset{\{i:C(i)=k\}}{\arg min}\sum _{C(i^{\prime })=k}D(x_i,x_{i^{\prime }})\end{array}$$ 则 $m_k=x_{i_k^{\ast }},k=1,2,\dots ,K$ 为当前步骤的簇中心估计。
2. 给定一组当前的簇中心 $\{m_1,\dots ,m_K\}$，将每个观测点分配至（当前步骤）距离它最近的簇中心，从而使总误差最小化： $$\begin{array}{}\text{(14.36)}& C(i)=\underset{1\le k\le K}{\arg min}D(x_i,m_k)\end{array}$$
3. 迭代进行步骤（1）和步骤（2），直到分配规则不再改变。

在每个临时的簇 $k$ 中求解式 14.32 需要的计算量与这个簇中分配的观测点数量成正比，然而求解式 14.35 的计算量则增加到了 $O(N_k^2)$。给定一组簇“中心”，$\{i_1,\dots ,i_K\}$，输出新分配规则：

所需要的计算量和之前一样与 $K\cdot N$ 成正比。因此，K 中心点要比 K 均值需要的计算量大的多。

用式 14.37 替换式 14.35 后（算法 14.2）就成为了求解下面这个最小化问题的一个特定的启发式搜索策略（heuristic search strategy）。

Kaufman and Rousseeuw (1990) 提出了另一个直接求解式 14.38 的策略：暂时地用一个非中心点的观测点替换每一个中心 $i_k$，然后选择那个使式 14.38 中的准则函数值降低程度最大的替换操作。重复这样的操作，直到无法找到使准则优化的替换操作。Massart et al. (1983) 推导出了一个可找出式 14.38 的全局最小化解的分支定界（branch and bound）的组合方法，但最只在比较小的数据集是可行的。

例：国家差异性

表 14.3：一个政治学的调研数据：一些政治学的学生的调查问卷中国家之间两两差异性的得分平均值。

这个示例来自 Kaufman and Rousseeuw (1990)，在一项研究中让一些政治学的学生给出 12 个国家的两两差异度：比利时（BEL）、巴西（BRA）、智利（CHI）、古巴（CUB）、埃及（EGY）、法国（FRA）、印度（IND）、以色列（ISR）、美国（USA）、苏联（USS）、南斯拉夫（YUG）、扎伊尔（ZAI）。表 14.3 列出了差异度评分的平均值。在这个差异矩阵上使用 3 中心点聚类。请注意在这里无法使用 K 均值，因为只有距离矩阵而没有原始的观测点。图 14.10 的左图展示了根据 3 中心点聚类进行分块并且重新排列的差异矩阵。右图为二维的多维尺度分析图，其中的颜色标记了三个中心点的聚类结果（关于多维尺度分析请阅读第 14.8 节）。两个图中都可看到三个簇的分隔清晰，但从多维尺度分析图中看“埃及”落在了两个簇中间的位置。

14.3.11 实践问题

应用 K 均值或 K 中心点聚类算法，需要选择簇的个数 $K^{\ast }$ 和初始化条件。后者可通过多种性质来定义：初始中心向量 $\{m_1,\dots ,m_K\}$、初始中心点标识 $\{i_i,\dots ,i_K\}$、或初始簇分配编码器 $C(i)$。通常指定中心向量或中心点的方式更方便。备选的方法有简单的随机选取，也有基于前向分布分配的策略。后者是在每一步中，假设 $i_1,\dots ,i_{k-1}$ 已在上一步中给出，通过对准则 14.33 或 14.38 最小化来选择一个新的中心点 $i_k$。进行 $K$ 次这样的步骤就可得到 $K$ 个初始化中心点，可以用于聚类的最优化算法中。

簇个数 $K$ 的选择依赖于具体的聚类目的。在数据分割问题中 $K$ 通常由问题中一部分背景而决定。例如一个公司可能想要雇用 $K$ 个销售人员，则目标是将一个客户的数据库分割成 $K$ 个部分，每个销售员对应着一个部分，并使每个销售员分配到的客户尽可能相似。不过也有很多时候，聚类分析是为了明确观测到的数据是否自然地落在若干不同组而提供一个描述性的统计论证。这里这些组的个数 $K^{\ast }$ 是未知的，簇的个数和分配规则都需要从数据中估计得来。

数据驱动的 $K^{\ast }$ 估计方法一般是考察簇内差异性 $W_k$ 对簇个数 $K$ 的函数。对 $K\in \{1,2,\dots ,K_{max}\}$ 各自得出聚类的解。对应的差异值 $\{W_1,W_2,\dots ,W_{K_{max}}\}$ 一般会随着 $K$ 的增加而降低。即便是在一个独立的测试集上，这也是成立的。因为更多的簇中心点可以更密集地填充在特征空间上，从而所有的数据点都会接近于簇中心。所以在监督学习的模型选择中非常有用的交叉验证方法，在这里却不可用了。

这个方法底层的直觉思想如下。如果真实的数据模型中存在 $K^{\ast }$ 个（由差异性度量定义的）不同的组，那么当 $K<K^{\ast }$ 时算法得出的每个簇中都包含了一个或多个真实的组。也就是说，聚类算法的解不会将同一个真实分组中的观测点分配到两个不同的簇中。在这种情况下，随着指定的簇个数逐个增加，真实的组被逐个地分配到独立的簇中，聚类解的准则取值会逐个大幅地下降， $W_{K+1}\ll W_K$。而当 $K>K^{\ast }$ 时，肯定有一个（或多个）估计的簇将至少一个真实的组分割成了两部分。这时随着 $K$ 的进一步增加而带来的准则下降幅度会更小。与将两个可分隔真实组的集合恰当地分成两部分相比，将一个内部观测点都很接近的真实组分开所引起的准则取值下降幅度会更下。

那么可以理解在这个情况下，在 $K=K^{\ast }$ 处准则取值的逐个差分 $W_K-W_{K+1}$ 会有一个大幅度的下降。也就是说：

那么就可以通过找到 $W_k$ 对 $K$ 的函数曲线中的一个“转折点”来得出一个 $K^{\ast }$ 的估计值 ${\hat{K}}^{\ast }$。与聚类计算过程的其他方面类似，这个方法也是有一些“启发性”的（heuristic）。

近期 (Tibshirani et al., 2001b) 提出了 间隙（Gap）统计量。在包含了样本数据的一个长方体中生成均匀分布的数据，然后得出一个 $\log W_K$ 曲线。将样本数据的 $\log W_K$ 曲线与生成的曲线做对比。它将两个曲线之间的差距最大位置的 $K$ 作为最优簇个数的估计值。这基本上是一个可以自动化地定位上文提到的“转折点”的方法。它在真实数据只有一个组时依然相当有效，这时会倾向于估计最优簇个数为一。而在这种场景下大多数其他计算方法会失效。

图 14.11 展示了在图 14.4 模拟数据上使用了间隙统计量的结果。左图为 $K=1,2,\dots ,8$ 个簇的 $\log W_K$ 的曲线（绿色），以及 20 次模拟的均匀分布数据的 $\log W_K$ 期望值曲线（蓝色）。右图展示了间隙统计量曲线，即为期望值曲线减去观测样本曲线。图中也标记了误差区间，区间的半个宽度为 $s_K^{\prime }=s_K\sqrt{1+1/20}$，其中的 $s_K$ 为 $\log W_K$ 在 20 次模拟数据中的标准差。在 $K=2$ 个簇时，间隙曲线最大。若将 $K$ 个簇的间隙曲线记为 $G(K)$，则 $K^{\ast }$ 估计值的正式规则为：

这个结果是 $K^{\ast }=2$，在图 14.4 上看也是合理的。

14.3.12 层次聚类

K 均值或 K 中心点聚类算法的结果依赖于簇个数的选择和初始化的分配规则配置。与之相反，层次聚类（hierarchical clustering） 方法不需要这些配置。然而它需要指定一个基于两组观测点之间的两两差异性的两个（不相交）组之间的差异性度量。正如它的名称，这些方法可得出层次结构，其中每一级中的簇都是由下一级中的簇合并而来的。在最低一级，每个簇只包含一个观测样本。在最高一级，只有一个包含了所有数据的簇。

层次聚类的策略可分为两个基本的模式：自下而上的 凝聚（agglomerative），和自上而下的 分裂（divisive）。凝聚方法从最底部开始，递归地在每一级选择一对簇将其合并起来。这样在上一个级别中簇的个数就变少了一个。被选择合并的一对簇是所有簇中簇间差异性最小的两个。分裂方法从最顶部开始，递归地在每一级选择一个簇将其分成两个簇。分裂的选择是可以使得到的两个新簇之间的差异性最大。两种模式的层次结构中都会有 $N-1$ 个等级。

层次结构中的每一级都代表了将观测数据划分到不相交簇中的一个特定的数据分组。层次结构的整体就形成了一个这些数据分组的有序序列。至于哪一个级别能最好地表现一个“自然”聚类，也就是在哪个级别上每个组内部的观测点足够相似而不同组之间的观测点相差足够大，则需要使用者自行决定。之前介绍的间隙统计量也可在这里用到。

递归二元分裂/凝聚可以通过一个有根二叉树表现出来。树中的节点代表了分组。根节点代表了整个数据样本集合。$N$ 个终节点每个都代表了一个单独的观测点（单元素簇）。每个非终节点（父节点）都有两个子节点。在分裂聚类中，两个子节点是从父节点分裂而形成的两个组；在凝聚聚类中，两个子节点是合并成父节点的两个组。

大部分的凝聚方法和一些分裂方法（从下向上观察）具有一个单调的性质。也就是说，被合并的两个簇之间的差异性随着合并级别的增加而单调递增。因此，在绘制二叉树时可以使每个节点的高度与其子节点之间的差异性取值成正比。只包含一个观测点的终节点在图中高度为零。这种类型的图形称之为 树状图（dendrogram）。

树状图以图像的形式提供了一个可解释性很高的对层次聚类的完整描述。这也是层次聚类方法非常流行的主要原因之一。

图 14.12 展示了微阵列数据上使用平均链接的凝聚聚类得到的树状图；本节稍后会更详细介绍关于凝聚聚类和这个示例的内容。若在一个特定的高度水平地将树状图截断，数据就会被分隔成与水平线相交的垂直线所代表的不相交簇。如果当最优组间差异性超过某个阈值点（树状图的高度）后就终止合并的过程，这些就是所得到的簇。

较高一级合并的（组间差异）值会高于该节点下面较低级别合并的值，前者是自然（真实）聚类的备选结果。可见这个过程可以发生在不同的级别上，也就形成了一个聚类的层次结构，即为簇嵌套着簇。

这样的树状图通常被看作是对数据本身的一个图像的概括，而不是对算法结果的一个描述。不过应该小心地看待这种理解方式。首先，不同的层次聚类方法（下文会详述），或者数据中的一些小变更，都可能会得出非常不同的树状图。其次，只有当观测点两两之间的差异性真的存在可被算法捕捉到的层次结构时，这样的概括图才有意义。其实无论数据中是否真的存在这样的层次结构，层次聚类方法都会得出层次结构。

树状图表现出的层次结构是否真的表达了数据本身的结构，可以通过 共表相关系数（cophenetic correlation coefficient） 来做判断。这个系数是输入给算法的 $N(N-1)/2$ 个观测点之间两两的差异性 $d_{i{i}^{\prime }}$与它们从树状图推导出对应的共表差异性 $C_{i{i}^{\prime }}$ 之间的相关系数。两个观测点 $(i,i^{\prime })$ 之间的共表差异性 $C_{i{i}^{\prime }}$ 是这两个观测点 $i$ 和 $i^{\prime }$ 第一次合并到同一个簇中时两个组间的差异性。

共表差异性是一个约束性很强的差异性度量。首先，观测样本上的 $C_{i{i}^{\prime }}$ 必定会有很多重复的值，因为在所有的 $N(N-1)/2$ 个差异性取之中（至多）只有 $N-1$ 个唯一值。其次，这些差异性值需要满足 超度量不等式（ultrametric inequality），即对任意三个观测点 $(i,i^{\prime },k)$：

从几何上做个示例，假设数据为一个欧式坐标系统中的点。为了使数据点之间的距离能满足式 14.40，所有的三个点形成的三角形都必须是等腰三角形，而且三角形的底要小于它的腰（Jain and Dubes, 1988）。因此一般一个数据集的差异性度量，尤其是当其中没有太多重复值时，都不太可能与从树状图中得出相应的共表差异性比较相似。所以树状图应该被更多地看作是对应用某个特定的算法后所得出的层次聚类结构的一个描述。

凝聚聚类

凝聚聚类（agglomerative clustering）算法从每个观测点代表一个单元素簇开始。在 $N-1$ 个步骤的每一步中，将两个（差异性最小的）簇合并为一个簇，使高一级中簇的个数少一个。所以需要先定义两个簇（观测点的集和）之间的差异性度量。

令 $G$ 和 $H$ 代表两个簇（组）。那么 $G$ 和 $H$ 之间的差异性 $d(G,H)$ 是基于那些观测点两两之间的差异性 $d_{i{i}^{\prime }}$ 计算来的，距离中的一个 $i$ 来自 $G$，另一个 $i^{\prime }$ 来自 $H$。单链接（single linkage，SL） 凝聚聚类取其中最近（差异性最小）的一个作为组间差异性：

这也通常被称为 最近邻（nearest-neighbor） 方法。全链接（complete linkage，CL） 凝聚聚类，或 最远邻（furthest-neighbor） 方法，则取其中最远（差异性最大）的一个作为组间差异性：

组平均（group average，GA） 凝聚聚类则使用组中的平均差异性：

其中的 $N_G$ 和 $N_H$ 为各自组中观测点的个数。尽管对凝聚聚类中的组间差异性的定义还有很多其他的方案，以上三个是最常用的。图 14.13 展示了这三种方法。

如果数据的差异性 $\{d_{i{i}^{\prime }}\}$ 有较明显的簇结构，即每个簇都很紧致（compact）并且与其他簇分隔清晰，那么三个方法都会得到相似的结果。紧致的簇，意味着簇内所有的观测点都相对接近（差异性小），而不同簇之间观测点的差异性大。如果事实并非如此，则三个方法的结果可能不同。

单链接的定义 14.41 只需要两个组中的一个差异性 $d_{i{i}^{\prime }},i\in G,i^{\prime }\in H$ 比较小，就可以认定这两个组 $G$ 和 $H$ 是相近的，而不需要考虑两组内其他观测点之间的差异性。因此它会倾向于在相对低的阈值下合并中间地带被一些观测点连接起来的两个组。这种现象被称为 连锁（chaining），通常被看作是这个方法存在的一个缺陷。单链接得到的簇可能会违背“紧致”的性质，即让每个簇内所有的观测点都彼此相似。如果将一组观测值的 直径（diameter）$D_G$ 定义为组中元素之间的最大差异性：

那么单链接可能会得出直径非常大的簇。

全链接的定义 14.42 则是另一个极端。只有当两个组的并集中所有的观测点都比较相似时，两个组 $G$ 和 $H$ 才会被认为是相近的。它会倾向于得出直径（14.44）很小的紧致簇。不过它得到的簇可能会违背“相近”的性质。也就是说，一个簇中的元素与其他簇中的元素可能要比自己簇中的其他元素更相近。

组平均 14.43 则是单链接和全链接两个极端的折中。它试图让得出的簇既是相对紧致的同时又相对彼此远离。不过它的结果依赖于观测点之间差异性 $d_{i{i}^{\prime }}$ 的测量数量尺度。如果对 $d_{i{i}^{\prime }}$ 进行一个严格单调递增的变换 $h(\cdot )$，$h_{i{i}^{\prime }}=h(d_{i{i}^{\prime }})$，那么使用 14.43 得出的结果可能会改变。与之相反，单链接（14.41）和全链接（14.42）只依赖于 $d_{i{i}^{\prime }}$ 的顺序，因此对这样的单调变换有不变性。这种不变性有时是对比组平均方法而偏好单链接或全链接的一个原因。

有观点认为组平均聚类具有统计学上一致性的性质，而单链接和全链接则不具备。假设数据的属性变量为 $X^T=(X_1,\dots ,X_p)$，每个簇 $k$ 是某个总体联合概率密度 $p_k(x)$ 的一个随机抽样。整个数据集就是 $K$ 个这种概率密度的混合分布的一个随机样本。组平均差异性 $d_{GA}(G,H)$（14.43）是下式的估计量：

其中的 $d(x,x^{\prime })$ 为在属性空间上两个点 $x$ 和 $x^{\prime }$ 之间的差异性。随着样本大小 $N$ 趋向于无穷大，$d_{GA}(G,H)$（式 14.43）趋向于式 14.45，这是两个密度函数 $p_G(x)$ 和 $p_H(x)$ 关系的性质。对于单链接，$d_{SL}(G,H)$（式 14.41）随着 $N\to \mathrm{\infty }$ 会趋向于零，这对任意的 $p_G(x)$ 和 $p_H(x)$ 都成立。对于全链接，$d_{CL}(G,H)$（14.42）随着 $N\to \mathrm{\infty }$ 会趋向于无穷大，这也对任意的密度函数都成立。所以无法确定 $d_{SL}(G,H)$ 和 $d_{SC}(G,H)$ 是对总体概率分布的哪种统计量的估计。

例：人类肿瘤微阵列数据（续）

图 14.13 中左图展示了微阵列样本数据上使用了平均链接的凝聚聚类的树状图结果。中图和右图分别为使用了全链接和单链接。平均链接和全链接给出了相似的结果，而单链接得到了不平衡的组、细长的簇。以下讨论聚焦在平均链接聚类上。

和 K 均值聚类一样，层次聚类成功将简单的癌症类型聚类在一起。但它还有其他一些不错的特性。首先，通过在不同高度截断树状图可以得出不同数量的簇，并且这些簇存在嵌套关系。其次，它可以给出关于样本的一些部分排序信息。在图 14.14 中，我们按层次聚类中得出的顺序排列了基因（行）和样本（列）的表达矩阵。

注意若在每次合并时调整树状图中分支的定位，得到的树状图依然与层次聚类的计算过程是一致的。所以为了确定叶节点的排序，必须（在合并时）要添加一个约束。为生成图 14.14 的行排序，我们使用了 S-PLUS 中的默认规则：在每次合并时，更紧凑的簇被排放在左侧（图中旋转后的树状图则为下侧）。所有单个基因的簇是最紧凑的簇，在对两个单个基因簇合并时，按它们的观测点序号排序。在列的排序中也使用了一样的规则。也有其他可能的规则，例如根据基因的多维尺度排序，可参考第 14.8 节。

图 14.14 中双向的重排列得出了一个信息丰富的基因和样本图片。这个图比第一章图 1.3 中的行与列随机排列提供了更多信息。此外，树状图本身也有使用价值，例如生物学家可以从生物过程的角度解释基因聚类。

分裂聚类

分裂聚类（divisive clustering）算法从所有数据集作为一个单独的簇开始，自上而下地在每次迭代中将一个已有的簇分割成两个子簇。这个方法在聚类研究领域中的研究远没有凝聚方法深入和广泛。而它在工程领域中的压缩应用场景上有一些探索（Gersho and Gray, 1992）。在聚类问题中，当目标是聚焦于将数据分割成相对少的簇时，分裂方法比凝聚方法可能存在潜在的优势。

分裂方法可以通过迭代地应用任意一种组合方法来进行分割，例如令 $K=2$ 的 K 均值（第 14.3.6 节）或 K 中心点（第 14.3.10 节）。然而这样的方法会依赖与每一步中指定的初始化配置。此外，它分裂出的簇序列不一定会具备树状图所需要的单调性质。

Macnaughton Smith et al. (1965) 提出了一个可以避免这些问题的分裂算法。它首先将所有的观测点放在一个簇 $G$ 中。然后从中选择一个与所有其他观测点的平均差异性最大的观测点。这个观测点就成为了第二个簇 $H$ 中的第一个对象。在后续的每一步中，从 $G$ 中选出一个观测点，使它与 $G$ 中其余观测点的平均距离减去它与 $H$ 中其他剩余观测点的平均距离最大，然后将这个观测点转移至 $H$。重复这个操作，直到其中的平均距离之差为负数。也就是说，$G$ 中不再有观测点平均来说距离 $H$ 比 $G$ 更近了。结果就是将原来的簇分成了两个子簇，被转移到 $H$ 的观测点，以及 $G$ 中剩余的观测点。这两个簇就形成了层次结构的第二级。后续的每一级都是对前一级中的某个簇进行这样的分割过程。Kaufman and Rousseeuw (1990) 建议选择每一级中直径（式 14.44）最大的簇进行分割。另一个方式是选择内部之间平均差异性最大的簇进行分割：

重复这个迭代分割操作，直到所有的簇或者成为单元素簇，或者每个簇内的所有元素彼此间差异性都为零。

本节练习

练习 14.1

Weights for clustering

Show that weighted Euclidean distance

satisfies

where

Thus weighted Euclidean distance based on x is equivalent to unweighted Euclidean distance based on z.

练习 14.2

Consider a mixture model density in p-dimensional feature space,

where $g_k=N({\mu }_k,\mathbf{L}\cdot {\sigma }^2)$ and ${\pi }_k\ge 0$ $\mathrm{\forall }k$ with $\sum _k{\pi }_k=1$. Here $\{{\mu }_k,{\pi }_k\}$, $k=1,\dots ,K$ and ${\sigma }^2$ are unknown parameters.

Suppose we have data $x_1$, $x_2$, …, $x_N\sim g(x)$ and we wish to fit the mixture model.

1. Write down the log-likelihood of the data
2. Derive an EM algorithm for computing the maximum likelihood estimates (see Section 8.1).
3. Show that if σ has a known value in the mixture model and we take $\sigma \to 0$, then in a sense this EM algorithm coincides with K-means clustering.