一、这个题目长什么样

输入是一堆点，输出是一个编号：第 1 个点属于第 2 簇，第 2 个点属于第 0 簇…… 簇的名字没有含义，簇的个数 K 通常要你自己定。 注意「簇」这个词没有数学定义。 同一份数据，按「挨得近」分是一种分法，按「密度高」分是另一种， 按「能被一条直线分开」分又是另一种。 每种分法对应一族算法，选哪族等于先声明「我认为什么算一簇」。 K 均值的声明是：一簇就是一个球，球心是簇内点的平均，球的大小尽量小。 这个声明很强，好处是算起来极快，坏处是它看不见形状。14.2 和 14.3 两节都是在放宽这个声明—— 一个放宽「球」的形状（核），一个放宽「中心」这个概念（图）。 最后说清楚评价指标。ARI 的取值范围大致是 0 到 1，1 表示两个划分完全一致， 0 表示和随机划分一样好（可能略微为负）。 它不是准确率：簇的编号可以任意置换，所以不能逐点比对标签。

二、算法只有两步

第一步，给定 K 个中心，把每个点分给最近的那个中心。 第二步，给定每个点的归属，把每个中心挪到它那群点的平均位置。 来回做，直到中心不再动。 就这么两步，没有学习率、没有迭代步长、没有梯度。 因为每一步都是「给定对方求自己最优」，这类做法叫交替优化，和第 5 章坐标下降是同一个思路。 为什么每一步都让目标函数下降？ 第一步：中心固定时，把点分给最近的中心，正是让「点到所属中心的距离平方」最小的分法。 第二步：归属固定时，让「簇内平方和」最小的中心位置就是平均值——这是求导令零的直接结果。 所以目标函数单调不增，而且只有有限个划分，一定会停。 但「会停」和「停对」是两件事。 停下来的状态叫局部最优：每个点都已经属于最近的中心，每个中心都已经是簇的平均， 再动任何一步都会变差。 问题在于这样的状态往往不止一个，落在哪一个由初始中心决定。 实测的迭代轨迹（k-means++ 初始化，三个团）： 1527.26、445.25、356.85、356.55、356.55，五步就停。 前三步降得很凶，第四步只降了 0.30， 第五步一点没动——实际工程里通常在「变化小于某个阈值」时就提前收工。

三、初始化决定成败

同一份数据、同一个 K=3，只换随机种子，跑十个：

四、K 该取几：目标函数不能当标准

K 是必须由人给的东西，那能不能「从小到大试，取目标函数最小的」？ 看实测数据就知道不行：

五、它擅长什么，怕什么

K 均值擅长的是：球状、大小相近、彼此隔开的簇。 这种数据上它又快又准，而且结果稳定、容易解释——中心就是簇的「代表点」。 它怕的东西有三个，都来自「球」这个声明： 第一，形状不是球的。 同心圆、月牙、长条，这三种它都处理不了——因为它只能用「到中心的距离」做判断， 而这类形状根本没有一个像样的中心。 14.3 节会实测：两个月牙上它的 ARI 只有 0.3741。 第二，簇的大小差很多。 一个大簇和一个小簇挨着时，分界面会被大簇「推」向小簇， 因为界面的位置由「到两个中心的距离相等」决定，而中心位置受点数影响。 第三，离群点。 中心是平均位置，一个离得很远的点会把中心拽过去一大截，连带把整个分界面推歪。 这也是为什么第 6 章的鲁棒思想在这里同样适用：把平均换成中位数， 就得到 K 中位数（K-medians），抗离群点能力强得多。 一句话总结：K 均值不是「聚类算法」，它是「用球去铺满数据」的算法。 数据能被球铺好，它就好用；铺不好，换 14.2 或 14.3。

六、动手 Demo：初始化决定成败，目标函数不能用来选 K

这段代码在三个高斯团上跑三组实验：换初始化看结果有多不稳、看目标函数怎么一步步降、 再从小到大试 K 看目标函数和 ARI 怎么分道扬镳。 重点看两处：随机初始化那一行的「目标值最差 1307.99」和 k-means++ 的 356.55； 以及 K 从 1 到 7 时目标函数单调降、ARI 却在 K=3 之后掉头。

import numpy as np

def d2(A, B):
    return np.maximum(((A[:, None, :] - B[None, :, :]) ** 2).sum(-1), 0.0)

