kmeans之於模式識別，如同“hello world”之於C、之於任何一門高階語言。

演算法的規格（specification）

在聚類問題（一般非監督問題）中，給定訓練樣本X={x(1),x(2),…,x(N)}，每個x(i)∈Rd。kmeans演算法的職責在於將這N個樣本聚類成k個簇（cluster, μ1,μ2,…,μk），流程如下：

1. 隨機選取k個聚類中心(cluster centroids)為μ1,μ2,…,μk
   C = X(randperm(m*n, k), :); # 程式語言

2. 重複一下過程直至收斂
   {
   對於每一個樣本i

   ，根據最近鄰（歐氏距離度量）計算其所屬分類
   
   c(i):=argminj∥x(i)−μj∥2
   對於每一個類j，重新計算該類的質心（centroids）
   μj:=∑mi=11{c(i)=j}x(i)∑mi=11{c(i)=j}
   }

演算法的規格：

• 一個引數k，聚類中心的數目，當然也有一些常規的引數，比如最大迭代次數epochs，容忍度tol
• 一個迴圈，判斷目標函式是否變化足夠小，以F範數（Frobenius norm）為度歸。

while true,
    ...
    if norm(J_cur-J_prev, 'fro'
 
) < tol,
        break;
    end
    J_prev = J_cur;
end

• 一條更新語句，更新各個類的聚類中心，根據每個樣本應屬的類別（歐式距離最小表徵）

這個公式看似高大上，實則不值一提，翻譯過來就是新的聚類中心（centroid）在該類別空間的中心處。

    dist = sum(X.^2, 2)*ones(1, k) + (sum(C.^2, 2)*ones(1, m*n))'...
        - 2*X*C';
    [~, idx] = min(dist, []
 
, 2) ;
    for i = 1:k,
       C(i, :) = mean(X(idx == i , :)); # 對應於這樣一條語句
    end

matlab實現

客戶端（client）程式

clear all; close all;
I = imread('./lena.bmp');
[m, n, p] = size(I);
k = 7;
[C, label, J] = kmeans(I, k);
I_seg = reshape(C(label, :), m, n, p);
figure
subplot(1, 2, 1), imshow(I, []), title('原圖')
subplot(1, 2, 2), imshow(uint8(I_seg), []), title('聚類圖')
figure
plot(1:length(J), J), xlabel('#iterations')

kmeans函式

function [C, label, J] = kmeans(I, k)
[m, n, p] = size(I);
X = reshape(double(I), m*n, p);
rng('default');
C = X(randperm(m*n, k), :);
J_prev = inf; iter = 0; J = []; tol = 1e-2;
while true,
    iter = iter + 1;
    dist = sum(X.^2, 2)*ones(1, k) + (sum(C.^2, 2)*ones(1, m*n))' - 2*X*C';
    [~, label] = min(dist, [], 2) ;
    for i = 1:k,
       C(i, :) = mean(X(label == i , :));
    end
    J_cur = sum(sum((X - C(label, :)).^2, 2));
    J = [J, J_cur];
    display(sprintf('#iteration: %03d, objective fcn: %f', iter, J_cur));
    if norm(J_cur-J_prev, 'fro') < tol,
        break;
    end
    J_prev = J_cur;
end

實驗結果

目標函式收斂情況

目標函式

J_cur = sum(sum((X - C(label, :)).^2, 2));

效果圖