k-means简介

k-means 算法在不带标签的多维数据集中寻找确定数量的簇。最优的聚类结果需要符合以
下两个假设。

• “簇中心点”（cluster center）是属于该簇的所有数据点坐标的算术平均值。
• 一个簇的每个点到该簇中心点的距离，比到其他簇中心点的距离短。

这两个假设是k-means 模型的基础，后面会具体介绍如何用该算法解决问题。先通过一个
简单的数据集，看看k-means 算法的处理结果。
首先，生成一个二维数据集，该数据集包含4 个明显的簇。由于要演示无监督算法，因此
去除可视化图中的标签：

from sklearn.datasets.samples_generator import make_blobs
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()

X, y_true = make_blobs(n_samples=300, centers=4,
        cluster_std=0.60, random_state=0)
plt.scatter(X[:, 0], X[:, 1], s=50);

E:\anaconda3\lib\site-packages\sklearn\utils\deprecation.py:143: FutureWarning: The sklearn.datasets.samples_generator module is  deprecated in version 0.22 and will be removed in version 0.24. The corresponding classes / functions should instead be imported from sklearn.datasets. Anything that cannot be imported from sklearn.datasets is now part of the private API.
  warnings.warn(message, FutureWarning)

通过肉眼观察，可以很轻松地挑选出4 个簇。而k-means 算法可以自动完成4 个簇的识别
工作，并且在Scikit-Learn 中使用通用的评估器API：

from sklearn.cluster import KMeans
kmeans = KMeans(n_clusters=4)
kmeans.fit(X)
y_kmeans = kmeans.predict(X)

下面用带彩色标签的数据来展示聚类结果。同时，画出簇中心点，这些簇中心点是由
k-means 评估器确定的：

plt.scatter(X[:, 0], X[:, 1], c=y_kmeans, s=50, cmap='viridis')
centers = kmeans.cluster_centers_
plt.scatter(centers[:, 0], centers[:, 1], c='black', s=200, alpha=0.5);

告诉你个好消息，k-means 算法可以（至少在这个简单的例子中）将点指定到某一个类，
就类似于通过肉眼观察然后将点指定到某个类。但你可能会好奇，这些算法究竟是如何快
速找到这些簇的，毕竟可能存在的簇分配组合方案会随着数据点的增长而呈现指数级增
长趋势，如果要做这样的穷举搜索需要消耗大量时间。好在有算法可以避免这种穷举搜
索——k-means 方法使用了一种容易理解、便于重复的期望最大化算法取代了穷举搜索。

k-means算法：期望最大化

期望最大化（expectation-maximization，E-M）是一种非常强大的算法，应用于数据科学
的很多场景中。k-means 是该算法的一个非常简单并且易于理解的应用，下面将简单介绍
E-M 算法。简单来说，期望最大化方法包含以下步骤：

(1) 猜测一些簇中心点。

(2) 重复直至收敛。

• 期望步骤（E-step）：将点分配至离其最近的簇中心点。
• 最大化步骤（M-step）：将簇中心点设置为所有点坐标的平均值。

期望步骤（E-step 或Expectation step）不断更新每个点是属于哪一个簇的期望值，最大
化步骤（M-step 或Maximization step）计算关于簇中心点的拟合函数值最大化对应坐标
（argmax 函数）——在本例中，通过简单地求每个簇中所有数据点坐标的平均值得到了簇
中心点坐标。
关于这个算法的资料非常多，但是这些资料都可以总结为：在典型环境下，每一次重复
E-step 和M-step 都将会得到更好的聚类效果。
将这个算法在图5-112 中可视化。[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-9JuKlaxM-1638085197886)(attachment:%E5%9B%BE%E7%89%87.png)]
如图所示，数据从初始化状态开始，经过三次迭代后收敛。下图显示的是聚类的交互式可视
化版本，详情请参见GitHub 在线附录（https://github.com/jakevdp/PythonDataScienceHandbook）
中的代码。
k-means 算法非常简单，只要用几行代码就可以实现它。以下是一个非常基础的k-means
算法实现：

