K-Means聚类算法（二）：算法实现及其优化

清雨影

2 年前

首先介绍一下题图，这个是个小萝莉的照片（如下），每一个点都可以视作是一个三维向量（点？）（RGB三通道图片），那么，使用K-Means算法对这些点进行聚类，我们就很容易得到几个中心点和几类，把同一类的数据点（像素点）用中心点表示就可以得到压缩后的图片，以上分别是把3通道×每通道8bit=24bit的像素点压缩为1bit，2bit，3bit和4bit的效果，可以看到，虽然信息损失了6倍，但是对画质的影响没有想象的大。

忘了说了，我把图像压缩的代码传上去了，但是图片没办法传送，因为本人的GitHub抽风了（准确的说是网络抽风）

DataScienceNote/KMeans_Img.m at master · TsingJyujing/DataScienceNote · GitHub

请准备好loli.jpg文件，并且配合k_means.m一起食用~

k_means.m代码放在了GitHub上，大家自己下载着玩:P
DataScienceNote/k_means.m at master · TsingJyujing/DataScienceNote · GitHub

其次修正一下上次的说法：并不只是欧氏空间可以用这个方法，只要是对于提供了以下几种算子且满足以下性质的空间都能用：
首先有个计算任意两点距离（或者叫不相似度）的算符记作
且有交换性：

其次有求任意一堆点的平均值的算法：
求出以后使得：

就可以了。
如果这一段有错误，求指出哈~

这一篇依旧不会涉及数学原理，我在纠结中，因为要从头说清楚EM算法是一件很费工夫的事情：首先是含有隐变量(Latent Variables)的模型去估计实际的数据，然后从概率论出发，到极大似然估计确定参数和E（Expectation）步和M（Maximization）步。以上是我清楚的，还有我到现在有云里雾里的收敛性证明。
并不是说费工夫就不写了，逃了，而是用最少的公式和步骤说清楚和K-Means有关的部分。

下一章会比较高能，说是高能，其实也就是求个偏导数，极大似然估计啥的而已……而已……

Part Ⅰ 编程实战的时候那些坑：
实战的时候是会遇到很多坑的，首先说第一个：
一、如何初始化中心点
事实上，初始化的不好，是会震荡的（当然，有了后面的解决震荡的办法要好得多）！
初始化的时候，不妨先估计一下数据的中心在哪里，数据的尺度（叫Scale合适么）有多大，然后在生成数据的时候以中心点为中心，以尺度为因子乘以服从标准正态分布的数据。
大致就是这样的一段代码：

group_center = mean(data_set);
group_range = range(data_set);
centers = (randn(K,dim).*repmat(group_range,K,1)./3+repmat(group_center,K,1));

至于为什么要除以3，是因为这里，一般数据不会越过。
这里根据每一个维度分别做了Scale，使得其能充填在数据中心点附近，然而又不太远。

二、如何计算没有“争取”到任意一点的中心点？
有的时候某些中心点太强势，直接把点抢没了：主要出现在点少类多的情况下。
这个时候，就需要忽略这些中心点。但是中心点不调整也不好，有可能就在那个地方注定孤独一生了。
这个时候我们把这个点初始化到中心点附近的地方：

centers(i,:) = (randn(1,dim).*group_range./3+group_center);

让他重新来过。
这一点在有权重参与计算的时候特别重要么，因为没有争取到点，所以这个时候权重会是0，就出现除零错误了。

三、我怎么知道有没有收敛呢？
这个在书中已经说了，使用代价函数：

其中：

的解释：如果第j个数据点属于第i类，那就记作1，否则记作0的一个 大小的矩阵。代价函数的差分值小于一定数值的时候（N次越不过最小值点）即可认为收敛了。

（以下一段算是题外废话）
这里需要补充另外两种评价聚类算法的方法：F-Measure和信息熵
但是使用这两个代价函数以后就不能用原来提出的迭代算法了，因为原来的EM算法是针对上面的那个代价函数的。
关于这两个代价函数和如何推导会在下一篇讲清楚。

四、代价函数不收敛，怎么破？
其实有一些论文和偏理论的书讲了代价函数的收敛性和初值敏感性的证明，然而不管是本文还是下一篇文章我都不准备讨论这些证明，主要原因是：我！不！太！懂！……（逃

（如果有希望知道这些但是不知道哪里找资料的大犇，我推荐你《统计学习方法》中关于EM算法的收敛性讨论的部分，只不过该书没有说K-Means算法，直说了硬币游戏模型和Gaussian混合模型，真正找EM算法在K-Means上的推导在Pattern Recognition and Machine Learning上，不过没有求偏导和导出求中心点的详细过程，直说了句we can easily solve，翻译过来就是“易证”、“易知”、“易求”或者“我们把这个简单有趣的证明留给读者”之类的鸟话）

