经典算法（5）：K-均值算法（K-Means）

以下内容摘自百度百科。

K-means算法是硬聚类算法，是典型的基于原型的目标函数聚类方法的代表，

它是数据点到原型的某种距离作为优化的目标函数，

利用函数求极值的方法得到迭代运算的调整规则。

k-means 算法缺点
① 在 K-means 算法中 K 是事先给定的，这个 K 值的选定是非常难以估计的。很多时候，事先并不知道给定的数据集应该分成多少个类别才最合适。这也是 K-means 算法的一个不足。有的算法是通过类的自动合并和分裂，得到较为合理的类型数目 K，例如 ISODATA 算法。关于 K-means 算法中聚类数目K 值的确定在文献中，是根据方差分析理论，应用混合 F统计量来确定最佳分类数，并应用了模糊划分熵来验证最佳分类数的正确性。在文献中，使用了一种结合全协方差矩阵的 RPCL 算法，并逐步删除那些只包含少量训练数据的类。而文献中使用的是一种称为次胜者受罚的竞争学习规则，来自动决定类的适当数目。它的思想是：对每个输入而言，不仅竞争获胜单元的权值被修正以适应输入值，而且对次胜单元采用惩罚的方法使之远离输入值。
② 在 K-means 算法中，首先需要根据初始聚类中心来确定一个初始划分，然后对初始划分进行优化。这个初始聚类中心的选择对聚类结果有较大的影响，一旦初始值选择的不好，可能无法得到有效的聚类结果，这也成为 K-means算法的一个主要问题。对于该问题的解决，许多算法采用遗传算法（GA），例如文献 中采用遗传算法（GA）进行初始化，以内部聚类准则作为评价指标。
③ 从 K-means 算法框架可以看出，该算法需要不断地进行样本分类调整，不断地计算调整后的新的聚类中心，因此当数据量非常大时，算法的时间开销是非常大的。所以需要对算法的时间复杂度进行分析、改进，提高算法应用范围。在文献中从该算法的时间复杂度进行分析考虑，通过一定的相似性准则来去掉聚类中心的侯选集。而在文献中，使用的 K-means 算法是对样本数据进行聚类，无论是初始点的选择还是一次迭代完成时对数据的调整，都是建立在随机选取的样本数据的基础之上，这样可以提高算法的收敛速度。

K-means算法以欧式距离作为相似度测度，它是求对应某一初始聚类中心向量V最优分类，

使得评价指标J最小。算法采用误差平方和准则函数作为聚类准则函数。

n是待聚类的样本个数，

c是类别数（猜，or根据经验，先验知识确定），

u是聚类中心。

先放一个我写的一个最简单的K-means算法。效果明显不是很好……

%% K-means聚类
% 作者：qcy

% 【问题】 现在的问题是：sigma大的时候，反而可以分开，sigma小的时候，分类要分错！

clear;
close all;
clc

%% 产生N个样本，每个样本的特征2维（平面上的点）
N = 1000;
% randn('seed',1);
m_pattern = [];
sigma = 0.2;
% 数据本身就是4类
% 中心
init_m1_center_x = 0;
init_m1_center_y = 0;
init_m2_center_x = 0;
init_m2_center_y = 1;
init_m3_center_x = 1;
init_m3_center_y = 1;
init_m4_center_x = 1;
init_m4_center_y = 0;

% 第1类
for k = 1 : N/4
    m_pattern(k).category = 0; % 初始化为第0类，即还没有分类
    m_pattern(k).feature = [sigma * randn(1,1) + init_m1_center_x; ...
        sigma * randn(1,1) + init_m1_center_y];
end

% 第2类
for k = N/4 +1 : N/2
    m_pattern(k).category = 0; % 初始化为第0类，即还没有分类
    m_pattern(k).feature = [sigma * randn(1,1) + init_m2_center_x; ...
        sigma * randn(1,1) + init_m2_center_y];
end

