机器学习练习之k均值

k-means属于聚类分析的其中一种算法，聚类分析在机器学习、数据挖掘、模式识别、决策支持和图像分割中有广泛的应用。聚类是无监督的分类方法，所谓无监督就是没有给定训练数据的标签信息，所以聚类出来的结果的类别是未定义的，而分类的目标是把数据分到已知的类别中。聚类是在给定的数据集合中寻找数据子集合，每个子集合形成一个类簇，簇内间的相似性高，而簇间的相似性低。通俗地说，k-means算法就是“物以类聚人以群分”，把相似性较大的事物聚在一起，相似性较低的事物分开来，从而得到一个个群组。比如在市场营销中，我们通过对市场上的顾客信息进行k-means划分成不同的顾客群体，因而可以分析不同群体性质从而制定不同的营销策略。

k-means的k是指聚类结果把数据划分成k个簇，mean表示簇的质心或中心是簇内数据的均值，当数据无法求得均值时，k-means算法就不适用。

算法过程如下：

算法首先初始化k个簇的中心点或质心。根据样本集中的每一样本与k个簇中心的距离，把它赋给最近的簇，然后重新计算簇中心，重复这个过程直到收敛。k个质心可以从样本集中随机挑选k个样本来代表。

下图表示k-means的收敛过程：

如果k个簇中心的初始值选取得不好，k-means会达到局部最优而不是全局最优，局部最优时的结果可能跟预料的结果大相径庭，所以k-means对初始值选取敏感。k-means还对“噪声”和异常值敏感，因为这些点对均值的计算结果影响很大。

k-means算法是一定能够收敛的，我们可以定义如下的损失函数：

函数J衡量的是每个训练样本到它被划分到的簇的中心点的平方距离的和，J不是关于c和u的凸函数，但是函数J可以通过坐标下降的方法来求出局部最小值，通过固定u来更新c即固定了中心点，把每个样本赋给离它最近的簇中心来减少J，然后通过固定c来更新u，因为使用的是平方距离损失函数，所以求出簇的均值就可以使得J最小，对于这两个步骤，J都是减少的，因此J得到收敛，但由于J不是凸函数，不一定收敛到全局最优解。

k-means算法也是EM思想的体现，EM是最大期望算法，也是机器学习和数据挖掘常用的一种方法。E步可对应把每个样本点划分到离它最近的簇中，M步可对应到求出簇的均值作为簇的中心点。但是，k-means的E步是硬的指定，把样本点划分到最近的簇中，而不像混合高斯模型中的E步是软的指定，样本点可以以概率分到不同的高斯分布中。

该练习是要在图片中把像素分为16个颜色，k-means的算法过程如下：

Matlab代码如下：

%% Exercise: k-means
% k-means算法对鸟图像素聚类
% 先在小鸟图上k-means得出k个中心点颜色，然后在大鸟图上把各个像素的颜色赋值为最近中心点的颜色

% 载入小鸟图像
% A是一个三维数组，A(y,x,1)、A(y,x,2)、A(y,x,3)分别是位置(i,j)的r、g、b像素
A = double(imread('bird_small.tiff'));

dim = size(A,1); %图片宽度或高度
k = 16; %显示的颜色种类

%随机选择k个初始中心点
means = zeros(k,3); %mean存储中心点(r,g,b)颜色
x_rand = ceil(dim*rand(k,1));
y_rand = ceil(dim*rand(k,1));
for i=1:k
    means(i,:) = A(y_rand(i), x_rand(i), :);
end

%k-means算法过程
MAX_ITR = 100;
for itr=1:MAX_ITR
    %这一轮新的中心点
    new_means = zeros(size(means));
    %分到每个中心点的像素个数
    num_assigned = zeros(k,1);
    
    for i=1:dim
        for j=1:dim
            %计算像素点到每个中心的距离
            r = A(j,i,1); g = A(j,i,2); b = A(j,i,3);
            diff = ones(k,1)*[r,g,b] - means;
            distance = sum(diff.^2, 2);
            %选择最近的中心点
            [val ind] = min(distance);
            %该像素点加到最近的中心点
            new_means(ind,1) = new_means(ind,1) + r;
            new_means(ind,2) = new_means(ind,2) + g;
            new_means(ind,3) = new_means(ind,3) + b;
            
            num_assigned(ind) = num_assigned(ind)+1;
        end
    end
    
    %计算新的中心点
    for i=1:k
        if (num_assigned(i) > 0)
            means(i,:) = new_means(i,:) ./ num_assigned(i);
        end
    end
    
    %检测是否收敛
    d = sum(sqrt(sum((means - new_means).^2, 2)));
    if d < 1e-5
        break
    end
end

disp(itr);

means = round(means);

large_image = double(imread('bird_large.tiff'));
large_dim = size(large_image,1);

for i=1:large_dim
    for j=1:large_dim
        r = large_image(j,i,1); g = large_image(j,i,2); b = large_image(j,i,3);
        diff = ones(k,1)*[r,g,b] - means;
        distance = sum(diff.^2, 2);
        [val ind] = min(distance);
        large_image(j,i,:) = means(ind,:);
    end
end

%显示
imshow(uint8(round(large_image))); hold off

%保存图像
imwrite(uint8(round(large_image)), 'bird_kmeans.tiff');

figure; hold on
for i=1:k
    col = (1/255).*means(i,:);
    rectangle('Position', [i,0,1,1], 'FaceColor', col, 'EdgeColor', col);
end
axis off

原始大鸟图：