首先说一下什么时候容易发生震荡：在数据点个数比较少而且比较稀疏的时候容易发生这种事情，发生的原因大约有两种常见的：
1、陷入某个环里，然后开始震荡，你就会看见某个中心点绕了一圈一圈的……低频振荡。
2、两个点交换来交换去，每次交换不改变J的值就收敛了，如果交换以后不幸影响了其它的点，就出现了高频振荡（多一句嘴，为什么称之为高频呢？因为达到了奈奎斯特采样频率最大值啊）。
这个时候给出一种简单的解决方案：阻尼。
简而言之，就是更新自己位置的时候考虑一下原来的位置，一般阻尼比（在0~1之间取值）决定收敛速度，收敛的慢了也就不容易震荡，也就越容易陷入局部极小值，也就是说，不震荡的情况下我们应该把阻尼比尽可能取小一点

其中：
是最后中心点的取值，是当前集合的中心点，是原来的中心点坐标。

代码中的体现就是：

dumper = 0.1;
...
centers = (1-dumper).*centers + dumper.*lst_cnt;

五、如何给某些点设置权重
首先明确一个问题：为什么要设置权重？
设置权重一般发生在两种情况下：
1，数据集中有大量的重复点，而且数据量比较大计算非常的烧CPU。
2，数据本身有轻重缓急之分，例如只买一两条裤子还扭扭捏捏挑挑选选的我这样的小屌丝客户和乔布斯那种为了自己喜欢的黑色套头衫让停产的生产线再运行的情怀客户的权重肯定是不一样的（请一口气读完）。
可以看看第一篇文章的标题图（封面）就是带权聚类的一个结果：

那么问题来了：如何设置权重？
设置权重其实不困难，首先要明确代价函数是什么：是点距离的和啊！那么更新策略是什么？是找中心点啊！那么现在我们就需要引入加权平均值的算法：

每次移到这个地方就可以了，其中满足范数是1，这是文艺说法。
普通说法就是
如果不满足也懒得归一化可以：

这就是我程序中使用的公式。
具体在代码中的表现就是：

centers(i,:) = sum(data_set(cls_idx,:).*repmat(weight(cls_idx),1,dim))..../sum(weight(cls_idx));

Part Ⅱ Matlab编程技巧：

一、如何愉快的求取范数/距离？
首先，我们一起确认一下范数的公式：
对于任意一个向量：
其范数是：

随后，我们确认一下求范数的几个步骤：
对于已知的一点pnt和行向量组data_set，首先确认行向量组的向量个数：

[N,~]=size(data_set )

随后用repmat把pnt复制N遍，变成和data_set同样大小的一个矩阵。
然后那么一减啊，再求个平方，就可以开始求和了（就是那个Σ），但是是横向拍扁求和，因为我们是行向量啊。
所以给sum函数第二个参数：2
最后开平方即可：

dist_list = sqrt(sum((data_set - repmat(pnt,N,1)).^2,2))

二、如何愉快的计算某一个点属于哪一类：
我们刚才求出了很多的列（dist_list ），这些就是到每个点的距离的列表，这里一共有k个点，所以做成一个k*N的矩阵：

    for i = 1:K%求到每一个中心的距离向量，且不用张量求更加节约空间dist_vec(:,i) = sqrt(sum((data_set - repmat(centers(i,:),N,1)).^2,2));end

这个时候，到第一个点距离最短那么就是1，第二个就是2，如果还用For循环一个个判断你就out了，min函数的第二个返回值是min值所在的位置：

[~,cls_vec] = min(dist_vec,[],2);

只要你乐意，其实还可以生成{r}_{n,k}矩阵（PRML书里用的，偏理论）

rnk = zeros(k,N);
rnk = rnk((0:(N-1)).*k+cls_vec)';

其实没必要。

三、用什么做返回值
这是个仁者见仁丹、智者见智障的问题，Matlab中用刚才计算的cls_vec作为返回值，其实还可以设计一些更加实用的返回值（看你怎么用了）

not_empty_cls = unique(cls_vec);
k_means_result = cell(length(not_empty_cls),1);
for i = 1:length(not_empty_cls)cls = not_empty_cls(i);sub_res = cell(3,1);cls_idx = find(cls_vec==cls);sub_res{1} = centers(cls,:);sub_res{3} = cls_idx;sub_res{2} = data_set(cls_idx,:);k_means_result{i} = sub_res;
end

首先看一下有那几类不空的，好、假设有N类。
然后生成一个的Cell向量，每一个Cell里面放3个Cell：
1，中心点，2，属于这一类的原始数据，3，刚才那些数据在输入数据里的Index。
这样的好处是：信息比较全面，调用方便。
缺点是不够直观么，层次复杂。
你只要开心可以两种都返回，没人拦着你（就好像我代码里示范的那样）