目录

TSNE + K-Means:可视化与聚类

Info

每个样本只能属于一个类,因此 K-Means 算法是硬聚类算法

最近在工作中做数据分析又用到了 K-Means 算法,忽然发现自己记录下的 K-Means 算法笔记,于是进行了一番整理发了出来 :)

K-Means 算法是一种聚类算法,即给定 $n$ 个样本的集合 $X={\mathbf x_1,\mathbf x_2,…,\mathbf x_n}$ 的情况下,K-Means 会将这 $n$ 个样本分到 $K$ 个不同的簇里面,每一个簇对应一个类别

因为 K-Means 算法用到了距离度量(后面会进行解释),所以需要做标准化。这也很简单,使用 Scikit-Learn 即可,代码如下

from sklearn import preprocessing
X_normalized = preprocess.normalzie(X)

输入:簇的数量(用 $K$ 表示)

算法

  1. 在数据中随机初始化 $K$ 个簇质心(cluster centroids)为 $\mu_1,\mu_2,…,\mu_k$,每一个表示不同类别数据的中心位置
  2. 重复下面的步骤直到收敛,即簇质心的位置不再显著变化
    1. 遍历数据集的每一个样本:看离 $K$ 个簇质心点哪一个最近就把这个样本分配给对应的簇质心
    2. 现在每个样本都被分配到了对应的簇里面(可能是新的,也可能是原有的),遍历每一个簇,计算新的簇质心

输出:$K$ 个簇,每个样本都会被分配到一个已有的簇上

可以看到,每一步都需要计算样本和簇质心之间的距离,因此做标准化才不会导致量纲不同的时候某些维度的特征对距离的计算影响过大

从前面的算法流程来看还是挺简单的,自己实现也不是什么难事,那么要如何度量 K-Means 算法的聚类效果呢?主要有下面 2 个指标

Inertia 指的是每个样本到所属的簇的质心的距离的平均值 1

$$ \texttt{inertia}=\frac{1}{m}\sum_m\texttt{distance}(\mathbf x_i,\mu_i) $$

其中:

  • $m$ 表示这个簇的点点数量
  • $\mu_i$ 表示样本 $\mathbf x_i$ 所属的质心
  • $\texttt{distance}$ 表示距离函数,比如欧式距离
Tip

Scikit-LearnKMeans 模型在 fit 之后就可以访问它的 inertia_ 拿到这个指标

轮廓系数用于衡量 K-Means 算法的聚类效果,它等于所有节点的轮廓系数的平均值。每个节点的轮廓系数的计算公式如下

$$ \texttt{Sihouette}(\mathbf x)=\frac{B - A}{max(A, B)} $$

其中

  • $A$ 表示点 $\mathbf x$ 到同簇的其他点的平均距离
  • $B$ 表示点 $\mathbf x$ 到相邻的簇的其他点的平均距离

理想情况是 $A$ 小 $B$ 大,极限情况下 $B > A, A=0$,那么此时轮廓系数就近似为 $B/B=1$,所以轮廓系数的最大值是 $1$。反之也可以推导出轮廓系数的最小值是 $-1$

如果你要解决的问题的 $K$ 是已知的,那么直接用就行

否则,可以用前面提到的 2 个度量指标看看大概 $K$ 等于多少的时候收敛,选取的原则是

  • Inertia 尽量低。一个比较出名的办法是看 Inertia 曲线的“肘部”法则
  • 轮廓系数尽量高
Info

你可以点击这里查看该 Demo

Note

在实际的项目中,每个样本的特征向量长度可能比 4 多得多,样本数量也多得多,但思想是类似的

这里我们用 Iris 数据集做演示,主要考虑到

  • 只有 150 个样本 => 比较容易观察聚类效果如何
  • 特征向量长度是 4 而不是 2 => 还可以顺便演示一下如何用 TSNE 做可视化

这是一个很出名的数据集,所以你可能已经知道了它实际上有多少个类别,但请你装作不知道 :)

首先导入需要的库

import pandas as pd
import altair as alt
from sklearn.datasets import load_iris
from sklearn.manifold import TSNE
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
from sklearn import preprocessing

然后导入数据并做标准化

