K-means clustering

Question

In MATLAB, there is a command kmeans() that divides an array into $k$ clusters and calculates the centroid of each cluster. Is there any command in Mathematica to perform the same action?

For example:

x = {{1,2,3,4,5}, {6,7,8,9,10}, {11,12,13,14,15}, {16,17,18,19,20}, {21,22,23,24,25}}

The MATLAB kmeans command does this:

[b, c] = kmeans(x, 2)

It divides x into two clusters and calculates the centroids of these two clusters, and indicates for each element in the array which cluster it is a member of

b = {1,1,2,2,2}

and

c={{3.5000,4.5000,5.5000,6.5000,7.5000},{16.0000,17.0000,18.0000,19.0000,20.0000}}

Answer 1

ClusteringComponents is indeed the function to go for. To get the same results as MATLAB you need to do the following:

x = {{1, 2, 3, 4, 5}, {6, 7, 8, 9, 10}, {11, 12, 13, 14, 
    15}, {16, 17, 18, 19, 20}, {21, 22, 23, 24, 25}};

cc = ClusteringComponents[x, 2, 1, Method -> "KMeans", 
            "DistanceFunction" -> SquaredEuclideanDistance, "RandomSeed" -> 1]

{1, 1, 2, 2, 2}

The arguments x and 2 are the same as with MATLAB. The '1' is used to indicate the level of the nested array to consider the data points. In this case we're looking at the top level, so we're considering 5D points.

Mathematica can use various methods. "KMeans" is default, so it isn't necessary to provide it here. The default "DistanceFunction" is EuclideanDistance. MATLAB's is SquaredEuclideanDistance, so we have to explicitly use that.

Since clustering uses a process with random initializations the results may differ depending on the RNG state. I used "RandomSeed" -> 1 to initialize the RNG to a state that yields the results you showed. {1,1,1,2,2} is a possible output too.

Now to the centroids:

Mean /@ {Pick[x, cc, 1], Pick[x, cc, 2]} // N

{{3.5, 4.5, 5.5, 6.5, 7.5}, {16., 17., 18., 19., 20.}}

I added //N since you seemed to want machine precision results. Leave it away for exact results.

---

On a side note: KMeans may sometimes yield disastrous results. It's a well-known property of this algorithm.

x1 = RandomVariate[MultinormalDistribution[{0, 0}, {{1, 0}, {0, 20}}],500];
x2 = RandomVariate[MultinormalDistribution[{6, 0}, {{1, 0}, {0, 20}}],500];
Graphics@{Red, Point@x1, Green, Point@x2}

xx = Join[x1, x2];
cc = ClusteringComponents[xx, 2, 1, 
        Method -> "KMeans", 
        "DistanceFunction" -> SquaredEuclideanDistance, "RandomSeed" -> 1];
{c1, c2} = {Pick[xx, cc, 1], Pick[xx, cc, 2]};
Graphics@{Red, Point@c1, Green, Point@c2}

In this case, one of the additional three clustering methods that Mathematica knows (Method -> "PAM") works wonders.

Answer 2

Since Mr.Wizard mentions that ClusteringComponents is unavailable in Mathematica 7, here's an implementation of Lloyd's algorithm for k-means clustering (can also be interpreted as an Expectation-Maximization approach) that will run on version 7.

Clear[kmeans]
kmeans[list_, k_, opts : OptionsPattern[
    {DistanceFunction -> SquaredEuclideanDistance, "RandomSeed" -> {}}]] :=

BlockRandom[SeedRandom[OptionValue["RandomSeed"]];
    Module[{m = RandomSample[list, k], update, partition, clusters},
        update[] := m = Mean /@ clusters;
        partition[_] := (clusters = GatherBy[list, RandomChoice@
            Nearest[m, #, (# -> OptionValue[#] &@DistanceFunction)] &]; update[]);
        FixedPoint[partition, list];
        {clusters, m}
    ]
]

Use it as the following:

{clusters, centroids} = kmeans[x, 2] // N;
(* clusters = {{{1., 2., 3., 4., 5.}, {6., 7., 8., 9., 10.}}, 
               {{11., 12., 13., 14., 15.}, {16., 17., 18., 19., 20.}, {21., 22., 23., 24., 25.}}}

   centroids = {{3.5, 4.5, 5.5, 6.5, 7.5}, {16., 17., 18., 19., 20.}} *)

You can also use a custom distance function using the DistanceFunction option (default is SquaredEuclideanDistance) or use a particular seed for the PRNG using the "RandomSeed" option to get reproducible results. Note that this implementation does not do any kind of error checking. Typically, k-means algorithms are run for a several different initial states to verify that the global optimum is reached (exceptions exist) and the above can be easily modified, if you so choose to, to incorporate these.

Answer 3

The code below is very instructive. I get from Sjoerd C. de Vries's comment here.

GraphicsGrid[

 Table[
  data1 = Table[{RandomReal[]/5, RandomReal[]}, {100}];
  data2 = Table[{RandomReal[]/5 + distance, RandomReal[]}, {100}];
  dataM = Flatten[{data1, data2}, 1] // RandomChoice[#, Length[#]] &;

  Table[clusters = 
    ClusteringComponents[dataM, 2, 1, 
     DistanceFunction -> EuclideanDistance, Method -> myMethod];
   Graphics[{Green, Point@Pick[dataM, clusters, 1], Red, 
     Point@Pick[dataM, clusters, 2]}, 
    PlotLabel -> myMethod], {myMethod, {"Agglomerate", "Optimize", 
     "KMeans", "PAM"}}]

  , {distance, 0.2, 0.8, 0.1}]

 , ImageSize -> 500, Frame -> All]

I miss k-means in FindClusters function, maybe in Mathematica 9.

Answer 4

In MATLAB, the command kmeans() also returns distances from data points to centroids, i.e.

[idx,C,sumd,D] = kmeans(___)

So, a more complete answer would be:

kmeans[data_, k_] :=
  Module[
   {idx , cntr, result, noOfCols, noOfRows, c, i, y, sumd, dist},
   idx = ClusteringComponents[data, k, 1, Method -> "KMeans",
        "DistanceFunction" -> SquaredEuclideanDistance,
        "RandomSeed" -> 1];
   noOfRows = Dimensions[data][[1]]; 
   noOfCols = Dimensions[data][[2]];
   cntr = ConstantArray[0, {k, noOfCols}];
   sumd = ConstantArray[0, {k, 1}];
   For[c = 1, c <= k, c++,
    y = Pick[data, idx, c];
    If[ y != {},
     cntr[[c]] = Mean[y];
     ];
    For[i = 1, i <= Length[y], i++,
     sumd[[c]] += SquaredEuclideanDistance[ y[[i]], cntr[[c]]];
     ];
    ];
   dist = ConstantArray[0, {noOfRows, k}];
   For[i = 1, i <= noOfRows, i++,
    For[c = 1, c <= k, c++,
      dist[[i, c]] = SquaredEuclideanDistance[
         data[[i]], cntr[[c]] ;
         ];
      ];
    ];
   result["idx"] = idx;
   result["cntr"] = cntr;
   result["sumd"] = sumd;
   result["dist"] = dist;
   Return[result] 
   ];