% 第3类
for k = N/2 +1 : 3*N/4
    m_pattern(k).category = 0; % 初始化为第0类，即还没有分类
    m_pattern(k).feature = [sigma * randn(1,1) + init_m3_center_x; ...
        sigma * randn(1,1) + init_m3_center_y];
end

% 第4类
for k = 3*N/4 +1 : N
    m_pattern(k).category = 0; % 初始化为第0类，即还没有分类
    m_pattern(k).feature = [sigma * randn(1,1) + init_m4_center_x; ...
        sigma * randn(1,1) + init_m4_center_y];
end

% rand('seed',1); % 伪随机打乱
m_pattern = m_pattern(randperm(length(m_pattern)));

% 画图
x_row = zeros(N,1);
x_col = zeros(N,1);
for k = 1 : N
    x_row(k) = m_pattern(k).feature(1);
    x_col(k) = m_pattern(k).feature(2);
end

figure(1);
scatter(x_row,x_col,10,'ko');
axis([-0.5 1.5 -0.5 1.5]);
axis square
grid on;
title('原始数据');

% 保存数据
save m_pattern.mat m_pattern

%% K均值聚类

for k = 1:N
    m_pattern(k).dist2center = inf;
end

% 先验知识：认为有4类
MAX_CNETER_NUMBER = 4; 
MAX_ITER_NUMBER = 1000; % 最大迭代次数

% Begin.
% 1. 随机挑4个样本，作为聚类中心
% 在这里，就用前4个
for k = 1:MAX_CNETER_NUMBER
    m_pattern(k).category = k;
    m_center(k).feature = m_pattern(k).feature;
    m_center(k).new_feature = inf(size(m_pattern(k).feature));
    m_center(k).patternNum = 0; % 这里先设为0，后面统计的时候再累加
    m_center(k).index = k;
end

n_iter = 1;
my_eps = 1e-6; % 聚类中心的前后变化的距离，如果不超过my_eps，就算收敛
isConverge = 0;
center_offset_dist = 1; % 相邻两次迭代，聚类中心移动的距离之和

% 循环条件：
% (1). 当前迭代次数 < 最大允许的迭代次数
% (2). 聚类中心还没有收敛

while (n_iter < MAX_ITER_NUMBER) && (center_offset_dist > my_eps)
    
    % 每个聚类中心所拥有的样本数置为0
    for k = 1:MAX_CNETER_NUMBER
        m_center(k).patternNum = 0;

    end

    % 把每一个点分到距离最近的聚类中心
    for k = 1:N
        feature_temp = m_pattern(k).feature;
        % 与每一个聚类中心进行距离比较
        dist_min = inf;
        idx_dist_min = 0; % 离哪个聚类最近 --> 最小欧式距离判决
        for k_center = 1:MAX_CNETER_NUMBER
            dist = norm(feature_temp - m_center(k_center).feature);
            if dist_min > dist
                dist_min = dist;
                idx_dist_min = k_center;
            end
        end
        
        m_pattern(k).category = idx_dist_min; % 归类
        % 属于 idx_dist_min 这一类的样本总数 + 1
        m_center(idx_dist_min).patternNum = m_center(idx_dist_min).patternNum + 1;
    end
    
    % 计算新的聚类中心
    for k = 1:MAX_CNETER_NUMBER
        feature_temp_sum = [0;0];
        for k_pattern = 1:N
            m_pattern_temp = m_pattern(k_pattern);
            if m_pattern_temp.category == k
                feature_temp_sum = feature_temp_sum + m_pattern_temp.feature;
            end
        end
        m_center(k).new_feature = feature_temp_sum / m_center(k).patternNum;
    end
    
    % 计算聚类中心改变了多少
    center_offset_dist = 0;
    for k = 1:MAX_CNETER_NUMBER
        center_offset_dist = center_offset_dist + ...
            norm(m_center(k).new_feature - m_center(k).feature);
    end
    
    for k = 1:MAX_CNETER_NUMBER  % 上一次的feature需要改一下
        m_center(k).feature = m_center(k).new_feature;
        
    end
      
    n_iter = n_iter + 1;