from sklearn.metrics import pairwise_distances_argmin
def find_clusters(X, n_clusters, rseed=2):
    # 1.随机选择簇中心点
    rng = np.random.RandomState(rseed)
    i = rng.permutation(X.shape[0])[:n_clusters]
    centers = X[i]
    while True:
        # 2a.基于最近的中心指定标签
        labels = pairwise_distances_argmin(X, centers)
        # 2b.根据点的平均值找到新的中心
        new_centers = np.array([X[labels == i].mean(0)
                        for i in range(n_clusters)])
        # 2c.确认收敛
        if np.all(centers == new_centers):
            break
        centers = new_centers
    return centers, labels
centers, labels = find_clusters(X, 4)
plt.scatter(X[:, 0], X[:, 1], c=labels,
s=50, cmap='viridis');

虽然大部分可用的聚类算法底层其实都是对上述示例的进一步扩展，但上述函数解释了期
望最大化方法的核心内容。
使用期望最大化算法时的注意事项
在使用期望最大化算法时，需要注意几个问题。

可能不会达到全局最优结果

• 首先，虽然E–M 算法可以在每一步中改进结果，但是它并不保证可以获得全局最优的解决方案。例如，如果在上述简单的步骤中使用一个随机种子（random seed），那么某些初始值可能会导致很糟糕的聚类结果：

centers, labels = find_clusters(X, 4, rseed=0)
plt.scatter(X[:, 0], X[:, 1], c=labels,
s=50, cmap='viridis');

虽然E–M 算法最终收敛了，但是并没有收敛至全局最优配置。因此，该算法通常会用
不同的初始值尝试很多遍，在Scikit-Learn 中通过n_init 参数（默认值是10）设置执
行次数。

簇数量必须事先定好

k-means 还有一个显著的问题：你必须告诉该算法簇数量，因为它无法从数据中自动学
习到簇的数量。如果我们告诉算法识别出6 个簇，它将很快乐地执行，并找出最佳的6
个簇：

labels = KMeans(6, random_state=0).fit_predict(X)
plt.scatter(X[:, 0], X[:, 1], c=labels,
s=50, cmap='viridis');

结果是否有意义是一个很难给出明确回答的问题。有一个非常直观的方法，但这里不
会进一步讨论，该方法叫作轮廓分析（http://scikit-learn.org/stable/auto_examples/cluster/
plot_kmeans_silhouette_analysis.html）。

不过，你也可以使用一些复杂的聚类算法，有些算法对每个簇的聚类效果有更好的度
量方式（例如高斯混合模型，Gaussian mixture models，详情请参见5.12 节），还有一些
算法可以选择一个合适的簇数量（例如DBSCAN、均值漂移或者近邻传播，这些都是
sklearn.cluster 的子模块）。

k-means 算法只能确定线性聚类边界

k-means 的基本模型假设（与其他簇的点相比，数据点更接近自己的簇中心点）表明，
当簇中心点呈现非线性的复杂形状时，该算法通常不起作用。

k-means 聚类的边界总是线性的，这就意味着当边界很复杂时，算法会失效。用下面的
数据来演示k-means 算法得到的簇标签，如图5-116 所示：

from sklearn.datasets import make_moons
X, y = make_moons(200, noise=.05, random_state=0)
labels = KMeans(2, random_state=0).fit_predict(X)
plt.scatter(X[:, 0], X[:, 1], c=labels,
s=50, cmap='viridis');

这个情形让人想起5.7 节介绍的内容，当时我们通过一个核变换将数据投影到更高维的
空间，投影后的数据使线性分离成为可能。或许可以使用同样的技巧解决k-means 算法
无法处理非线性边界的问题。

这种核k-means 算法在Scikit-Learn 的SpectralClustering 评估器中实现，它使用最近
邻图（the graph of nearest neighbors）来计算数据的高维表示，然后用k-means 算法分配
标签：

from sklearn.cluster import SpectralClustering
model = SpectralClustering(n_clusters=2,
affinity='nearest_neighbors',
assign_labels='kmeans')
labels = model.fit_predict(X)
plt.scatter(X[:, 0], X[:, 1], c=labels,
s=50, cmap='viridis');