def blobs(n=600, seed=3, std=0.55, sep=3.2):
    rg = np.random.default_rng(seed)
    cs = np.array([[0.0, 0.0], [sep, 0.0], [sep / 2.0, sep * 0.87]])
    per = n // 3
    X, y = [], []
    for i, c in enumerate(cs):
        X.append(c + rg.normal(0.0, std, (per, 2)))
        y.append(np.full(per, i))
    return np.vstack(X), np.concatenate(y)

def pp_init(X, K, rg):     # k-means++ 初始化
    n = len(X)
    C = [X[rg.integers(n)]]
    for _ in range(K - 1):
        d = d2(X, np.array(C)).min(axis=1)
        C.append(X[rg.choice(n, p=d / d.sum())])
    return np.array(C)

def kmeans(X, K, seed=0, iters=100, plusplus=False, ret_path=False):
    rg = np.random.default_rng(seed)
    C = (pp_init(X, K, rg) if plusplus
         else X[rg.choice(len(X), K, replace=False)]).astype(float).copy()
    path = [float(((X - C[np.argmin(d2(X, C), axis=1)]) ** 2).sum())]
    for _ in range(iters):
        lab = np.argmin(d2(X, C), axis=1)
        newC = C.copy()
        for k in range(K):
            m = lab == k
            if m.any():
                newC[k] = X[m].mean(axis=0)
        path.append(float(((X - newC[np.argmin(d2(X, newC), axis=1)]) ** 2).sum()))
        if np.allclose(newC, C):
            C = newC
            break
        C = newC
    lab = np.argmin(d2(X, C), axis=1)
    return C, lab, path if ret_path else None

def ari(a, b):
    ua, ub = np.unique(a), np.unique(b)
    tab = np.zeros((len(ua), len(ub)))
    for i, x in enumerate(ua):
        for j, y in enumerate(ub):
            tab[i, j] = np.sum((a == x) & (b == y))
    c2 = lambda z: z * (z - 1.0) / 2.0
    sij, sa, sb = c2(tab).sum(), c2(tab.sum(1)).sum(), c2(tab.sum(0)).sum()
    exp, mx = sa * sb / c2(float(len(a))), 0.5 * (sa + sb)
    return float((sij - exp) / (mx - exp))

X, y = blobs(600, 3)
for tag, pp in (('随机', False), ('++', True)):
    a, it = [], []
    for sd in range(10):
        C, lab, _ = kmeans(X, 3, seed=sd, plusplus=pp)
        a.append(ari(lab, y))
        it.append(float(((X - C[lab]) ** 2).sum()))
    print('%s ARI 最小 %.3f 最大 %.3f 目标最差 %.2f 最好 %.2f'
          % (tag, min(a), max(a), max(it), min(it)))

C, lab, path = kmeans(X, 3, seed=0, plusplus=True, ret_path=True)
print('迭代轨迹', np.round(path, 2))

for K in range(1, 8):
    C, lab, _ = kmeans(X, K, seed=0, plusplus=True)
    print('K=%d 簇内平方和 %.2f ARI %.4f'
          % (K, float(((X - C[lab]) ** 2).sum()), ari(lab, y)))

实跑输出：

--------------------------------------------------------------------
K 均值：换个初始化，ARI 从 0.43 跳到 0.99，目标值还更低
--------------------------------------------------------------------
  数据                     600 个二维点，三个各向同性的高斯团，摆成正三角形
  真值                     簇内平方和 356.66

  先看初始化有多要命。同一个数据、同一个 K=3，只换随机种子：
  初始化          ARI 最小    中位    最大   目标值最差     最好
  -------------------------------------------------------------
  随机挑 3 个点      0.432   0.985   0.985      1307.99   356.55
  k-means++          0.985   0.985   0.985       356.55   356.55
  **随机初始化的十个种子里，ARI 从 0.432 到 0.985——它有时候能找到三个团，有时候只能找到两个半。**
  **k-means++ 把最差的那一次从 0.432 抬到 0.985。**
  为什么？随机挑三个点当中心，很容易两个点落在同一个团里。
  **一旦两个中心都在同一个团，迭代就不容易出来了——因为每个点只往最近的中心跑。**
  k-means++ 的做法是：第一个中心随机，之后每个中心按「离已有中心的平方距离」抽，
  **离得越远越容易被选中，于是三个中心自然散开。**

  再看迭代过程。目标函数（簇内平方和）每一步都在降：
  每一步的目标值           1527.26  445.25  356.85  356.55  356.55
  步数                     5 步（含初值那一步）
  **它是单调下降的，这一点可以证明：先固定中心重新分配，目标降；再固定分配重算中心，目标还降。**
  **所以 K 均值一定停在一个「再也动不了」的状态——但这不等于停在对的地方。**
  **它停在局部最优，而局部最优有几个，取决于初始化。**

  现在换一个常见的错误做法：拿目标函数去选 K。
  K   簇内平方和      ARI
  ----------------------
  1      2356.92   0.0000
  2      1324.34   0.5608
  3       356.55   0.9850
  4       312.21   0.8533
  5       283.64   0.8097
  6       240.50   0.6636
  7       201.70   0.5090
  **簇内平方和一路从 2356.92 降到 201.70，永远在降，永远不会回升。**
  **因为 K 每加一，就多给模型一个中心，它总能找到一点点「用得上」的地方——极端情况是每个点自己一个簇，目标值直接归零。**
  **所以「目标函数最小」这个标准在聚类里根本不能用。**
  **ARI 才是真正该看的：它在 K=3 达到 0.9850，然后一路掉。**
  注意 ARI 只在有标签时才算得出来。真实场景里没有标签，
  **于是问题变成：没有答案的时候，怎么知道 K 该取几？这是 14.4 节的事。**

