K-Mean Clustering Algorithm

Question

I have implemented the k-Means Clustering Algorithm using the Wolfram Language. However, I think it can be more efficient. Do you have any idea how to make it more efficient (e.g., by removing the Table[], Do[], AppendTo[])?

Many thank!

kMeans[data_, k_] := Module[{x = N[data], c, cn, n = Length[data], min, clusters, sum, 
  t = 0},
 (* Randomly choose k centriods from the data points *)
 
 c = RandomSample[x, k];
 Do [
  t++;
  Print["iteration # ", t];
  (* Initialize the clusters as empty *)
  clusters = Table[{}, k];
  (* Assign each point to the cluster of the closest centriod *)
  
  Table [
   d = Map[EuclideanDistance[#, x[[i]]] &, c];
   min = Ordering[d, 1][[1]];
   AppendTo[clusters[[min]], x[[i]]];
   , {i, 1, n}]; (* End do *)
  Print["Centroids = ", c];
  (* Calculate the new centriods *)
  cn = Mean /@ clusters;
  (* If the centriods are the "almost" same then terminate*)
  
  If [ Select[
     Flatten[
      MapThread[Abs[#1 - #2] &, { c, cn}]], # > 0.001 &] == {} , 
   Break[], c = cn];
  , {10 (* Maximum number of iterations *)}];
 {clusters, c}
 ]

Answer 1

Compiled, parallelized ordering plus a more functional approach. The switch between DistanceMatrix and Outer (for dm) may be system dependent.

ordC = Compile[{{x, _Real, 1}},
   First@Ordering[x, 1],
   RuntimeAttributes -> {Listable}, Parallelization -> True];

kMeans2 // ClearAll;
kMeans2[data_, k_] := Module[{x = Developer`ToPackedArray@N[data], dm},
  If[ArrayDepth[x] > 1 || k*Length[x] < 10^5,
   dm =(*Sqrt@*)Total[Outer[Subtract, ##, 1]^2, {3}] &,
   dm = DistanceMatrix;
   ];
  NestWhile[
   With[
     {clusters = 
       With[{minpos = ordC[dm[x, Last@#]]}, 
        Pick[x, minpos, #] & /@ Range@k]},
     {clusters, Mean /@ clusters}
     ] &,
   (* initial {clusters, means} *)
   {{}, RandomSample[x, k]},
   Norm[Last[#1] - Last[#2], Infinity] - 0.001 > 0 &,
   2,
   100 (*Maximum number of iterations*)]
  ];

Example:

foo = RandomReal[10, {10^3, 3}];
(SeedRandom[0];
  km1 = kMeans[foo, 8];) // RepeatedTiming
(SeedRandom[0];
  km2 = kMeans2[foo, 8];) // RepeatedTiming
km1 == km2
(*
{0.212839, Null}
{0.0588948, Null}
True
*)