算法笔记：Kmeans聚类算法简介

算法笔记：Kmeans聚类算法简介
- 1. Kmeans算法简介
- 2. Kmeans算法细节
- 3. Kmeans算法收敛性证明
- 4. Kmeans算法的变体
  - 1. cosine距离变体
  - 2. 点积距离版本
- 5. Kmeans算法实现
  - 1. 基于sklearn的kmeans算法
  - 2. python自实现
- 6. 参考链接

1. Kmeans算法简介

Kmeans算是非常经典的一个聚类算法了，早已经被写到教科书里面了，不过很不幸的是，最近干活遇到了这个，然后我发现我已经忘得差不多一干二净了……

所以这里就过来挖个坟，考个古，把这玩意拉出来复习一下。

如前所述，Kmeans算法是一个聚类算法，具体来说，我们输入一个包含NNN个点的点集，我们的目的是要将这NNN个点分为KKK个簇，使得每个点到各自的簇的中心距离之和最小。

用公式来表达的话就是：

s=∑i=1Nminj∈{1,...,K}(d(xi,uj))s = \sum_{i=1}^{N} \mathop{min}\limits_{j \in \{1, ..., K\}}(d(x_i, u_j)) s=i=1∑Nj∈{1,...,K}min(d(xi,uj))

要找到一组uju_juj使得sss最大。

其中，d(x,y)d(x, y)d(x,y)表示x,yx,yx,y两点间的距离，一般我们在这里使用欧氏距离。

2. Kmeans算法细节

Kmeans算法的核心思路是迭代。

首先，我们随机从NNN个点当中选出KKK个点作为簇的中心点。

然后，根据全部的NNN个点到这KKK个中心点之间的距离，我们就可以将这全部的NNN个点进行分类，分配到这KKK个簇当中。

而后，我们更新这KKK个簇的中心，具体来说，我们取这KKK个点的均值点作为这KKK个簇的新的中心。

我们不断地重复上述两个步骤，直到达到迭代上限或者簇的中心点不再发生变化即可。

具体的，我们可以给出上述Kmeans算法的算法整理如下：

step 1: 从NNN个给定点当中随机KKK个点作为KKK个簇的中心点；
step 2: 计算每一个点到这KKK个簇的中心点之间的欧式距离，将其分配到最小的那个簇当中，从而对所有的点进行聚类；
step 3: 对于2中得到的每一个簇，更新其中心点为所有点的均值，即u=∑ixin\bold{u} = \frac{\sum_i \bold{x}_i}{n}u=n∑ixi；
step 4: 重复上述2-3两步，直到迭代次数达到上限或者簇的中心不再发生变化。

而Kmeans的算法的优缺点因此也就比较明显：

优点
- 易实现，易debug
缺点
- 迭代非常耗时，对于大数据量尤其明显；
- 较依赖于初始化中心的选择，不同初始化中心点的选择会带来较大的结果差异；

3. Kmeans算法收敛性证明

现在，给出了kmeans聚类算法之后，我们来考察一下kmeans算法的收敛性，也就是说，为什么kmeans算法的迭代是有效的。

我们使用原始的kmeans算法进行说明，即是说，使用欧式距离来对两点间的距离进行描述，此时，前述提到的loss函数就可以表达为：

s=∑i=1Nminj∈{1,...,K}∣∣xi,uj∣∣s = \sum_{i=1}^{N} \mathop{min}\limits_{j \in \{1, ..., K\}} ||x_i, u_j|| s=i=1∑Nj∈{1,...,K}min∣∣xi,uj∣∣

具体到第kkk次迭代上，即有：

sk=∑i=1Nminj∣∣xi,ujk∣∣s^{k} = \sum_{i=1}^{N} \mathop{min}\limits_{j} ||x_i, u_j^k|| sk=i=1∑Njmin∣∣xi,ujk∣∣

显然，sks^{k}sk是一个大于0的数列，因此，我们只需要证明sks^{k}sk递减，那么数列sks^{k}sk必然收敛。

因此，我们只需要证明sk+1≤sks^{k+1} \leq s^{k}sk+1≤sk即可。

我们考察第kkk次迭代，它分为两步：