第一张表的第一行就是这次实验的核心：随机初始化十个种子，ARI 从 0.432 到 0.985， 目标值从 1307.99 到 356.55。 同一份数据、同一个算法、同一个 K，只因为一开始挑的三个点不同，结果差四倍。 第二行 k-means++ 全是 0.985、目标值全是 356.55，十个种子一个不差。 这说明「初始化」不是玄学，它就是这一节的胜负手。 注意 k-means++ 的最好目标值 356.55 比真划分的 356.66 还略低一点点。 这不是算错了：真划分是按「谁生成的点」定的，而 K 均值的最小值是「按簇内平方和」定的， 两者不是同一个东西。 所以「和真值一致」和「目标函数最小」并不always等同——这个区别在第 14.4 节会变得更刺眼。 迭代轨迹那一行是 1527.26、445.25、356.85、356.55、356.55。 它严格不增，这就是交替优化的保证；但停下时只说明「动不了了」，不说明「停对了」。 最后七行是这一节的第二个落点： 簇内平方和从 2356.92 单调降到 201.70，而 ARI 在 K=3 达到 0.9850 后一路掉到 0.5090。 一个永远在降、一个先升后降。选 K 的时候该看哪个，一目了然。

30 秒实验

1. 把 blobs 的 std 从 0.55 改成 1.2，三个团会重叠起来。这时 K 均值的 ARI 会明显下降， 而且 K=3 和 K=4 的目标函数差距会变小——簇的边界本来就模糊了。
2. 把 kmeans 里的 iters 改成 1，跑一次就停。 看 ARI 掉到多少，再想想为什么实际工程里从不会这么干。
3. 把随机初始化的种子数从 10 加到 50，统计「ARI 低于 0.9 的比例」。 这个比例就是「单次运行失败率」，它决定了工程上该跑几次。

常见误区

1. 「跑出来的目标函数更小，说明聚类更好」。 实测随机初始化最差的一次目标值 1307.99，最好的一次 356.55， 而两者的划分质量差别巨大——但在不知道答案的时候，你只看到目标值，看不到 ARI。 目标值小只说明「这批中心配这个划分比较紧」，不说明「这个划分是真的」。
2. 「K 均值一定能找到最优解，因为有收敛保证」。 收敛保证的是「一定停」，不是「停在对的地方」。 它停在局部最优，而局部最优有几个取决于初始化——这一点和梯度下降没有区别。
3. 「先用 K 均值跑一遍，再用簇内平方和选 K」。 簇内平方和对 K 单调，永远选最大的那一个，等于什么都没选。要选 K， 得用带复杂度惩罚的判据（14.4 节），或者用轮廓系数这类内部指标。
4. 「K 均值对数据尺度不敏感」。它算的是欧氏距离， 哪个坐标的数值范围大，哪个坐标就主导分簇——这和 13.2 节 PCA 的问题一模一样。 不同量纲的坐标必须先标准化。

K 均值只有一个声明：一簇就是一个球。下一节不换算法，只换空间——把数据先映射到一个高维空间， 在那里球是球，回到原始空间就成了任意弯的形状。

• ☐ （本节完，下一节 14.2 14.2）

相关笔记

• 核 K 均值聚类
• 谱聚类
• 核函数主成分分析
• 🤖 图解机器学习 · 索引