What if k-means starts in a local minimum? I have to find 10 clusters of 100 samples with dimension 100. I have access to two k-means implementations. Both of them initialize the means with 10 randomly picked samples. When I run these algorithms on my samples, I can see that they immediately converge: After one iteration, the solution is stable with the following outcome:


*

*The samples which coincide with the means get assigned to these clusters

*All other (90) samples are assigned to one cluster, e.g. cluster no 8.


How can this be? Is there an obvious (common) bug in these two implementations or is it simply the consequence of some special property of the samples? I know that the samples are all orthogonal to each other ($X = V \cdot \Lambda^{-1/2}$, where $V$ and $\Lambda$ are eigenvectors and -values of another matrix), but I don't see a connection.
Here is the more readable implementation (in MATLAB). Note that on a simple 2-dimensional toy example, both implementations work fine.
function [z, means] = kmeans(X, k)
    % input: X: (d x n) n observations in columns of size d
    %        k: the number of clusters

    % output: z: (k x 1) assignment vector
    %         means: (d x k) matrix whose columns are mean vectors
    [d n] = size(X);

    max_iter = 300;
    tol = 1e-4;

    means = X(:, randperm(n)(1:k));
    z = zeros(n, 1);

    iter = 0;
    change = inf;
    err = inf;
    while (change > tol && iter < max_iter)
        % assignment step
        for i = 1:n
            diff = means - repmat(X(:,i), 1, k);
            [_, best_mean] = min(norm(diff, 'cols'));
            z(i) = best_mean;
        end

        % update step
        for j = 1:k
            ind = (z == j) ./ sum(z == j);
            means(:, j) =  X * ind;
        end

        % update error
        old_err = err;
        err = 0;
        for i = 1:n
            err = err + norm(means(:, z(i)) - X(:, i), 'cols');
        end

        change = old_err - err;
        iter = iter + 1;
    end
end

 A: In the first iteration :
Let $y$ be one of the points considered to be the mean in the first iteration. Let $x$ be one of the remaining 90 points. 
Now, $||x - y||^2 = ||x||^2 + ||y||^2 - 2\langle x, y \rangle$. However, since $x$ and $y$ are orthogonal, $||x - y||^2 = ||x||^2 + ||y||^2$. So, when I try to assign $x$ to one of the cluster means, the point will simply be assigned to the cluster mean with the smallest norm. This assignment is independent of which $x$ is under consideration. Therefore, all the remaining 90 points will be assigned to the mean with the smallest norm.
So, in the first iteration, all the 90 points are simply assigned to the cluster mean with smallest norm. Let us call this cluster, cluster8 (as in your example).
All the remaining cluster means remain unchanged because no other points are assigned to them, but cluster8's mean will change as it has 91 points in it. Let $z$ be the recomputed mean of cluster8. Note that, $$z = \frac{1}{91}\sum_{i = 1}^{91} x_i$$ where $x_i$ are the points assigned to cluster8.
Now, in the second iteration:
Note that,
$$\langle x, z \rangle = \frac{1}{91}\sum_{i = 1}^{91} \langle x, x_i \rangle$$
Since, all the points are orthogonal to each other, the above expression either will evaluate to 0, if $x \neq x_i \forall i$ or it will evaluate to $\frac{||x||^2}{91}$ if $\exists i$ s.t. $x = x_i$.
The $x$'s which are not equal to any $x_i$ are exactly the other cluster means and they will be simply be assigned to their own cluster. All the other $x$ will equal some $x_i$ and therefore,$$||x - z||^2 = ||x||^2 + ||z||^2 - 2\langle x,z \rangle = ||x||^2 + ||z||^2 - \frac{2}{91}||x||^2 = \frac{89}{91}||x||^2 + ||z||^2$$.
Therefore, all those $x$'s that have been assigned to cluster8 will only change their mean to some other cluster if $z$ does not continue to be the smallest mean amongst all cluster means. Actually, even if $z$ is not the smallest but slightly larger than the smallest, even then the $x$'s will continue to remain in $z$ because of the $\frac{89}{91}$ factor.
A: Are you by chance using high-dimensional data, or data where Euclidean distance might not be appropriate? If so, stay away from k-means.
A good example is binary data. The distance between two instances will be the square root of the number of bit differences. So assuming you have random vectors, the expected distance will be $\sqrt{\frac{d}{2}\cdot 1}$. At the same time, all of them will be at a distance of $\sqrt{d\cdot\frac{1}{4}}$ from the center vector at $\{0.5\}^d$.
Now if you use MacQueen style iteration, the first 10 objects will be the seeds. The 11th object will be assigned to one of the others, moving this mean somewhere into the middle. From now on, this mean will likely be closer to any of the remaining objects than any of the other means and likely collect all the remaining points. Lloyd style iteration is not much better, it may just take a bit longer for the "most central" mean to surface.
To a certain amount this also applies to any discretized data. Say Fisher's iris data set, which has a measurement resolution of $0.1$. You can actually compute the estimated artifacts from this. I'm not going to say it breaks k-means on Fisher's iris data, but it does have an effect. And you need to be aware of this.
Bottom line is: don't use k-means on data that has a small number of discrete values. Because the mean will yield inappropriate values for this. You might want to try k-medians, which should work slightly better here. However, IMHO, you might still be using an inappropriate combination of data set, distance function and algorithm. Definitely make sure that the distance function returns useful values!
You can try ELKI which has MacQueen and Lloyd style iterations, and also k-means++ seeding, if it performs any better for you. But if you are using inappropriate data, this will not help either. I'm not sure if they also have k-medians.
Oh, and you might want to update your question, and share some more information on your data set. Because all of this reply is based on guessing what your data might be like.