对于上一次分类完成的簇，更新簇的中心从uku^{k}uk到uk+1u^{k+1}uk+1；
sk+1′=∑i=1N∣∣xi,ujk+1∣∣s^{{k+1}'} = \sum_{i=1}^{N} ||x_i, u_j^{k+1}|| sk+1′=i=1∑N∣∣xi,ujk+1∣∣
使用新的簇中心uk+1u^{k+1}uk+1对所有的点进行更新；
sk+1=∑i=1Nminj∣∣xi,ujk+1∣∣s^{k+1} = \sum_{i=1}^{N} \mathop{min}\limits_{j} ||x_i, u_j^{k+1}|| sk+1=i=1∑Njmin∣∣xi,ujk+1∣∣

其中，对于步骤二，显然有sk+1≤sk+1′s^{k+1} \leq s^{{k+1}'}sk+1≤sk+1′。因此，我们只要说明步骤一当中的聚类中心变换之后获得的新的sk+1′s^{{k+1}'}sk+1′小于等于sks^{k}sk即可。

而在这步骤一当中，由于簇的成员都没有发生改变，因此，我们要证明的问题也就是：

对一系列点x1,...,xni\bold{x}_1, ..., \bold{x}_{n_i}x1,...,xni，s=∑j=1ni∣∣xj−μ∣∣s = \sum\limits_{j=1}^{n_i} ||\bold{x}_j - \bold{\mu}||s=j=1∑ni∣∣xj−μ∣∣在μ=1ni∑j=1nixj\bold{\mu} = \frac{1}{n_i}\sum\limits_{j=1}^{n_i} \bold{x}_jμ=ni1j=1∑nixj时取到最小值。

而这个问题的解答也是比较简单的，我们求一下s2s^2s2对于μ\bold{\mu}μ的导数在值为000时的μ\bold{\mu}μ即可，而这个证明是简单的，这里就不做展开了。

4. Kmeans算法的变体

这里，我们对kmeans算法进行一点延申。

如前所述，kmeans本质上是一个迭代算法，他的算法有效性本质上来源于kmeans的迭代的收敛性。而进一步的，从上述kmeans算法在欧式距离下的收敛性证明当中我们看到，迭代算法的收敛性其实本质上由来源于对每一个簇更新簇的中心时整体loss的递减关系。

因此，我们可以给出一个更为一般的结论：

如果对于某一个距离度量，我们可以找到簇当中一个恒定的极小值点表达式，且这个表达式仅与簇中的点坐标有关，那么我们就可以用这个极小值点来进行簇中心的迭代，此时kmeans算法总是收敛的。

基于此，我们可以给出kmeans算法的一些变体：

cosine距离版本
点积距离版本

1. cosine距离变体

首先，我们给出cosine距离下的kmeans变体：

step 1: 从NNN个给定点当中随机KKK个点作为KKK个簇的中心点；
step 2: 计算每一个点到这KKK个簇的中心点之间的cosine距离，将其分配到cosine距离最大的那个簇当中，从而对所有的点进行聚类；
step 3: 对于2中得到的每一个簇，更新其中心点为u=1n∑ixi∣∣xi∣∣\bold{u} = \frac{1}{n} \sum\limits_{i} \frac{\bold{x}_i}{||\bold{x}_i||}u=n1i∑∣∣xi∣∣xi；
step 4: 重复上述2-3两步，直到迭代次数达到上限或者簇的中心不再发生变化。

可以看到，这里事实上主要也就是第三步当中的簇中心更新与之前的欧式距离版本的原版kmeans有所区别，其他基本上是完全一模一样的。

因此，我们下面只需要证明在cosine距离下，1n∑icos(xi,u)\frac{1}{n}\sum\limits_{i}cos(\bold{x}_i, \bold{u})n1i∑cos(xi,u)总在u=1n∑ixi∣∣xi∣∣\bold{u} = \frac{1}{n} \sum\limits_{i} \frac{\bold{x}_i}{||\bold{x}_i||}u=n1i∑∣∣xi∣∣xi时取到极大值即可。

这个证明事实上也是相对简单的：

f(μ)=∑i=1ncos(xi,μi)=∑i=1nxi∣∣xi∣∣⋅μ∣∣μ∣∣=(∑i=1nxi∣∣xi∣∣)⋅μ∣∣μ∣∣\begin{aligned} f(\bold{\mu}) &= \sum\limits_{i=1}^{n} \mathop{cos}(\bold{x}_i, \bold{\mu}_i) \\ &= \sum\limits_{i=1}^{n} \frac{\bold{x}_i}{||\bold{x}_i||} \cdot \frac{\bold{\mu}}{||\bold{\mu}||} \\ &= (\sum\limits_{i=1}^{n} \frac{\bold{x}_i}{||\bold{x}_i||}) \cdot \frac{\bold{\mu}}{||\bold{\mu}||} \end{aligned} f(μ)=i=1∑ncos(xi,μi)=i=1∑n∣∣xi∣∣xi⋅∣∣μ∣∣μ=(i=1∑n∣∣xi∣∣xi)⋅∣∣μ∣∣μ

可以看到，显然当 μ∣∣μ∣∣∼(∑i=1nxi∣∣xi∣∣)\frac{\bold{\mu}}{||\bold{\mu}||} \sim (\sum\limits_{i=1}^{n} \frac{\bold{x}_i}{||\bold{x}_i||})∣∣μ∣∣μ∼(i=1∑n∣∣xi∣∣xi)时，即两者同向时，f(μ)f(\bold{\mu})f(μ)可以取到极大值。

因此，我们只要更新μ\bold{\mu}μ为如下表达式即可：

μ=1n(∑i=1nxi∣∣xi∣∣)\bold{\mu} = \frac{1}{n} (\sum\limits_{i=1}^{n} \frac{\bold{x}_i}{||\bold{x}_i||}) μ=n1(i=1∑n∣∣xi∣∣xi)

2. 点积距离版本

同样的，我们可以给出点积版本的kmeans迭代的变体：

step 1: 从NNN个给定点当中随机KKK个点归一化后作为KKK个簇的中心点；
step 2: 计算每一个点到这KKK个簇的中心点之间的点积距离，将其分配到点积距离最大的那个簇当中，从而对所有的点进行聚类；
step 3: 对于2中得到的每一个簇，更新其中心点为u=∑ixi/∣∣∑ixi∣∣\bold{u} = \sum\limits_{i} \bold{x}_i/||\sum\limits_{i} \bold{x}_i||u=i∑xi/∣∣i∑xi∣∣；
step 4: 重复上述2-3两步，直到迭代次数达到上限或者簇的中心不再发生变化。

而要说明上述迭代的收敛性，我们同样只需要说明点积距离下u=∑ixi/∣∣∑ixi∣∣\bold{u} = \sum\limits_{i} \bold{x}_i/||\sum\limits_{i} \bold{x}_i||u=i∑xi/∣∣i∑xi∣∣总可以令簇内部点的点积之和最大即可。

仿上，我们可以快速地写出簇内部所有的点到某个向量之间的点积之和：

f(μ)=∑i=1nxi⋅μ=(∑i=1nxi)⋅μ\begin{aligned} f(\bold{\mu}) &= \sum\limits_{i=1}^{n} \bold{x}_i \cdot \bold{\mu} \\ &= (\sum\limits_{i=1}^{n} \bold{x}_i) \cdot \bold{\mu} \end{aligned} f(μ)=i=1∑nxi⋅μ=(i=1∑nxi)⋅μ

显然，μ\bold{\mu}μ的模长越长，其与向量∑i=1nxi\sum\limits_{i=1}^{n} \bold{x}_ii=1∑nxi的方向越一致，整体的f(μ)f(\bold{\mu})f(μ)就会越大。

而方向方面，我们是很方面就可以控制的，但是模长方面我们却不太可控，因此这里对此进行了一定的限制，将μ\bold{\mu}μ限定为一个单位向量，此时，我们就能找到一个μ\bold{\mu}μ使得f(μ)f(\bold{\mu})f(μ)取到最大值了。

因此，此时的迭代关系即为：

μ=∑inxi∣∣∑inxi∣∣\bold{\mu} = \frac{\sum\limits_{i}^n \bold{x}_i}{||\sum\limits_{i}^n \bold{x}_i||} μ=∣∣i∑nxi∣∣i∑nxi

5. Kmeans算法实现

最后，我们来看一下Kmeans算法的实现。

这里，我们首先基于sklearn库给出一个简易的kmeans实现，一般情况下这也就够用了。

然后我们手撸一个标准版本的基于欧氏距离的kmeans代码，用于加深理解以及方便于后续的各类定制版本的kmeans算法实现。

我们统一定义输入数据格式如下：

import numpy as npn = 100000
dim = 32data = np.random.random(size=(n, dim)) # [n, dim]

1. 基于sklearn的kmeans算法

下面，我们首先来看一下基于sklearn的kmeans算法实现。

这个其实过于简单几乎是无脑调用就行：

from sklearn.cluster import KMeanscluster_num = 100
model = KMeans(n_clusters=cluster_num, max_iter=300)
model.fit(data)

然后，如果我们需要保存下这个kmeans模型的聚类中心，只需要做如下操作即可：

model_path = "model/kmeans_centers.npy"
cluster_centers = model.cluster_centers_
np.save(model_path, cluster_centers)

2. python自实现

下面，我们来手撸一个经典的kmeans算法。

import numpy as np
from collections import defaultdictclass Kmeans:def __init__(self, n_cluster):self.n_cluster = n_clusterdef fit(self, data, max_iter=300):self.cluster_centers_ = self.init_cluster(data)for _ in range(max_iter):# clusteringclusters = defaultdict(list)distances = self.distance_fn(data, self.cluster_centers_)cluster_ids = np.argmin(distances, axis=-1)for i, cid in enumerate(cluster_ids):clusters[cid].append(i)# update cluster centernew_cluster_centers = self.update_cluster(clusters, data)max_shift = np.max(np.abs(new_cluster_centers - self.cluster_centers_))self.cluster_centers_ = new_cluster_centersif max_shift < 1e-3:breakreturndef init_cluster(self, data):n = data.shape[0]ids = np.random.choice(n, self.n_cluster)return np.take(data, ids, axis=0)def update_cluster(self, clusters, data):new_cluster_centers = []for idx in range(self.n_cluster):cluster_data = np.take(data, clusters[idx], axis=0)cluster_center = np.mean(cluster_data, axis=0)new_cluster_centers.append(cluster_center)return np.array(new_cluster_centers)def distance_fn(self, x, y):x = np.expand_dims(x, axis=1) # [i, 1, d]y = np.expand_dims(y, axis=0) # [1, j, d]d = np.sqrt(np.sum((x - y)**2, axis=-1)) # [i, j]return d

6. 参考链接

What is K-Means algorithm and how it works