How do you test an implementation of k-means? Disclaimer: I posted this question on Stackoverflow, but I thought maybe this is better suited for this platform.
How do you test your own k-means implementation for multidimensional data sets?
I was thinking of running an already existing implementation (i.e., Matlab) on the data and compare the results with my algorithm. But this would require both algorithms to work more than roughly the same, and the mapping between the two results probably is no piece of cake.
Do you have a better idea?
 A: The mapping between two sets of results is easy to compute, because the information you obtain in a test can be represented as a set of three-tuples: the first component is a (multidimensional) point, the second is an (arbitrary) cluster label supplied by your algorithm, and the third is an (arbitrary) cluster label supplied by a reference algorithm.  Construct the $k$ by $k$ classification table for the label pairs: if the results agree, it will be a multiple of a permutation matrix.  That is, each row and each column must have exactly one nonzero cell.  That's a simple check to program.  It's also straightforward to track small deviations from this ideal back to individual data points so you can see precisely how the two answers differ if they differ at all.  I wouldn't bother to compute statistical measures of agreement: either there is perfect agreement (up to permutation) or there is not, and in the latter case you need to track down all points of disagreement to understand how they occur.  The results either agree or they do not; any amount of disagreement, even at just one point, needs checking.
You might want to use several kinds of datasets for testing: (1) published datasets with published k-means results; (2) synthetic datasets with obvious strong clusters; (3) synthetic datasets with no obvious clustering.  (1) is a good discipline to use whenever you write any math or stats program.  (2) is easy to do in many ways, such as by generating some random points to serve as centers of clusters and then generating point clouds by randomly displacing the cluster centers relatively small amounts.  (3) provides some random checks that potentially uncover unexpected behaviors; again, that's a good general testing discipline.
In addition, consider creating datasets that stress the algorithm by lying just on the boundaries between extreme solutions.  This will require creativity and a deep understanding of your algorithm (which presumably you have!).  One example I would want to check in any event would be sets of vectors of the form $i \mathbb{v}$ where $\mathbb{v}$ is a vector with no zero components and $i$ takes on sequential integral values $0, 1, 2, \ldots, n-1$.  I would also want to check the algorithm on sets of vectors that form equilateral polygons.  In either situation, cases where $n$ is not a multiple of $k$ are particularly interesting, including where $n$ is less than $k$.  What is common to these situations is that (a) they use all the dimensions of the problem, yet (b) the correct solutions are geometrically obvious, and (c) there are multiple correct solutions.
(Form random equilateral polygons in $d \ge 2$ dimensions by starting with two nonzero vectors $\mathbb{u}$ and $\mathbb{v}$ chosen at random.  (A good way is to let their $2d$ components be independent standard normal variates.)  Rescale them to have unit length; let's call these $\mathbb{x}$ and $\mathbb{z}$.  Remove the $\mathbb{x}$ component from $\mathbb{z}$ by means of the formula
$$\mathbb{w} = \mathbb{z} - ( \mathbb{z} \cdot \mathbb{x} ) \mathbb{x}.$$  
Obtain $\mathbb{y}$ by rescaling $\mathbb{w}$ to have unit length.  If you like, uniformly rescale both $\mathbb{x}$ and $\mathbb{y}$ randomly.  The vectors $\mathbb{x}$ and $\mathbb{y}$ form an orthogonal basis for a random 2D subspace in $d$ dimensions.  An equilateral polygon of $n$ vertices is obtained as the set of $\cos(2 \pi k / n) \mathbb{x} + \sin(2 \pi k / n) \mathbb{y}$ as the integer $k$ ranges from $0$ through $n-1$.)
A: One very simple 'naive' approach would be to use simple synthetic data, for that every implementation should result in the same clusters.
Example in Python with import numpy as np:
test_data = np.zeros((40000, 4))
test_data[0:10000, :] = 30.0
test_data[10000:20000, :] = 60.0
test_data[20000:30000, :] = 90.0
test_data[30000:, :] = 120.0

For n_clusters = 4 it should give you a permutation of [30, 60, 90, 120]
A: The k-means includes a stochastic component, so it is very unlikely you will get the same result unless you have exactly the same implementation and use the same starting configuration. However, you could see if your results are in agreement with well-known implementations (don't know about Matlab, but implementation of k-means algorithm in R is well explained, see Hartigan & Wong, 1979). 
As for comparing two series of results, there still is an issue with label switching if it is to be run multiple times. Again, in the e1071 R package, there is a very handy function (;matchClasses()) that might be used to find the 'best' mapping between two categories in a two-way classification table. Basically, the idea is to rearrange the rows so as to maximise their agreement with columns, or use a greedy approach and permute rows and columns until the sum of on the diagonal (raw agreement) is maximal. Coefficient of agreement like the Kappa statistic are also provided.
Finally, about how to benchmark your implementation, there are a lot of freely available data, or you can simulate a dedicated data set (e.g., through a finite mixture model, see the MixSim package).
A: Since k-means contains decisions that are randomly chosen (the initialization part only), I think the best way to try your algorithm is to select the initial points and let them fixed in your algorithm first and then choose another source code of the algorithm and fix the points in the same way. Then you can compare for real the results.
