k-means聚类算法

k-means聚类算法主要思想源自于EM，整个过程分为两步：

第一步，固定需要聚类的中心点，然后计算所有需要进行聚类操作的点到中心点的距离；

第二步，利用上面的最小距离更新聚类的中心点，并且重复第一步。

当不再需要做出调整或者调整的幅度很小的时候就表明聚类成功了。

function [Cluster Codebook] = cvKmeans(X, K, stopIter, distFunc, verbose)
if ~exist('stopIter', 'var') || isempty(stopIter)
    stopIter = .05;
end
if ~exist('distFunc', 'var') || isempty(distFunc)
    distFunc = @cvEucdist;
end
if ~exist('verbose', 'var') || isempty(verbose)
    verbose = false;
end
[D N] = size(X);
if K > N,
    error('K must be less than or equal to the number of vectors N');
end
Codebook = X(:, randsample(N, K));

improvedRatio = Inf;
distortion = Inf;
iter = 0;
while true
    d = distFunc(Codebook, X);
    [dataNearClusterDist, Cluster] = min(d, [], 1);
    old_distortion = distortion;
    distortion = mean(dataNearClusterDist);
    improvedRatio = 1 - (distortion / old_distortion);
    if verbose
        fprintf('%d: improved ratio = %f\n', iter, improvedRatio);
    end
    iter = iter + 1;
    if improvedRatio <= stopIter, break, end;
    for i=1:K
        idx = find(Cluster == i);
        Codebook(:, i) = mean(X(:, idx),2);
    end
end

上面是聚类函数的定义，函数需要最多五个参数，如果参数数目不够则会使用默认的参数数值。其中最主要的是第一个，这个是用于聚类的矩阵，而第二个是聚类之后类别的个数。函数的开头首先随机得到一个初始的中心点，这个随机点的实际上是一个大小为D*K的矩阵。其中前面的D表示这个中心点的坐标表示是在D维空间上的，而后面的K则是表示需要K个这样的中心点。经过上面的分析可以猜测经过distFunc函数计算的距离矩阵应该是K*N，其中K表示需要聚合的K个类，而N则表示需要被聚类的N个点。

所以下面的min函数在d的第一维里面进行比较，也就是每一列里面比较得到最小的值，并且将相应的列的索引存放在Cluster中。通过mean求得最小值队列的平均值，然后进行调整的更新，最后一步是更新中心点。

function d = cvEucdist(X, Y)
if ~exist('Y', 'var') || isempty(Y)
    %% Y = zeros(size(X, 1), 1);
    U = ones(size(X, 1), 1);
    d = abs(X'.^2*U).'; return;
end
V = ~isnan(X); X(~V) = 0; % V = ones(D, K); 
U = ~isnan(Y); Y(~U) = 0; % U = ones(D, N); 
d = abs(X'.^2*U - 2*X'*Y + V'*Y.^2);

在上面的默认参数里面，实际上distFunc是一类似于C语言里面函数指针的概念，这个值被赋值为cvEucdist函数，所以实际的调用实在cvEucdist函数。

根据之前的分析返回的d应该是一个K*N的矩阵，而欧氏距离又是（x-y）^2。如果对这个距离公式进行展开的话就是x^2-2*x*y+y^2。这里最后一句分开来看就是上面的这种展开。X'.^2将矩阵内的所有元素原地平方——也就是在不改变矩阵X的大小的情况下，对每一个X中的元素进行平方处理，不过到目前为止X的大小还是为K*D，而根据欧氏距离的扩展公式，应该是（x1-x2)^2+(y1-y2)^2。所以还需要一个对矩阵的调整使得X'^2可以表示第一项，也就是需要一个累加过程，这个过程是利用U进行的，所以矩阵变为K*N。对照这个推理，第二个和第三个矩阵的变化也就很好理解。