在Matlab实现Kmeans算法（每行代码带注释）

Qun.da

已于 2022-12-06 15:17:53 修改

阅读量8.7k

点赞数 7

文章标签： matlab kmeans 算法

于 2022-12-04 19:56:18 首次发布

本文链接：https://blog.csdn.net/qq_49665818/article/details/128175841

版权

一、前言

本人对机器学习、人工智能算法方面没什么研究，只是学习过程中恰好碰到了。

一开始看Kmeans算法只是为了图像（矩阵）的VQ（vector quantization），找了网上不少资料，跟VQ相关的比较多是LBG算法，单独找kmeans跟VQ联系不起来，后面研究了一下，得到这篇博客主要想表达的内容。

二、VQ概述

VectorQuantization (VQ)是一种基于块编码规则的有损数据压缩方法。事实上，在 JPEG 和 MPEG-4 等多媒体压缩格式里都有 VQ 这一步。它的基本思想是：将若干个标量数据组构成一个矢量，然后在矢量空间给以整体量化，从而压缩了数据而不损失多少信息。

在以前，VQ运用的一个难点在于它要要解决一个多维积分（multi-dimensional integration）的问题。后来，在1980年，Linde, Buzo和Gray（LBG，这个缩写也是LBG算法的命名）提出一种基于训练序列的VQ设计算法，对训练序列的运用绕开了多维积分的求解，使得世上又诞生了一种经典的被世人称为LBG-VQ的算法！它一直延绵至今，经典永不褪色。

三、Kmeans算法

K-means 有一个著名的解释：牧师—村民模型：

有四个牧师去郊区布道，一开始牧师们随意选了几个布道点，并且把这几个布道点的情况公告给了郊区所有的村民，于是每个村民到离自己家最近的布道点去听课。
听课之后，大家觉得距离太远了，于是每个牧师统计了一下自己的课上所有的村民的地址，搬到了所有地址的中心地带，并且在海报上更新了自己的布道点的位置。
牧师每一次移动不可能离所有人都更近，有的人发现A牧师移动以后自己还不如去B牧师处听课更近，于是每个村民又去了离自己最近的布道点……
就这样，牧师每个礼拜更新自己的位置，村民根据自己的情况选择布道点，最终稳定了下来。

K-means 的算法步骤为：

选择初始化的 k 个样本作为初始聚类中心 a=a1,a2,…ak ；
针对数据集中每个样本 xi 计算它到 k 个聚类中心的距离并将其分到距离最小的聚类中心所对应的类中；
针对每个类别 aj ，重新计算它的聚类中心 aj=1|ci|∑x∈cix （即属于该类的所有样本的质心）；
重复上面 2 3 两步操作，直到达到某个中止条件（迭代次数、最小误差变化等）。

四、Matlab代码实现过程

（代码由参考文献中代码修改而来）

function [W, E_in, V] = KMeans(data, K)
%W是k个中心点；E_in是聚合效果，显示所有点的平均距离;V用中心点表示的新的数据
[N, d] = size(data);             %d个n维的点
% init W
sampleIds = randsample(d, K);    %从d个点中随机选择k个点作为中心点
W = data(:,sampleIds);                   %以这k个点为中心形成簇类
labels_u = zeros(1,d);                   %初始换建立一个1行d列的零数组
stop = true;  
while stop                                %把true复制给stop，需要一直循环
    stop = false;
    for i = 1:d                           %从第1个点一直到第d个点,得到每个点与对应最近中心
        x = data(:,i);                   %读取第1个数据放到X里面
        % check label        
        label = 0;                        %初始化label为0，代表是第几个簇类
        dist = 0;                         %初始化dist距离为0
        for j = 1:K                       %计算到达三个中心点的距离，依次推断属于哪个簇类
            tmp_dist = norm(x-W(:,j));        %计算欧式距离
            if label == 0 || tmp_dist < dist       %如果是第一次计算lable=0或者此时的距离小于上一次计算出的距离
                label = j;                         %当前的点暂时属于第j个聚类
                dist = tmp_dist;                   %欧式距离更新为当前的更小值
            end                                    
        end                                        %循环结束
        if labels_u(i) == label                    %如果第个i点不等于label
            stop = stop | false;                          %继续循环
        else
            stop = stop | true;
            labels_u(i) = label;                       %第个i点属于第label个簇类
        end
    end
    if stop == false                                %退出循环
        break;
    end
    %update W                                      %更新中心点
    new_W = zeros(N, K);                           %初始化中心点，并全部清零
    labels_count = zeros(1,K);                    %统计不同簇类的个数
    for i = 1:d                                    %遍历所有点
        label = labels_u(i);                       %提取出簇类标志
        new_W(:,label) = new_W(:,label) + data(:,i);    %相同簇类data数据之和
        labels_count(label) = labels_count(label) + 1;     %属于相同簇类的点的个数
    end
    for i = 1:K   %
        new_W(:,i) = new_W(:,i)/labels_count(i);        %初始化的中心点除以每个聚类里面总的个数，求出新质心
    end
    W = new_W;                         %用新的W来代替
end
E_in = 0;
V = zeros(N,d);
for i = 1:d                            %d个点需要重新遍历
    label = labels_u(i);               %将label标签提取出来
    V(:,i) = W(:,label);                 %用中心表示新数据
    E_in = E_in + norm(data(i)-W(:,label));     %每一个点跟对应中心的距离，所有的距离应该是欧式距离的和
end
E_in = E_in/d;                         %欧式距离的和除以d,每个点距离的均值，表示聚合的效果
end

输入参数：data 为原始数据，对应一列为“一点”，列数即为“点数”；k 表示需要的聚类中心个数

输出参数：W 为得到的聚类中心；E_in 所有点到对应中心距离的平均值，表示聚合效果；V 用对应聚类中心代替原有数据得到的新数据。

由VQ的相关概念可知，V可用 V = WH 进行表示，即用V近似表示原有数据data，或用W表示data达到数据降维的效果。系数H可在上面代码中倒数第五行加入 H(label,i) = 1 ，或得到输出结果后用 H = W\V 得到，表示原来第i个数据的聚类中心为第label个。