Why isn't k-means optimized using gradient descent? I know k-means is usually optimized using Expectation Maximization. However we could optimize its loss function the same way we optimize any other! 
I found some papers that actually use stochastic gradient descent for large scale k-means, but I couldn't get my question answered.
So, do anyone know why is that? Is it because Expectation Maximization converges faster? Does it has any particular guarantee? Or is it a historical reason?
 A: K-means clustering is unsupervised, and the closest unsupervised technique which uses EM is model-based clustering (Gaussian mixture models, GMM).  An annoying problem with GMM model-based clustering occurs when many of the features are correlated, which causes near-singularity in the feature-based covariance(correlation) matrix.  In this situation, the likelihood function becomes unstable, with condition indexes reaching infinity, causing 
GMM to break down completely.  
Thus, drop the idea of EM and kNN -- since it's based on covariance (correlation) matrices for unsupervised analysis.  Your inquiry on optimization closely resembles Sammon mapping, and classical metric and non-metric multidimensional scaling (MDS).  Sammon mapping is derivative-iterative based, while various forms of MDS are commonly iterative or one-step eigendecompositions, which can nevertheless optimize during a one-step matrix operation.  
Looking again back at your request: the answer is: it's already been done in Sammon mapping.     
A: As the OP mentions, it's possible to solve k-means using gradient descent, and this may be useful in the case of large scale problems.
There are certainly historical reasons for the prevalence of EM style algorithms for solving k-means (i.e. Lloyd's algorithm). Lloyd's algorithm is so popular that people sometimes call it "the k-means algorithm", and may even be unaware that other approaches exist. But, this popularity is not undeserved.
Bottou and Bengio (1995) showed that Lloyd's algorithm is equivalent to optimizing the k-means cost function using Newton's method. In general optimization problems, second order methods like Newton's method can converge faster than first order methods like gradient descent because they exploit information about the curvature of the objective function (and first order methods don't). In an experiment on the well-known Iris dataset, they showed that Lloyd's algorithm did indeed converge faster than gradient descent. It would be interesting to see this comparison on a wider variety of datasets.
References:
Bottou and Bengio (1995). Convergence properties of the k-means algorithms.