end

if n_iter == MAX_ITER_NUMBER
    fprintf('超过最大迭代次数\n');
end


%% 聚类后画图
figure(2);hold on;
grid on;
axis([-0.5 1.5 -0.5 1.5]);
axis square
title('聚类后');

% 分好类的模式
m_modes_count = [1 1 1 1];
for k = 1 : N
    m_pattern_temp = m_pattern(k);
    category = m_pattern_temp.category;
    feature = m_pattern_temp.feature;
    
    if category == 1
        x1_row(m_modes_count(category)) = m_pattern_temp.feature(1);
        x1_col(m_modes_count(category)) = m_pattern_temp.feature(2);
        m_modes_count(category) = m_modes_count(category) + 1;
    elseif category == 2
        x2_row(m_modes_count(category)) = m_pattern_temp.feature(1);
        x2_col(m_modes_count(category)) = m_pattern_temp.feature(2);
        m_modes_count(category) = m_modes_count(category) + 1;
    elseif category ==3
        x3_row(m_modes_count(category)) = m_pattern_temp.feature(1);
        x3_col(m_modes_count(category)) = m_pattern_temp.feature(2);
        m_modes_count(category) = m_modes_count(category) + 1;
    else
        x4_row(m_modes_count(category)) = m_pattern_temp.feature(1);
        x4_col(m_modes_count(category)) = m_pattern_temp.feature(2);
        m_modes_count(category) = m_modes_count(category) + 1;
    end
end

scatter(x1_row,x1_col,10,'ro');
scatter(x2_row,x2_col,10,'go');
scatter(x3_row,x3_col,10,'bo');
scatter(x4_row,x4_col,10,'mo');

我只算了一次，聚类中心的初值是用的最开始的几个点，有可能根本就不收敛了，也没有多次尝试去计算。

方差小了，反而还出问题了！-_-！

据说是因为初值没有取好…哎…………

所以似乎应该要试很多个初值…

Matlab有自带的kmeans函数，效果非常好…

%% K-means聚类
% 调用matlab自带的kmeans函数
% qcy
% 时间：2016年6月16日21:05:19

clear;
close all;
clc

%% 产生N个样本，每个样本的特征2维（平面上的点）
N = 1000;
% randn('seed',1);
m_pattern = [];
sigma = 0.2;
% 数据本身就是4类
% 中心
init_m1_center_x = 0;
init_m1_center_y = 0;
init_m2_center_x = 0;
init_m2_center_y = 1;
init_m3_center_x = 1;
init_m3_center_y = 1;
init_m4_center_x = 1;
init_m4_center_y = 0;

% 第1类
for k = 1 : N/4
    m_pattern(k).category = 0; % 初始化为第0类，即还没有分类
    m_pattern(k).feature = [sigma * randn(1,1) + init_m1_center_x; ...
        sigma * randn(1,1) + init_m1_center_y];
end

% 第2类
for k = N/4 +1 : N/2
    m_pattern(k).category = 0; % 初始化为第0类，即还没有分类
    m_pattern(k).feature = [sigma * randn(1,1) + init_m2_center_x; ...
        sigma * randn(1,1) + init_m2_center_y];
end

% 第3类
for k = N/2 +1 : 3*N/4
    m_pattern(k).category = 0; % 初始化为第0类，即还没有分类
    m_pattern(k).feature = [sigma * randn(1,1) + init_m3_center_x; ...
        sigma * randn(1,1) + init_m3_center_y];
end

% 第4类
for k = 3*N/4 +1 : N
    m_pattern(k).category = 0; % 初始化为第0类，即还没有分类
    m_pattern(k).feature = [sigma * randn(1,1) + init_m4_center_x; ...
        sigma * randn(1,1) + init_m4_center_y];
