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

来源：互联网发布：如何卸载软件 mac 编辑：程序博客网时间：2024/05/01 12:46

以下内容摘自百度百科。

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);endfigure(1);scatter(x_row,x_col,10,'ko');axis([-0.5 1.5 -0.5 1.5]);axis squaregrid 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;endn_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;endif n_iter == MAX_ITER_NUMBER    fprintf('超过最大迭代次数\n');end%% 聚类后画图figure(2);hold on;grid on;axis([-0.5 1.5 -0.5 1.5]);axis squaretitle('聚类后');% 分好类的模式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;    endendscatter(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:19clear;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);endfigure(1);scatter(x_row,x_col,10,'ko');axis([-0.5 1.5 -0.5 1.5]);axis squaregrid 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 offaxis([-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);endfigure(1);scatter(x_row,x_col,10,'ko');axis([-0.5 1.5 -0.5 1.5]);axis squaregrid 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 squaretitle('聚类后');% 分好类的模式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;    endendscatter(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    endend    %% 选取初始点 方法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

效果非常爽啊……

0 0