data, target = load_iris(return_X_y=True, as_frame=True)
X, y = data.to_numpy(), target.to_numpy()
X_normalized = preprocessing.normalize(X)
df = pd.concat((pd.DataFrame(X_normalized, columns = df.columns[:-1]), target), axis=1)

可以用 df.head(5) 看一下前面 5 条数据

id sepal length (cm) sepal width (cm) petal length (cm) petal width (cm) target
0 0.803773 0.551609 0.220644 0.031521 0
1 0.828133 0.507020 0.236609 0.033801 0
2 0.805333 0.548312 0.222752 0.034269 0
3 0.800030 0.539151 0.260879 0.034784 0
4 0.790965 0.569495 0.221470 0.031639 0

因为特征向量长度是 4,无法直接做可视化,因此我们可以用 TSNE 做降维

model = TSNE(n_components=2)
X_proj = model.fit_transform(X)
Tip

注意,做 TSNE 降维的时候使用原始数据 X,而不是做了标准化后的数据

我这里选择的画图库是 altair,所以需要先将数据准备成 DataFrame 的格式

df_proj_without_label = pd.DataFrame(X_proj, columns=["x", "y"])
df_proj_with_label = pd.concat((df_proj_without_label, df['target']), axis=1)

然后就可以画图了

alt.Chart(df_proj_without_label).mark_circle(size=60).encode(
    x = "x",
    y = "y",
).interactive().properties(width=400, height=400)

图如下所示,可以看到,大概有 3 个簇的样子,所以理想的 $K$ 可能是 3 左右

接下来我们看一下 $K=3$ 的情况下 KMeans 的聚类结果

Tip

注意,训练 K-Means 模型的时候用的是做了标准化之后的数据 X_normalized

model = KMeans(n_clusters=3, random_state=42)
model.fit(X_normalized)

df_with_new_pred = pd.concat((df_proj_without_label, pd.Series(model.labels_)), axis=1)
df_with_new_pred.columns = ["x", "y", "target"]
alt.Chart(df_with_new_pred).mark_circle(size=60).encode(
    x = "x",
    y = "y",
    color='target:N',
).interactive().properties(width=400, height=400)

可以看到,$K=3$ 的效果还可以,那么会有更好的 $K$ 吗?让我们用前面的方法找一下看有没有

ks = range(1, 12)
inertias = []
silhouette_scores = []
models = []
for k in ks:
    model = KMeans(k, random_state=42)
    model.fit(X_normalized)

    models.append(model)
    
    inertias.append(model.inertia_)
    if k != 1:
        silhouette_scores.append(silhouette_score(X, model.labels_))
        
# For plotting
silhouette_scores = [silhouette_scores[0]] + silhouette_scores
curve = pd.DataFrame(
    zip(ks, inertias, silhouette_scores), columns=["K", "inertia", "silhouette"]
)
base = alt.Chart(curve).encode(x='K')
alt.layer(
    base.mark_line(color='blue').encode(y='inertia'),
    base.mark_line(color='red').encode(y='silhouette')
)
Warning

注意:轮廓系数的计算最少需要 2 个簇,所以我这里插入了一个 dummy value

画出来的图长这个样子(注意蓝色的是 Inertia、红色的是轮廓系数)

如果看 Inertia 的话,“肘部”的位置在 $K=3$,而轮廓系数则一直在下降,所以综合来看,$K=3$ 是比较好的,而这也是 Iris 数据集的真实类别数量

那么标准答案是咋样的呢?也可以画看看

alt.Chart(df_proj_with_label).mark_circle(size=60).encode(
    x = "x",
    y = "y",
    color='target:N',
).interactive().properties(width=400, height=400)

在本文我们复习了 KMeans 算法是如何工作的,并以 Iris 数据集展示了 K-Means 算法是如何工作的。主要的几个点包括

  • 如何处理高维数据的可视化 👉 用TSNE 对原始数据做降维
  • 如何聚类 👉 用 K-Means 算法对标准化后的数据聚类
  • 如何获取比较好的 $K$ 👉 如果能提前确定则使用指定的 $K$,否则尝试不同的 $K$ 并结合 Inertia 和轮廓系数综合判断,可以用“肘部”法则

  1. 《统计学习方法》第 2 版. 李航 ↩︎