E:\anaconda3\lib\site-packages\sklearn\manifold\_spectral_embedding.py:236: UserWarning: Graph is not fully connected, spectral embedding may not work as expected.
  warnings.warn("Graph is not fully connected, spectral embedding"

可以看到，通过核变换方法，核k-means 就能够找到簇之间复杂的非线性边界了。

当数据量较大时，k-means 会很慢

由于k-means 的每次迭代都必须获取数据集所有的点，因此随着数据量的增加，算法
会变得缓慢。你可能会想到将“每次迭代都必须使用所有数据点”这个条件放宽，例
如每一步仅使用数据集的一个子集来更新簇中心点。这恰恰就是批处理（batch-based）
k-means 算法的核心思想，该算法在sklearn.cluster.MiniBatchKMeans 中实现。该算法
的接口和标准的KMeans 接口相同，后面将用一个示例来演示它的用法[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-qchN4tqu-1638085197898)(attachment:%E5%9B%BE%E7%89%87.png)]

案例

了解了算法的限制条件之后，就可以将k-means 的优势发挥到合适的场景中了。下面来演
示一些示例。

1. 案例1：用k-means算法处理手写数字
   

首先，将k-means 算法用于5.8 节和5.9 节演示的手写数字数据。这次试试能不能不使用原
始的标签信息，就用k-means 算法识别出类似的数字。这个过程就好像是在事先没有标签
信息的情况下，探索新数据集的含义。

首先导入手写数字，再使用KMeans 聚类。手写数字数据集包含1797 个示例，每个样本有
64 个特征，其实就是8×8 图像中的每个像素：

from sklearn.datasets import load_digits
digits = load_digits()
digits.data.shape

(1797, 64)

kmeans = KMeans(n_clusters=10, random_state=0)
clusters = kmeans.fit_predict(digits.data)
kmeans.cluster_centers_.shape

(10, 64)

结果是在64 维中有10 个类。需要注意的是，这些簇中心点本身就是64 维像素的
点，可以将这些点看成是该簇中“具有代表性”（typical）的数字。这些簇中心点：

fig, ax = plt.subplots(2, 5, figsize=(8, 3))
centers = kmeans.cluster_centers_.reshape(10, 8, 8)
for axi, center in zip(ax.flat, centers):
    axi.set(xticks=[], yticks=[])
    axi.imshow(center, interpolation='nearest', cmap=plt.cm.binary)

我们会发现，即使没有标签，KMeans 算法也可以找到可辨识的数字中心，但是1 和8
例外。

由于k-means 并不知道簇的真实标签，因此0~9 标签可能并不是顺序排列的。我们可以将
每个学习到的簇标签和真实标签进行匹配，从而解决这个问题：

from scipy.stats import mode
labels = np.zeros_like(clusters)
for i in range(10):
    mask = (clusters == i)
    labels[mask] = mode(digits.target[mask])[0]

from sklearn.metrics import accuracy_score
accuracy_score(digits.target, labels)

0.7952142459654981

from sklearn.metrics import confusion_matrix
mat = confusion_matrix(digits.target, labels)
sns.heatmap(mat.T, square=True, annot=True, fmt='d', cbar=False,
xticklabels=digits.target_names,
yticklabels=digits.target_names)
plt.xlabel('true label')
plt.ylabel('predicted label');

正如我们之前看到的簇中心点图，混淆的地方主要是在数字8 和1。但仍然可以看出，通
过k-means 可以构建一个数字分类器，该数字分类器不需要任何已知的标签。
其实还可以更进一步，使用t- 分布邻域嵌入算法（详情请参见5.10 节）在执行k-means 之
前对数据进行预处理。t-SNE 是一个非线性嵌入算法，特别擅长保留簇中的数据点。下面
来看看如何实现：

from sklearn.manifold import TSNE
# 投影数据：这一步将耽误几秒钟
tsne = TSNE(n_components=2, init='pca', random_state=0)
digits_proj = tsne.fit_transform(digits.data)
# 计算类
kmeans = KMeans(n_clusters=10, random_state=0)
clusters = kmeans.fit_predict(digits_proj)
# 排列标签
labels = np.zeros_like(clusters)
for i in range(10):
    mask = (clusters == i)
    labels[mask] = mode(digits.target[mask])[0]
# 计算准确度
accuracy_score(digits.target, labels)

0.9398998330550918

同样在没有标签的情况下，它可以达到94% 的分类准确率，这就是合理使用无监督学习的
力量。无监督学习可以从数据集中抽取难以用手眼直接提取的信息。

2. 案例2：将k-means用于色彩压缩

聚类算法的另一个有趣应用是图像色彩压缩。设想你有一幅包含几百万种颜色的图像，但
其实大多数图像中的很大一部分色彩通常是不会被眼睛注意到的，而且图像中的很多像素
都拥有类似或者相同的颜色。

