k-means (after Lloyd1957)

Input: data set 𝒜 ⊂ (ℝ^n),  concatenation of centers  CC; Using LS distance-like function;
Output: Vector a, concatenation of optimal centers, objective function value;
If Ind≠0then the intermediate steps and results are shown
Submodules: Fmin[]

kMeansLoyd[A_,CC_,Ind_]:=Module[{m=Length[A],k=Length[CC],it=0, IT=50,cen=CC,a,C1,n1,eps=0.005},
F0=Fmin[A,CC];
Print[Style["it ",Blue,14],Style[it,Blue,14],": CC=",CC//N,"; Fmin(CC)=",F0//N];
While[
it=it+1; 
a=Table[0,{i,Length[A]}];
Do[
a[[i]]=1;
Do[
If[Norm[A[[i]]-cen[[j]]]< Norm[A[[i]]-cen[[a[[i]]]]],a[[i]]=j];
,{j,k}]
,{i,m}];
			(* Determining the new center *)
C1=n1=Table[0,{j,k}];
Do[
C1[[a[[i]]]]=C1[[a[[i]]]]+A[[i]];
n1[[a[[i]]]]=n1[[a[[i]]]]+1
,{i,m}];
Do[
C1[[j]]=C1[[j]]/n1[[j]]
,{j,k}];
F1=Fmin[A,C1]//N;
Err=(F0-F1)/F0;
If[Ind≠0,Print[it,". C1=",C1//N,"; Fmin=",F1,"; Err = ",Err]];
Err>eps && it<IT,
cen=C1;F0=F1;
];
{a,C1,F1,it}
]

Example

Modules

Needs["ComputationalGeometry`"]

kMeansLoyd[A_,CC_,Ind_]:=Module[{m=Length[A],k=Length[CC],it=0, IT=50,cen=CC,a,C1,n1,eps=0.005},
F0=Fmin[A,CC];
Print[Style["it ",Blue,14],Style[it,Blue,14],": CC=",CC//N,"; Fmin(CC)=",F0//N];
While[
it=it+1; 
a=Table[0,{i,Length[A]}];
Do[
a[[i]]=1;
Do[
If[Norm[A[[i]]-cen[[j]]]< Norm[A[[i]]-cen[[a[[i]]]]],a[[i]]=j];
,{j,k}]
,{i,m}];
			(* Determining the new center *)
C1=n1=Table[0,{j,k}];
Do[
C1[[a[[i]]]]=C1[[a[[i]]]]+A[[i]];
n1[[a[[i]]]]=n1[[a[[i]]]]+1
,{i,m}];
Do[
C1[[j]]=C1[[j]]/n1[[j]]
,{j,k}];
F1=Fmin[A,C1]//N;
Err=(F0-F1)/F0;
If[Ind≠0,Print[it,". C1=",C1//N,"; Fmin=",F1,"; Err = ",Err]];
Err>eps && it<IT,
cen=C1;F0=F1;
];
{a,C1,F1,it}
]

Fmin[A_,c_]:=Module[{m=Length[A],k=Length[c]},
Sum[Min[Table[Norm[A[[i]]-c[[j]]]^2,{j,k}]],{i,m}]
]

Data1

A={{1,9},{2,9},{2,6},{1,3},{5,3},{6,4},{4,6},{7,7},{8,6},{9,8}};
CC={{2,6},{5,4},{6,4}};
AG={{0,10},{0,10}};
vor0=DiagramPlot[CC];
Print[Style["Data set A, Initial centers and the Voronoi diagram:",Blue,14]];
slCC=ListPlot[CC,PlotStyle->{Red,PointSize[.06],Opacity[.5]}];
Print[Show[ListPlot[A,PlotStyle->{PointSize[.02]}],slCC,vor0,AspectRatio->Automatic,GridLines->Automatic,Frame->True,PlotRange->AG,ImageSize->Small]]

Data2

SeedRandom[13]
m=200; k=5;
A=RandomReal[{0,10},{m,2}];
CC=RandomReal[{0,10},{k,2}]; 
Print[Style["Data set A, Initial centers and the Voronoi diagram:",Blue,14]];
AG={{0,10},{0,10}};
vor0=DiagramPlot[CC];
slCC=ListPlot[CC,PlotStyle->{Red,PointSize[.06],Opacity[.5]}];
Print[Show[ListPlot[A,PlotStyle->{PointSize[.02]}],slCC,vor0,AspectRatio->Automatic,GridLines->Automatic,Frame->True,PlotRange->AG,ImageSize->Small]]

Execution

Ind=0;km=kMeansLoyd[A,CC,Ind];
Print[Style["it ",Blue,14],Style[km[[4]],Blue,14],Style[". Optimal centroids:\n",Blue,14],km[[2]]//N,Style["\nObjective LS-function value = ",Blue,14],km[[3]]//N];
				(* Optimal partition PI1 *)
PI1=Table[{},{j,Length[km[[2]]]}];
Do[
Do[
	If[km[[1,i]]==j,PI1[[j]]=Append[PI1[[j]],A[[i]]]] 
,{j,Length[km[[2]]]}]
,{i,Length[A]}];
					(* Figure *)
boja={Blue,Green,Brown,Red,Yellow,Orange,Purple};
vor1=DiagramPlot[km[[2]]];
slC1=ListPlot[km[[2]],PlotStyle->{Red,PointSize[.06],Opacity[.5]}];
(* ListPlot[PI1,PlotStyle->{PointSize[.035]}] *)
Print[Show[Table[ListPlot[PI1[[j]],PlotStyle->{boja[[j]],PointSize[.035]}],{j,Length[km[[2]]]}],slC1,vor1,AspectRatio->Automatic,GridLines->Automatic,Frame->True,PlotRange->AG,ImageSize->Small]]
If[Ind≠0,Print["PI1:\n",PI1]]