视觉机器学习20讲 K-means

from:http://m.blog.youkuaiyun.com/blog/kyu_saku/44495223

[原][几行Matlab代码之]kmeans

提示 
kmeans是什么不在此文中叙述，如有需要请先阅读维基百科或其他文献。本文主要记录litekmeans的代码优化原理。

kmeans计算过程 
随机给定k个初始数据中心点，交替执行下面两个步骤直到收敛： 
1. 将每个数据点分配到离它最近的中心点集合里面； 
2. 更新每个集合的中心点； 
收敛时，相邻两次迭代数据所属的集合不发生改变。

算法执行的结果是，数据点会根据距离被分配到k个集合中。 
这个算法的应用很广泛，在数据分析中，两个距离接近的样本点，通常认为其相似。

正打算用kmeans的时候，我读了一下matlab中的代码，不是优化的代码，于是找到了下面的它。

litekmeans 
两处优化： 
1.在分配数据点时，算法所关心的并非真正的距离值，而是对距离的排序。因而，在计算过程中可以去除那些并不影响排序的计算项。 
观察 ||xi−cj||2=||xi||2+||cj||2−2<xi,cj> 这条距离计算式可以发现，在kmeans迭代过程中， ||xi||2 不影响排序，因而可以不进行计算。 
问题转化为求:  min(||cj||2−2<xi,cj>)  
为了方便矩阵运算，将问题进一步改写为：  max(2<xi,cj>−||cj||2)  (个人理解：其实都是nkd的复杂法，只是为了方便matlab的内部快速运行能调用bsxfun,复制||Cj||而改写的)
其中 ||cj||2      的矩阵大小远小于 <xi,cj> 的矩阵大小，于是改变常数2的位置将进一步优化计算。最后该问题的求解进一步转换为： max(<xi,cj>−||cj||2/2)

2.将中心点的更新矢量化。这里主要引入了示性矩阵E，如果样本i属于集合k，则E(i,k)=1/n(k)，其它为0。之后中心点的计算就由矩阵运算表示（公式难打，请看着代码脑补）。

Code