end

% rand('seed',1); % 伪随机打乱
m_pattern = m_pattern(randperm(length(m_pattern)));

% 画图
x_row = zeros(N,1);
x_col = zeros(N,1);
for k = 1 : N
    x_row(k) = m_pattern(k).feature(1);
    x_col(k) = m_pattern(k).feature(2);
end

figure(1);
scatter(x_row,x_col,10,'ko');
axis([-0.5 1.5 -0.5 1.5]);
axis square
grid on;
title('原始数据');

% 保存数据
save m_pattern.mat m_pattern

%% K均值聚类

MAX_CNETER_NUMBER = 4;  
opts = statset('Display','final');
for k = 1:N
    X(k,:) = m_pattern(k).feature;
end
[idx,Center] = kmeans(X,MAX_CNETER_NUMBER,'Distance','cityblock',...
    'Replicates',5,'Options',opts);

figure;
plot(X(idx==1,1),X(idx==1,2),'r.','MarkerSize',12)
hold on;grid on;
plot(X(idx==2,1),X(idx==2,2),'b.','MarkerSize',12)
plot(X(idx==3,1),X(idx==3,2),'g.','MarkerSize',12)
plot(X(idx==4,1),X(idx==4,2),'m.','MarkerSize',12)
plot(Center(:,1),Center(:,2),'kx',...
     'MarkerSize',15,'LineWidth',3)
legend('Cluster 1','Cluster 2','Cluster 3','Cluster 4','Centroids',...
       'Location','NW')
title 'Cluster Assignments and Centroids'
hold off
axis([-0.5 1.5 -0.5 1.5]);
axis square

Matlab自带的kmeans，运行速度极快，估计用了并行？或者直接用C写的也有非常有可能…

另外，网上还有一段，写得很好很好很好！

1. 用空间换了时间；（e.g. 最小距离判别的时候，直接把一个点repmat，和四个中心相减…不像我还用了个for去一个一个比……别人一次性比完了！）

2. 用空间换了简单的思维；

3. 一开始多次随机地选择聚类中心，外层有个大循环，并把每次循环找出来的4个聚类中心都保存起来，最后去评价“哪个最好”；

4. 评价哪次循环找到的几个聚类中心最好，依据的是

a. 聚类后的每一类的所有样本点到该类的聚类中心的距离之和；

b. 把每一类算出来的距离之和再加起来，求个总误差

找到总误差最小的那次循环所找到的聚类中心，就是这N次大循环中，效果最好的一次聚类中心……

主函数。

clear;
close all;
clc


%% 产生N个样本，每个样本的特征2维（平面上的点）
N = 1000;
% randn('seed',1);
m_pattern = [];
sigma = 0.1;
% 数据本身就是4类
% 中心
init_m1_center_x = 0;
init_m1_center_y = 0;
init_m2_center_x = 0;
init_m2_center_y = 1;
init_m3_center_x = 1;
init_m3_center_y = 1;
init_m4_center_x = 1;
init_m4_center_y = 0;

% 第1类
for k = 1 : N/4
    m_pattern(k).category = 0; % 初始化为第0类，即还没有分类
    m_pattern(k).feature = [sigma * randn(1,1) + init_m1_center_x; ...
        sigma * randn(1,1) + init_m1_center_y];
end

% 第2类
for k = N/4 +1 : N/2
    m_pattern(k).category = 0; % 初始化为第0类，即还没有分类
    m_pattern(k).feature = [sigma * randn(1,1) + init_m2_center_x; ...
        sigma * randn(1,1) + init_m2_center_y];
end

% 第3类
for k = N/2 +1 : 3*N/4
    m_pattern(k).category = 0; % 初始化为第0类，即还没有分类
    m_pattern(k).feature = [sigma * randn(1,1) + init_m3_center_x; ...
        sigma * randn(1,1) + init_m3_center_y];
end

