阅读须知：
- 阅读本文需要有一定的Python及Numpy基础
本文将介绍：
- K-means算法实现步骤
- 使用Python实现K-means算法
- 借助Numpy的向量计算提升计算速度
- 使用Gap Statistic法自动选取合适的聚类中心数K
聚类是一个将数据集中在某些方面相似的数据成员进行分类组织的过程，聚类就是一种发现这种内在结构的技术，聚类技术经常被称为无监督学习。
k均值聚类是最著名的划分聚类算法，由于简洁和效率使得他成为所有聚类算法中最广泛使用的。给定一个数据点集合和需要的聚类数目k，k由用户指定，k均值算法根据某个距离函数反复把数据分入k个聚类中。
K-means原理
设有样本( x 1 , x 2 , . . . , x n ) (x_1, x_2, …, x_n)(x 1 ,x 2 ,…,x n )，其中每个样本有d个特征（由d维实向量组成），K-means聚类的目标是将n n n个样本聚类至k ( ≤ n ) k(\le n)k (≤n )个集合 S = { S 1 , S 2 , . . . , S k } S = {S_1, S_2, …, S_k}S ={S 1 ,S 2 ,…,S k } 中，使簇中样本的距离和达到最小。
即用公式表示为：
arg min S ∑ i = 1 k ∑ x ∈ S i ∥ x − μ i ∥ 2 = arg min S ∑ i = 1 k ∣ S i ∣ Var S i {\displaystyle {\underset {\mathbf {S} }{\operatorname {arg\,min} }}\sum {i=1}^{k}\sum {\mathbf {x} \in S_{i}}\left\|\mathbf {x} -{\boldsymbol {\mu }}{i}\right\|^{2}={\underset {\mathbf {S} }{\operatorname {arg\,min} }}\sum {i=1}^{k}|S_{i}|\operatorname {Var} S_{i}}S a r g m i n i =1 ∑k x ∈S i ∑∥x −μi ∥2 =S a r g m i n i =1 ∑k ∣S i ∣V a r S i
其中μ i \mu i μi 是S i S_i S i 中所有点的质心，因此上式相当于最小化簇中每两点的平方距离：
arg min S ∑ i = 1 k 1 ∣ S i ∣ ∑ x , y ∈ S i ∥ x − y ∥ 2 {\displaystyle {\underset {\mathbf {S} }{\operatorname {arg\,min} }}\sum {i=1}^{k}\,{\frac {1}{|S_{i}|}}\,\sum {\mathbf {x} ,\mathbf {y} \in S{i}}\left\|\mathbf {x} -\mathbf {y} \right\|^{2}}S a r g m i n i =1 ∑k ∣S i ∣1 x ,y ∈S i ∑∥x −y ∥2
K-means步骤
K-means算法根据此从距离入手，迭代找出（局部）最小时的聚类中心S i S_i S i 。具体步骤如下：
基本实现
- 需要用聚类中心数量k进行初始化，此外可以给定迭代次数与最小误差等终止条件。
- fit函数以若干n维特征的数据作为输入。执行后，通过classifications获取分类结果；centroids获取聚类中心。
- Predict函数以一个n维特征的数据作为输入，输出归属的聚类中心索引。
- 每步操作的作用详见代码注释
class Kmeans:
def __init__(self, k=2, tolerance=0.01, max_iter=300):
self.k = k
self.tol = tolerance
self.max_iter = max_iter
self.features_count = -1
self.classifications = None
self.centroids = None
def fit(self, data):
"""
:param data: numpy数组，约定shape为：(数据数量，数据维度)
:type data: numpy.ndarray
"""
self.features_count = data.shape[1]
self.centroids = np.zeros([self.k, data.shape[1]])
for i in range(self.k):
self.centroids[i] = data[i]
for i in range(self.max_iter):
self.classifications = [[] for i in range(self.k)]
for feature_set in data:
classification = self.predict(feature_set)
self.classifications[classification].append(feature_set)
prev_centroids = np.ndarray.copy(self.centroids)
for classification in range(self.k):
self.centroids[classification] = np.average(self.classifications[classification], axis=0)
for c in range(self.k):
if np.linalg.norm(prev_centroids[c] - self.centroids[c]) > self.tol:
break
else:
return
def predict(self, data):
distances = np.linalg.norm(data - self.centroids, axis=1)
return distances.argmin()
测试代码
blobs函数产生若干个数据点云。n_features是维度，centers是中心数目，random_state是随机种子。
from sklearn.datasets import make_blobs
def blobs(n_samples=300, n_features=2, centers=1, cluster_std=0.60, random_state=0):
points, _ = make_blobs(n_samples=n_samples,
n_features=n_features,
centers=centers,
cluster_std=cluster_std,
random_state=random_state)
return points
输出聚类图象的函数，支持2D和3D，8类别以内不同颜色区分。
def kmeans_plot(kmeans_model):
"""
又长又臭的函数，简单可视化2d或3d kmeans聚类结果，不是算法必须的，直接使用即可。
:param kmeans_model: 训练的kmeans模型
:type kmeans_model: Kmeans | FastKmeans
"""
style.use('ggplot')
colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'w']
if kmeans_model.features_count == 2:
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot()
for i in range(kmeans_model.k):
color = colors[i%len(colors)]
for feature_set in kmeans_model.classifications[i]:
ax.scatter(feature_set[0], feature_set[1], marker="x", color=color, s=50, linewidths=1)
for centroid in kmeans_model.centroids:
ax.scatter(centroid[0], centroid[1], marker="o", color="k", s=50, linewidths=3)
elif kmeans_model.features_count == 3:
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(projection='3d')
for i in range(kmeans_model.k):
color = colors[i%len(colors)]
for feature_set in kmeans_model.classifications[i]:
ax.scatter(feature_set[0], feature_set[1], feature_set[2], marker="x", color=color, s=50, linewidths=1)
for centroid in kmeans_model.centroids:
ax .scatter(centroid[0], centroid[1], centroid[2], marker="o", color="k", s=50, linewidths=3)
plt.show()
测试调用
model = Kmeans(k=3)
model.fit(blobs(centers=3, random_state=1, n_features=2))
kmeans_plot(model)
model.fit(blobs(centers=3, random_state=3, n_features=3))
kmeans_plot(model)
测试结果
计算效率问题
该实现虽然可以使用，但是在进行大规模数据处理时非常慢，通过性能监控找出耗时最大的过程：
分析后发现：距离计算的函数被调用了很多次。进而说明，主要函数（如距离）的计算单位是小规模点集，因此这些函数在python层面上调用了多次。众所周知python的执行效率比较慢，特别是耗时排第一的norm函数，对于小规模计算的开销是极大的。
最好的解决办法是计算向量化：设法一次计算一批、甚至全部数据。这样就可以把计算交给python底层的c语言实现；而且向量（或矩阵）的计算本身也是被优化过的。如此速度上会提升许多。
计算优化实现
计算优化借助了numpy库中 向量广播的特性。
如下图，被减数是(N, 1, D)维向量，表示N个D维数据的数据集；减数是(1, K, D)维向量，表示K个D维的聚类中心。相减时，由于向量维度不匹配，numpy会将向量广播为矩阵（结果为虚线表示的矩阵）再相减。而该例子中，所得结果矩阵恰好为：每个点与各个聚类中心的距离信息。
如此一来，甚至不需要显式写一个循环就可以完成一整轮的距离计算。这样做把计算函数交给了受优化的矩阵运算，交给了高效的底层c实现。
- compute_l2_distance利用了numpy广播
class FastKmeans:
def __init__(self, k=2, tolerance=0.01, max_iter=300):
self.k = k
self.tol = tolerance
self.max_iter = max_iter
self.features_count = -1
self.classifications = None
self.centroids = None
def fit(self, data):
"""
:param data: numpy数组，约定shape为：(数据数量，数据维度)
:type data: numpy.ndarray
"""
self.features_count = data.shape[1]
self.centroids = np.zeros([self.k, data.shape[1]])
for i in range(self.k):
self.centroids[i] = data[i]
points_centroid = np.zeros(data.shape[0])
for i in range(self.max_iter):
prev_centroid = np.ndarray.copy(points_centroid)
points_centroid = self.get_closest_centroid(data)
for c in range(self.centroids.shape[0]):
classification = data[points_centroid == c]
self.centroids[c] = classification.mean(axis=0)
if FastKmeans.compute_l2_distance(prev_centroid, points_centroid) < self.tol:
break
self.classifications = []
for c in range(self.centroids.shape[0]):
self.classifications.append(data[points_centroid == c])
def predict(self, data):
distances = np.linalg.norm(data - self.centroids, axis=1)
return distances.argmin()
@staticmethod
def compute_l2_distance(x, centroid):
dist = ((x - centroid) ** 2).sum(axis=x.ndim - 1)
return dist
def get_closest_centroid(self, x):
dist = FastKmeans.compute_l2_distance(x[:, None, :], self.centroids[None, :, :])
closest_centroid_index = np.argmin(dist, axis=1)
return closest_centroid_index
每组对比实验分别使用两种做法对n个数据进行从1到20类的聚类。记录两种做法分别的耗时。
500100050001000020000普通实现0.862.9123.7259.17129.69计算优化实现0.070.983.846.4128.34比率12.28572.96946.17719.23094.5762
可见在不同的数据规模下均有一定提升。两种做法的速度差距在之后选择k的过程中更加明显。
常用方法
其思想是：
随着聚类数k的增大，样本划分会更加精细，每个簇的聚合程度会逐渐提高，那么误差平方和SSE自然会逐渐变小。
当k小于真实聚类数时，由于k的增大会大幅增加每个簇的聚合程度，故SSE的下降幅度会很大；而当k到达真实聚类数时，再增加k所得到的聚合程度回报会迅速变小，所以SSE的下降幅度会骤减，然后随着k值的继续增大而趋于平缓，也就是说SSE和k的关系图是一个手肘的形状，而这个肘部对应的k值就是数据的真实聚类数。
但该做法的缺点是需要人工判断”手肘”位置到底在哪，在类别多时难以判断。
其定义为：
G a p ( K ) = E ( l o g D k ) − l o g D k Gap(K)=E(logD_k)-logD_k G a p (K )=E (l o g D k )−l o g D k
其中，E ( l o g D k ) E(logD_k)E (l o g D k )是l o g D k logD_k l o g D k 的期望，一般使用蒙特卡洛模拟产生。
模拟的基本过程是：首先在样本所在区域内按照均匀分布随机地产生和原始样本数一样多的随机样本，并用K-means模型对这个随机样本聚类，计算结果的误差和，得到一个D k D_k D k 。重复多次就可以近似计算出E ( l o g D k ) E(logD_k)E (l o g D k )。
实际上，计算E ( l o g D k ) E(logD_k)E (l o g D k )只是为了给评判实际误差l o g D k logD_k l o g D k 提供一个标准。当选取合适的K K K值时，D k D_k D k 应处于偏离（远低于）平均值的极端情况，因此通过与期望误差做差，得到最大G a p ( K ) Gap(K)G a p (K )所对应的K即是最好的选择。
设聚类簇C r C_r C r 中两点间距离和为：
D r = ∑ i , i ′ ∈ C r d i i ′ D_r=\sum_{i,i’\in C_r} d_{ii’}D r =i ,i ′∈C r ∑d i i ′
若d i i ′ d_{ii’}d i i ′为欧几里得距离，则可定义”每个簇内平方和均值”的总加和W k W_k W k 为：
W k = ∑ r = 1 k 1 2 n r D r W_k=\sum_{r=1}^{k}\frac{1}{2n_r}D_r W k =r =1 ∑k 2 n r 1 D r
其中因子2 2 2恰好使上式成立；数据规模n n n被约分掉了。
该方法通过与一个合适的零分布比较来标准化log ( W k ) \log(W_k)lo g (W k )，而最优聚类数则是使log ( W k ) \log(W_k)lo g (W k )最偏离于上述参考分布时的值k k k。因此定义：
G a p n ( k ) = E n ∗ { log ( W k ) } − log ( W k ) Gap_n(k)=E_n^{\log(W_k)}-\log(W_k)G a p n (k )=E n ∗{lo g (W k )}−lo g (W k )
其中E n ∗ E_n^E n ∗表示样本规模n n n的数据在参考分布中的数学期望。考虑抽样分布后，使G a p n ( k ) Gap_n(k)G a p n (k )最大化的值就是所估计聚类数k k k的最优值。
这个方法的操作是一般化的，可以用在以任意形式计算距离d i i ′ d_{ii’}d i i ′的任意的聚类方法。
Gap Statistic的思路启发是：假设有n个、k簇、均匀分布的p维数据点，设想它们在各自的簇中均匀分布，则log ( W k ) \log(W_k)lo g (W k )的期望近似是：
log ( p n 12 ) − ( 2 p ) log ( k ) + C o n s t a n t \log(\frac{pn}{12})-(\frac{2}{p})\log(k)+Constant lo g (1 2 p n )−(p 2 )lo g (k )+C o n s t a n t
如果数据确实有K个相互分离的聚类，对于k ≤ K k\le K k ≤K，log ( W k ) \log(W_k)lo g (W k )应比预期下降率( 2 p ) log ( k ) (\frac{2}{p})\log(k)(p 2 )lo g (k )下降的更快；当k > K k \gt K k >K，事实上是在加入不必要的聚类中心，此时由简单代数可得log ( W k ) \log(W_k)lo g (W k )应下降的比其预期速率更慢。因此G a p n ( K ) Gap_n(K)G a p n (K )会在k = K k=K k =K时取得最大值。
使用Gap Statistic选取最优的k
首先，该算法需要评估聚类误差D k D_k D k ，由于最后的标准是相对的，因此D k D_k D k 的求法并无过多约束，只要能体现其误差即可，因此还是选择最简便的做法：各点到聚类中心的距离为单个误差，将其加和作为最终的误差，由sum_distance处理这一过程。
虽然聚类中心K的最优解是不依赖算法客观存在的，但由于不同K-means实现会得出不同的D k D_k D k ，因此为了G a p ( K ) Gap(K)G a p (K )具有可比性，需要借助同一种K-means实现求解误差D k D_k D k ，因此sum_distance中统一使用FastKmeans。
import scipy
def sum_distance(data, k):
model = FastKmeans(k=k)
model.fit(data)
disp = 0
for m in range(len(model.classifications)):
disp += sum(np.linalg.norm(model.classifications[m] - model.centroids[m], axis=1))
return disp
gap函数需要给定k的测试范围ks以及蒙特卡洛模拟次数nrefs，负责产生随机样品估计E ( l o g D k ) E(logD_k)E (l o g D k )、计算l o g D k logD_k l o g D k ，然后返回范围内各个k k k值对应的G a p ( K ) Gap(K)G a p (K )值。
def gap(data, refs=None, nrefs=20, ks=range(1, 11)):
shape = data.shape
if refs == None:
tops = data.max(axis=0)
bots = data.min(axis=0)
dists = scipy.matrix(np.diag(tops - bots))
rands = scipy.random.random_sample(size=(shape[0], shape[1], nrefs))
for i in range(nrefs):
rands[:, :, i] = rands[:, :, i] * dists + bots
else:
rands = refs
gaps = np.zeros((len(ks),))
for (i, k) in enumerate(ks):
disp = sum_distance(data, k)
refdisps = np.zeros((rands.shape[2],))
for j in range(rands.shape[2]):
refdisps[j] = sum_distance(rands[:, :, j], k)
gaps[i] = np.lib.scimath.log(np.mean(refdisps)) - np.lib.scimath.log(disp)
return gaps
调用代码
先生成随机点云，此处生成了一组6个聚类中心的3维点云。
之后调用gap函数。
my_data = blobs(centers=6, random_state=121, n_features=3)
gaps = gap(my_data, nrefs=100)
结果输出
此处顺便做了之前两种实现的效率对比，下面两个输出分别对应K-means和Fast K-means实现。
[-0.05127935 -0.01656365 0.30313623 0.78059812 1.56439394 1.70980351
1.67235943 1.63859173 1.62538263 -0.84390957]
best k: 6
58.349995613098145
[-0.05141505 -0.01660111 0.30252902 0.78299894 1.56696863 1.70689603
1.67327937 1.63563849 1.62526949 -0.8487873 ]
best k: 6
4.088269472122192
结果分析
可以看到两种实现的gap值（两个列表输出）有略微不同，但都得到k = 6 k=6 k =6的结果，而之前生成的数据正是6 6 6个中心，说明算法成功检测出最优的k k k值。
此外，可以看到K-means版本的gap函数运行需要58秒，而Fast K-means只需要4秒，速度差距超过十倍，之前的优化还算是效果拔群的。
Gap Statistic总结
- 理论上可以自动化找出最优的K。
- 但由于实现上需要借助特定的K-means算法，而K-means具有局部最优的特点，因此该算法找出的K并不一定是最优的。
- 再加上期望使用蒙特卡洛方法，想得到稳定、靠谱的答案需要花费更多时间进行频率测试。
- 有时会出现多个峰值（设想3类别聚类完全可以对半分成6类），一般根据经验主义选取第一个。
- 无论如何，用于自动化估计较优的K，Gap statistic已经足够了。
[1] Tibshirani R , Hastie W T . Estimating the number of clusters in a data set via the gap statistic[J]. Journal of the Royal Statistical Society B, 2001, 63(2):411-423.
非文献参考
【机器学习】K-means（非常详细）
https://zhuanlan.zhihu.com/p/78798251
A Python implementation of the Gap Statistic
https://gist.github.com/michiexile/5635273
Nuts and Bolts of NumPy Optimization Part 2: Speed Up K-Means Clustering by 70x
https://blog.paperspace.com/speed-up-kmeans-numpy-vectorization-broadcasting-profiling/
Original: https://blog.csdn.net/mkr67n/article/details/125777929
Author: mkr67n
Title: 【机器学习】K-means算法Python实现教程
原创文章受到原创版权保护。转载请注明出处：https://www.johngo689.com/730504/
转载文章受原作者版权保护。转载请注明原作者出处！