# K-means clustering

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}}

• Have you tried typing "k means" into the documentation browser? The second hit is this page. Commented Jun 25, 2012 at 17:00
• @Szabolcs not on version 7. Commented Jun 25, 2012 at 18:49
• @Mr.Wizard You're right, that's a new in 8 function. But Google gives it as a first hit as well. Commented Jun 25, 2012 at 18:52
• @Szabolcs I'm subtly suggesting that a version-7 solution would make this question more interesting. ;-) Commented Jun 25, 2012 at 19:05
• @Mr.Wizard Sure, a question is what you make out of it! :-) Commented Jun 25, 2012 at 19:28

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.

• Thank you for this nice answer. Unfortunately, you were a bit faster than me (I had problems to recognize why I got in my first attempts {1,1,1,2,2} and not {1,1,2,2,2}). Great that you give an explanation for this too! Commented Jun 25, 2012 at 18:33
• Thank you very much dear Sjoerd C. de Vries Commented Jun 25, 2012 at 19:41

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.

• Cool simple implementation of K-means. Here is a nice post about K-means Lloyd's algorithm, using Mathematica. Commented May 22, 2014 at 2:11

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.

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]
];