% 第4类
for k = 3*N/4 +1 : N
    m_pattern(k).category = 0; % 初始化为第0类，即还没有分类
    m_pattern(k).feature = [sigma * randn(1,1) + init_m4_center_x; ...
        sigma * randn(1,1) + init_m4_center_y];
end

% rand('seed',1); % 伪随机打乱
m_pattern = m_pattern(randperm(length(m_pattern)));

% 画图
x_row = zeros(N,1);
x_col = zeros(N,1);
for k = 1 : N
    x_row(k) = m_pattern(k).feature(1);
    x_col(k) = m_pattern(k).feature(2);
end

figure(1);
scatter(x_row,x_col,10,'ko');
axis([-0.5 1.5 -0.5 1.5]);
axis square
grid on;
title('原始数据');

% 保存数据
save m_pattern.mat m_pattern

%% K均值聚类

MAX_CNETER_NUMBER = 4;  
opts = statset('Display','final');
for k = 1:N
    X(:,k) = m_pattern(k).feature;
end
 [best_Label best_Center best_ind label] = KM(X,MAX_CNETER_NUMBER,'KMeans');

 for k = 1:N
     m_pattern(k).category = best_Label(k);
 end
 
 
%% 聚类后画图
figure(2);hold on;
grid on;
axis([-0.5 1.5 -0.5 1.5]);
axis square
title('聚类后');

% 分好类的模式
m_modes_count = [1 1 1 1];
for k = 1 : N
    m_pattern_temp = m_pattern(k);
    category = m_pattern_temp.category;
    feature = m_pattern_temp.feature;
    
    if category == 1
        x1_row(m_modes_count(category)) = m_pattern_temp.feature(1);
        x1_col(m_modes_count(category)) = m_pattern_temp.feature(2);
        m_modes_count(category) = m_modes_count(category) + 1;
    elseif category == 2
        x2_row(m_modes_count(category)) = m_pattern_temp.feature(1);
        x2_col(m_modes_count(category)) = m_pattern_temp.feature(2);
        m_modes_count(category) = m_modes_count(category) + 1;
    elseif category ==3
        x3_row(m_modes_count(category)) = m_pattern_temp.feature(1);
        x3_col(m_modes_count(category)) = m_pattern_temp.feature(2);
        m_modes_count(category) = m_modes_count(category) + 1;
    else
        x4_row(m_modes_count(category)) = m_pattern_temp.feature(1);
        x4_col(m_modes_count(category)) = m_pattern_temp.feature(2);
        m_modes_count(category) = m_modes_count(category) + 1;
    end
end

scatter(x1_row,x1_col,10,'ro');
scatter(x2_row,x2_col,10,'go');
scatter(x3_row,x3_col,10,'bo');
scatter(x4_row,x4_col,10,'mo');

KM函数。

function [best_Label best_Center best_ind label] = KM(P,K,method)
%%%%-----------------------------------------------------------------------
%   Version 1.0 
%   Author: feitengli@foxmail.com   from DUT
%   CreateTime: 2012-11-29
%%%------------------------------------------------------------------------
%KM   K-Means Clustering or K-Medoids Clustering
%    P is an d-by-N data matrix 
%    K is the clustering number
%    method = KMeans    :K-Means Clustering
%           = KMedoids  :K-Medoids Clustering
%References：
%        1.The Elements of Statistical Learning 2nd Chapter14.3.6&&14.3.10
%%%%-----------------------------------------------------------------------

[d N] = size(P); 
%% 本算法要求数据矩阵P的每列代表一个数据点，如果不是 需要转置矩阵
if d > N
    ButtonName = questdlg('数据维数小于点的个数,是否转置矩阵',...
                          'MATLAB quest','Yes','No','Yes');
    if  strcmp(ButtonName, 'Yes')
        P = P';
        [d N] = size(P); 
%     else
%         return
    end
end
    
