Subgradient Method for K-Means Like Problem
You should use something like Robust Convex Clustering:

Assume that we are given $ n $ data points and each data point
is described by a $ d $ dimensional feature vector. We collectively
represent the data by a data matrix $ X \in \mathbb{R}^{d \times n} $.
The convex clustering method clusters data points into groups via solving
a convex optimization problem:

$$ \arg \min_{P} \frac{1}{2} \left\| X - P \right\|_{F}^{2} + \alpha \sum_{i &lt; j} {w}_{i, j} \left\| {P}_{i} - {P}_{j} \right\|_{p} $$

Where $ \left\| \cdot \right\|_{p} $ is the $ {L}_{p} $ norm, $ P \in \mathbb{R}^{d \times n} $ is the matrix consisting cluster centers and assignments, $ {P}_{i} $ is the $ i $ -th column of the matrix $ P $ and represents the centroid of cluster that the $ i $ -th data point is assigned, $ \alpha $ is a non-negative regularization parameter that controls the number of clusters. $ {w}_{i, j} $ is the weight between $ i $ -th and $ j $ -th data point specified by the user. Clusters are implicitly given by $ P $ matrix and $ {P}_{i} = {P}_{j} $ means the two data points are in the same cluster. When $ \alpha $ is zero, each data point forms a cluster since $ P $ equals to $ X $ when Eq. (1) reaches the
optimum. As $ \alpha $ increases, $ \left\| {P}_{i} − {P}_{j} \right\|_{p} $ becomes smaller for all pairs so that the whole function can reach the minimal value. Hence more data points will be grouped into a cluster.

Another great reference is Splitting Methods for Convex Clustering.

As the formulation you chose is not convex and of top of that I don't think anyone has ever calculates its Sub Gradient.

Update

I encountered another simple yet effective of K-Means - Revisiting K-Means: New Algorithms via Bayesian Nonparametrics.
It formulates the problem in as, intuitively, the Lagrange form of the K-Means with continuous parameter $ \lambda $ instead of the discrete parameter $ k $.