function [label, centroids] = litekmeans(X, k)
n = size(X, 2);
last = 0;
label = X(:, randsample(n,k)); %随机初始化
while any(label ~= last)
    E = sparse(1:n, label, 1, n, k, n); %构建示性矩阵
    m = X*(E*spdiags(1./sum(E,1)',0,k,k)); %计算每个集合的均值
    last = label; %保存上一次计算结果
    [~,label] = max(bsxfun(@minus,m'*X,dot(m,m,1)'/2,[],1); %分配样本所属集合
end
end

Trick总结 
* 仔细研究计算过程，去掉不必要的计算； 
* 运用矩阵运算性质，将计算矢量化（本文中运用了示性矩阵和对角矩阵）； 
* 在Matlab中稀疏矩阵乘法运算速度快于稠密矩阵，将必要矩阵存储为稀疏矩阵形式。

友情声明 
* 原始的代码优化在多年前就被人提出并且广泛使用，我通过在国外的朋友看到了原作者的原稿。 
* 本文中的代码为最简形式，Deng Cai提供了更完善的版本方便使用，”Litekmeans: the fastest matlab implementation of kmeans,” Available at:http://www.zjucadcg.cn/dengcai/Data/Clustering.html

相关的查阅：

浅入浅出：从Kmeans到Kmeans++
http://www.letiantian.me/2014-03-15-kmeans-kmeans-plus-plus/
matlab函数 bsxfun浅谈
http://blog.sina.com.cn/s/blog_9e67285801010ttn.html
litekmeans
http://www.cad.zju.edu.cn/home/dengcai/Data/code/litekmeans.m

对比matlab自带kmeans，机器学习20讲里kmeans，litekmeans(k = 5; % 聚类个数 max_iter = 100; %最大迭代次数)：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%功能：演示Kmeans聚类算法在计算机视觉中的应用
%实现如何利用Kmeans聚类实现图像的分割；
%环境：Win7，Matlab2012b
%Modi: NUDT-VAP
%时间：2014-10-17
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function kmeans_demo1()
clear;close all;clc;
%% 读取测试图像
im = imread('city.jpg');
imshow(im), title('Imput image');
%% 转换图像的颜色空间得到样本
cform = makecform('srgb2lab');
lab = applycform(im,cform);
ab = double(lab(:,:,2:3));
nrows = size(lab,1); ncols = size(lab,2);
X = reshape(ab,nrows*ncols,2)';
figure, scatter(X(1,:)',X(2,:)',3,'filled');  box on; %显示颜色空间转换后的二维样本空间分布
%print -dpdf 2D1.pdf
%% 对样本空间进行Kmeans聚类
k = 5; % 聚类个数
max_iter = 100; %最大迭代次数
mode = 3;%1:matlab kmeans;2:kmeans in book;3:lite kmeans;
fprintf('data size:%d\n',length(X));
tic;
if mode == 1
    fprintf('Matlab k-means\n');
    [labels,centroids] = litekmeans(X', k, 'MaxIter',max_iter); 
    centroids = centroids';
    labels = labels';
elseif mode == 2
    fprintf('book k-means\n');
    [centroids, labels] = run_kmeans(X, k, max_iter); 
elseif mode == 3
    %%litekmeans
    fprintf('litekmeans\n');
    [labels,centroids] = litekmeans(X', k, 'MaxIter',max_iter); 
    centroids = centroids';
    labels = labels';
end
fprintf('Elapsed time is %d seconds\n',toc);
%% 显示聚类分割结果
figure, scatter(X(1,:)',X(2,:)',3,labels,'filled'); %显示二维样本空间聚类效果
hold on; scatter(centroids(1,:),centroids(2,:),60,'r','filled')
hold on; scatter(centroids(1,:),centroids(2,:),30,'g','filled')
box on; hold off;
%print -dpdf 2D2.pdf

pixel_labels = reshape(labels,nrows,ncols);
rgb_labels = label2rgb(pixel_labels);
figure, imshow(rgb_labels), title('Segmented Image');
saveas(gcf,sprintf('mode_%d',mode), 'png');
%print -dpdf Seg.pdf
end

function [centroids, labels] = run_kmeans(X, k, max_iter)
% 该函数实现Kmeans聚类
% 输入参数：
%                   X为输入样本集，dxN
%                   k为聚类中心个数
%                   max_iter为kemans聚类的最大迭代的次数
% 输出参数：
%                   centroids为聚类中心 dxk
%                   labels为样本的类别标记

%% 采用K-means++算法初始化聚类中心
  centroids = X(:,1+round(rand*(size(X,2)-1)));
  labels = ones(1,size(X,2));
  for i = 2:k
        D = X-centroids(:,labels);
        D = cumsum(sqrt(dot(D,D,1)));
        if D(end) == 0, centroids(:,i:k) = X(:,ones(1,k-i+1)); return; end
        centroids(:,i) = X(:,find(rand < D/D(end),1));
        [~,labels] = max(bsxfun(@minus,2*real(centroids'*X),dot(centroids,centroids,1).'));
  end
  
%% 标准Kmeans算法
  for iter = 1:max_iter
        for i = 1:k, l = labels==i; centroids(:,i) = sum(X(:,l),2)/sum(l); end
        [~,labels] = max(bsxfun(@minus,2*real(centroids'*X),dot(centroids,centroids,1).'),[],1);
  end
  
end

对比结果：

data size:120000
Matlab k-means
Elapsed time is 3.985365e-01 seconds

data size:120000
book k-means
Elapsed time is 1.437269e+00 seconds

data size:120000
litekmeans
Elapsed time is 8.787930e-01 seconds

%1:matlab kmeans;

2:kmeans in book;

3:lite kmeans;