%% 选取初始点 方法2 
max_Initial = max(20,N/(5*K));
label = zeros(max_Initial,N);
center = zeros(d,K,max_Initial);
C = zeros(1,N);
%% 主循环
for initial_Case = 1:max_Initial
    
    pointK = Initial_center(P,K);    
    iter = 0;
%     max_iter = 1;
    max_iter = 1e+3;
    % xK = pointK;
    disp(['------------KM进行第 ' num2str(initial_Case) ' 次重新选择初始中心-----------'])
    %% 每次初始化K个中心点后，进行的循环
    while iter < max_iter
        iter = iter+1;
        if mod(iter,50)==0
            disp(['  内部循环进行第 ' num2str(iter) ' 次迭代'])
        end
        %%%根据数据矩阵P中每个点到中心点的距离(最小)确定所属分类
        for i = 1:N
            dert = repmat(P(:,i),1,K)-pointK;
            distK = sqrt(diag(dert'*dert));
            [~,j] = min(distK);
            C(i) = j;
        end
        %%%重新计算K个中心点  
        xK_ = zeros(d,K);
        for i = 1:K
            Pi = P(:,find(C==i));
            Nk = size(Pi,2);
            % K-Means K-Medoids唯一不同的地方：选择中心点的方式
            switch lower(method)
                case 'kmeans'  
                    xK_(:,i) = sum(Pi,2)/Nk;
                case 'kmedoids'
                    Dx2 = zeros(1,Nk);
                    for t=1:Nk
                       dx = Pi - Pi(:,t)*ones(1,Nk);
                       Dx2(t) = sum(sqrt(sum(dx.*dx,1)),2);
                    end
                    [~,min_ind] = min(Dx2);
                    xK_(:,i) = Pi(:,min_ind);
                otherwise
                    errordlg('请输入正确的方法：kmeans-OR-kmedoids','MATLAB error');
            end
        end
        
        % 判断是否达到结束条件
        if xK_==pointK   % & iter>50  %【qcy批注】 我认为，当二者之差小于某个Eps时，认为收敛
            disp(['###迭代 ' num2str(iter) ' 次得到收敛的解'])
            label(initial_Case,:) = C;
            center(:,:,initial_Case) = xK_;
          % plot_Graph(C);
            break
        end
        
        pointK = xK_;
        %xK = xK_;
    end
    if iter == max_iter
         disp('###达到内部最大迭代次数1000,未得到收敛的解') 
         label(initial_Case,:) = C;
         center(:,:,initial_Case) = xK_;
        % plot_Graph(C);
         % break
    end
    
end

%%%%增加对聚类结果最优性的比较  
%距离差
 dist_N = zeros(max_Initial,K);
 for initial_Case=1:max_Initial     
     for k=1:K
         tem = find(label(initial_Case,:)==k);
         dx = P(:,tem)-center(:,k,initial_Case)*ones(1,size(tem,2));
         dxk = sqrt(sum(dx.*dx,1));
         dist_N(initial_Case,k) = sum(dxk);       
         % dist_N(initial_Case,k) = dxk;   
     end     
 end
 
 %%%%对于max_Initial次初始化中心点得到的分类错误
 %%%%取错误最小的情况的Label作为最终分类
 dist_N_sum = sum(dist_N,2); %求K类总的误差
 [~,best_ind] = min(dist_N_sum);
 best_Label = label(best_ind,:); 
 best_Center = center(:,:,best_ind);

Initial_center初始化聚类中心的函数

function center = Initial_center(X,K)  
        N = size(X,2);
        rnd_Idx = randperm(N);  
        center = X(:,rnd_Idx(1:K));  
end   

效果非常爽啊……

【补充一些东西】

1. K-means聚类，除了一开始要指定多少类，其他的“几乎都是盲的”，所谓的非监督嘛……

2. 实验室的卢老板，用K-means解决了方型QAM调制时，星座点歪了，怎么判决的问题，还发表了一篇PJ的文章。祝贺祝贺。

求卢老板以后带飞。