k-means推导和python实现

---

一、算法流程

首先有2点需要注意：

1. k-means和 k-nn(k邻近法) 是不一样的，不要混淆。
2. k-means本身理解起来不难，讲这个是为了引出后面的EM算法，两者有共通之处。

算法流程如下所示，这里以$k=3$, 数据维数$D=2$, 数据样本个数$N=500$为例：

1. 给定K个聚类中心，用 $\mu$ 来表示， 并且适当进行初始化

2. 现在对于给定的$\mu =({\mu }_1,{\mu }_2,{\mu }_3)$, 找出500个样本数据中距离${\mu }_k$最近的一批数据标注为类别$k$

3. 对于属于k类的所有数据进行求平均，将这个平均值作为新的${\mu }_k$使用，同时也得到新的$\mu =({\mu }_1,{\mu }_2,{\mu }_3)$

4. 对比更新前后的$\mu$值，如果其差量收束的话停止更新，否则重复第2步

二、推导

• $x$: $D$维的数据
• $X=\{x_1,x_2,..,x_N\}$: $N$个数据样本
• $K$: 聚类数，已知
• ${\mu }_k(k=1,2,..,K)$: $D$ 维的聚类中心(centriod)
• $r_{nk}$: 第$n$个样本属于$k$类的话该值为1，否则为0

损失函数定义如下：
$$J=\sum _{n=1}^N\sum _{k=1}^K{r}_{nk}||x_n-{\mu }_k|{|}^2$$

该损失函数的最优化流程如下：

1. 固定${\mu }_k$关于$r_{nk}$求偏微分，将其最小化
   令$d_{nk}=||x_n-{\mu }_k|{|}^2$
   $$J=\sum _{n=1}^N(r_{n1}{d}_{n1}+r_{n2}{d}_{n2}+...+r_{nk}{d}_{nk})$$
   为了实现$J$的最小化，等同于将$(r_{n1}{d}_{n1}+r_{n2}{d}_{n2}+...+r_{nk}{d}_{nk})$最小化，因为$r_{nk}$取值为$\{0,1\}$，所以该最小值等同于最小的$d_{nk}$，因此除了最小的$d_{nk}$的$r_{nk}$为1，其他$r_{nk}$均为0
   $$r_{nk}=\{\begin{array}{l}1(k=argmin||x_n-{\mu }_k|{|}^2)\\ 0(otherwise)\end{array}$$

2. 固定$r_{nk}$关于${\mu }_k$求变微分，将其最小化
   $$\begin{gathered} \frac{\partial J}{\partial {\mu }_k}=\frac{\partial }{\partial {\mu }_k}\sum _{n=1}^N\sum _{k=1}^K{r}_{nk}||x_n-{\mu }_k|{|}^2 \\ =\frac{\partial }{\partial {\mu }_k}\sum _{n=1}^N\sum _{k=1}^K{r}_{nk}(-2x_n^T{\mu }_k+{\mu }_k^T{\mu }_k) \\ =\sum _{n=1}^N\sum _{k=1}^K{r}_{nk}(-2\frac{\partial }{\partial {\mu }_k}{x}_n^T{\mu }_k+\frac{\partial }{\partial {\mu }_k}{\mu }_k^T{\mu }_k) \\ =-2\sum _{n=1}^N{r}_{nk}(x_n-{\mu }_k)=0 \end{gathered}$$
   最后再展开一下
   $$\sum _{n=1}^N{r}_{nk}{x}_n=\sum _{n=1}^N{r}_{nk}{\mu }_k$$
   $${\mu }_k^{\ast }=\frac{{\sum }_{n=1}^N{r}_{nk}{x}_n}{{\sum }_{n=1}^N{r}_{nk}}$$

最后这个最优的${\mu }_k^{\ast }$也很好理解，就是$k$类所有数据样本的平均值。

三、python代码实现

具体的代码实现可以参考我的github，这里只贴出一部分核心代码

# k-means算法实现
from collections import Counter
iterations = 100
for iter in range(iterations):
  r = np.zeros(N)
  for i in range(N):
    r[i] = np.argmin([np.linalg.norm(data[i] - mu[k]) for k in range(K)])
  
  if iter % 10 == 0:
    print(mu)
    plt.figure()
    for i in range(N):
      plt.scatter(data[i, 0], data[i, 1], s = 30, c = color_dict[r[i]], alpha = 0.5, marker = "+")
    
    for i in range(K):
      ax = plt.axes()
      ax.arrow(mu[i, 0], mu[i, 1], mu_true[i, 0] - mu[i, 0], mu_true[i, 1] - mu[i, 1], lw = 0.8, head_width = 0.02)
      plt.scatter(mu[i,0], mu[i,1], c = color_dict[i], marker = "o", edgecolors = "k", linewidths=1)
      plt.scatter(mu_true[i, 0], mu_true[i, 1], c = color_dict[i], marker = "o", edgecolors = "k", linewidths=2)
  
    plt.title("iter:{}".format(iter))
    plt.show()
  
  mu_prev = mu.copy()
  # cnt = dict(Counter(r))
  # N_k = [cnt[k] for k in range(K)]
  mu = np.array([np.sum(data[r == k], axis = 0)/len(data[r == k]) for k in range(K)])
  diff = mu - mu_prev

本文参考了这篇日文博客

本文地址：https://blog.csdn.net/dhj_tsukuba/article/details/110392431