如图下所示，该图像来源于Scikit-Learn 的datasets 模块（在本例中，你将需要安装
Python 的pillow 图像程序包）：

from sklearn.datasets import load_sample_image
china = load_sample_image("china.jpg")
ax = plt.axes(xticks=[], yticks=[])
ax.imshow(china);

china.shape

(427, 640, 3)

可以将这组像素转换成三维颜色空间中的一群数据点。先将数据变形为[n_samples ×
n_features]，然后缩放颜色至其取值为0~1：

data = china / 255.0 # 转换成0~1区间值
data = data.reshape(427 * 640, 3)
data.shape

(273280, 3)

还可以在颜色空间中对这些像素进行可视化。为了演示方便，这里只使用包含前10 000 个
像素的子集：

def plot_pixels(data, title, colors=None, N=10000):
    if colors is None:
        colors = data
        # 随机选择一个子集
        rng = np.random.RandomState(0)
        i = rng.permutation(data.shape[0])[:N]
        colors = colors[i]
        R, G, B = data[i].T
        fig, ax = plt.subplots(1, 2, figsize=(16, 6))
        ax[0].scatter(R, G, color=colors, marker='.')
        ax[0].set(xlabel='Red', ylabel='Green', xlim=(0, 1), ylim=(0, 1))
        ax[1].scatter(R, B, color=colors, marker='.')
        ax[1].set(xlabel='Red', ylabel='Blue', xlim=(0, 1), ylim=(0, 1))
        fig.suptitle(title, size=20);

plot_pixels(data, title='Input color space: 16 million possible colors')

现在对像素空间（特征矩阵）使用k-means 聚类，将1600 万种颜色（255 × 255 × 255 =
16 581 375）缩减到16 种颜色。因为我们处理的是一个非常大的数据集，所以将使用
MiniBatchKMeans 算法对数据集的子集进行计算。这种算法比标准的k-means 算法速度更快

from sklearn.cluster import MiniBatchKMeans
kmeans = MiniBatchKMeans(16)
kmeans.fit(data)
new_colors = kmeans.cluster_centers_[kmeans.predict(data)]
plot_pixels(data, colors=new_colors,
title="Reduced color space: 16 colors")

用计算的结果对原始像素重新着色，即每个像素被指定为距离其最近的簇中心点的颜色。
用新的颜色在图像空间（427 × 640），而不是像素空间（273 280 × 3）里重新画图，展示
效果：

china_recolored = new_colors.reshape(china.shape)
fig, ax = plt.subplots(1, 2, figsize=(16, 6),
subplot_kw=dict(xticks=[], yticks=[]))
fig.subplots_adjust(wspace=0.05)
ax[0].imshow(china)
ax[0].set_title('Original Image', size=16)
ax[1].imshow(china_recolored)
ax[1].set_title('16-color Image', size=16);

虽然右图丢失了某些细节，但是图像总体上还是非常容易辨识的。右图实现了将近一百万
的压缩比！这就是k-means 的一个有趣的应用，当然还有很多更好的压缩图像的算法，但
是这个示例足以显示无监督算法（如k-means）解决问题的力量

另外，在平时的使用过程中，聚类只是一个小小的数据处理方式，并没有那么复杂。在确定簇类数目有一个小窍门。使用SSE与轮廓系数确定：

from sklearn.cluster import KMeans
SSE = []  # 存放每次结果的误差平方和，越小越好
for k in range(1,9):
    estimator = KMeans(n_clusters=k)  # 构造聚类器
    estimator.fit(x)
    SSE.append(estimator.inertia_)
X1 = range(1,9)
plt.style.use('Solarize_Light2')
plt.xlabel('k')
plt.ylabel('SSE')
plt.plot(X1,SSE,'o-',c='y')
plt.savefig('误差平方和.pdf')


from sklearn.metrics import silhouette_score

Scores = []  # 存放轮廓系数,越大越好
for k in range(2,9):
    estimator = KMeans(n_clusters=k)  # 构造聚类器
    estimator.fit(x)
    Scores.append(silhouette_score(x,estimator.labels_,metric='euclidean'))
X = range(2,9)
plt.style.use('Solarize_Light2')
plt.xlabel('k')
plt.ylabel('Silhouette Coefficient')
plt.plot(X,Scores,'o-',c='y')
# plt.savefig('轮廓系数.pdf')