coming from the computing science side rather than from the data analysis one, I studied the k-means clustering problem for a short time and noticed that the NP-hardness of the problem for $k=2$ seems to be separately proved in various places (publications as well as some online PDFs). For instance, the article NP-hardness of Euclidean sum-of-squares clustering (link found on Wikipedia) proves the result for $k=2$ only.
What I don't understand well is that proving such a result seems to lead to an easy one-sentence generalization (see below) to any value of $k$. Knowing very little about the k-means problem, I am pretty sure I am missing something but why not add the following claim to all these articles?
Having proved the NP-hardness of the k-means problem for $k=2$, we notice that any instance of the problem for some value of $k$ can be reduced to an instance of the same problem for the case $k+1$ (by adding one more point to the dataset far enough from all existing points), we easily prove the NP-hardness of the problem for any value of $k$.
For instance, say I want to solve the problem for some $k$ and $n$ points in $[0,1]$. I know the sum of the square is smaller than $n$. I now add the new point $1+\sqrt{2n}+\epsilon$ to the dataset, building an instance of the clustering problem for $k+1$. The solution must put this new point alone in a separate cluster (thus not changing the previous sum of squares); because in any other case the new sum of squares would be greater than $n$. Thus solving this new problem for $k+1$ clusters immediately solves the initial instance of the problem. Of course, the very same idea can be used for the plane or any euclidean space.