Any Solution for $k$-means with minimum and maximum cluster size constraint?

I am looking for an efficient approach to $$k$$-means clustering with minimum cluster size constraints. The clusters are non overlapping, so, one point can belong to only one cluster.

$$N$$ be the number of points

$$C$$ be the number of clusters

$$S_{\min}$$ be the minimum cluster size

$$S_{\max}$$ be the maximum cluster size

Let $$x_{nc}$$ be an indicator variable. If $$x_{nc}=1$$, then point $$n$$ belongs to cluster $$c$$, otherwise not.

So, the constraints we can express as

$$S_{\min}\le\sum_{n=1}^Nx_{nc}\le S_{\max},\forall c, c=1,2,\cdots,C$$

$$\sum_{c=1}^Cx_{nc}=1,\forall n, n=1,2,\cdots,N$$

$$\sum_{c=1}^C\sum_{n=1}^Nx_{nc}=N$$

What should be the objective here.

Or any other efficient mathematical formulations for this clustering problem?

We have the objective function for Kmean as

$$\min \sum_{c=1}^C\sum_{n\in \mathcal{S}_c}||d_n-\mu_c||^2$$

where, $$\mu_c$$ (a two-dimensional vector) is the mean of points in cluster $$c$$, $$\mathcal{S}_c$$

Here, $$d_n$$ is a two-dimensional vector containing the coordinates of point $$n$$.

• The last constraint is not mandatory. With the second constraint, it will always be satisfied. May 5, 2021 at 9:00
• @Kuifje, yes, of course. May 5, 2021 at 9:06
• @dipaknarayanan If you are interested in tackling your problem in a model-and-run fashion while keeping speed and scalability, have a look at localsolver.com/docs/last/exampletour/kmeans.html. Disclaimer: LocalSolver is commercial software, but trial and academic licenses are available for free. May 14, 2021 at 7:23

It depends on what you want to achieve exactly. You could, for example, want to minimize the average distance between a point and the "center" of the cluster. This requires that you first define this center. Assuming it is one of the points, then intoduce a binary variable $$y_j$$ that takes value $$1$$ if and only if point $$j$$ is the center of one of the clusters. And slightly modifiy your $$x_{ij}$$ variables, so that they now take value $$1$$ if point $$i$$ is in the cluster with center point $$j$$. For consistency you need to add the following constraints:

• You need $$C$$ centers, one for each cluster: $$\sum_{j} y_j = C$$
• If point $$i$$ is in cluster $$j$$, then variable $$y_j$$ needs to be activated: $$x_{ij} \le y_j \quad \forall i,j$$
• A point belongs to exactly one cluster: $$\sum_{j}x_{ij} = 1 \quad \forall i$$
• The size of each cluster is bounded: $$L_j \le \sum_{i}x_{ij} \le U_j\quad \forall j$$

You can now minimize the average distance from a point to the center of its cluster: $$\min\;\sum_{i}\sum_{j}d_{ij}x_{ij}$$ Or perhaps the maximum distance from a point to its center: $$\min\;\max_{i,j}\;\{d_{ij}x_{ij}\}$$

where $$d_{ij}$$ denotes the distance between points $$i$$ and $$j$$.

If you really don't want the centers to be one the points, then things get a little messier (linearity is sacrificed). You could start by plotting the points (and add the plot to your question), so that we can eyeball if it makes sense to restrict the centers to one of the points.

EDIT: since the distribution of points seems to be sparse, it definitely makes sense to consider centers which are not initial points. In this case, you can solve the following (non linear) problem. $$\min \; \sum_{i}\sum_{j} d_{ij}^2x_{ij}$$

subject to your constraints, plus the following ones:

• Define $$d_{ij}$$ as the Euclidean distance between point $$i$$ and center $$\mu_j$$: $$d_{ij}^2=(X_i - X_{\mu_j})^2+(Y_i - Y_{\mu_j})^2$$
• Define the coordinates of the centers $$\mu_{j}$$ as the weighted sum of all the points in cluster $$j$$: $$X_{\mu_j} = \frac{\sum_i X_i x_{ij}}{\sum_i x_{ij}} \\ Y_{\mu_j} = \frac{\sum_i Y_i x_{ij}}{\sum_i x_{ij}}$$

Again, this is not linear...

• thanks a lot. I will check this. But yes, I wanted artificial centers instead. May 5, 2021 at 9:43
• Thanks a lot for your supports. I have added an example image showing two different distribution of points. The left one is more irregular while the right one is kind of regular. May 5, 2021 at 9:54
• My system will have around 36 points with left-image like distribution. I need an efficient LP/ILP formulation. Looking forward to your support on this. May 5, 2021 at 10:07
• thanks a lot. Looks only the objective function is nonlinear. So, if linear approximation of Euclidean distance can be performed, the problem becomes ILP, right? May 5, 2021 at 10:36
• do I still need the constraints you devised for your previous solution? May 5, 2021 at 10:39

I had the following model lying around:

\begin{aligned} \min&\sum_{i,k}\color{darkred}d_{i,k}\\ & \color{darkred}d_{i,k} \ge \sum_c \left(\color{darkblue}p_{i,c}-\color{darkred}\mu_{k,c}\right)^2-\color{darkblue}M(1-\color{darkred}x_{i,k}) \\ & \sum_k\color{darkred}x_{i,k} = 1 && \forall i\\ & \color{darkred}n_k = \sum_i \color{darkred}x_{i,k} \\ & \color{darkred}n_k \in [\color{darkblue}\ell,\color{darkblue}u] \\ & \color{darkred}x_{i,k} \in \{0,1\}\\ & \color{darkred}\mu_{k,c} \in [\min_i\color{darkblue}p_{i,c},\max_i\color{darkblue}p_{i,c}] \\ & \color{darkred}d_{i,k} \ge 0 \\ & \color{darkblue}M=\max_{i,i'} \sum_c \left(\color{darkblue}p_{i,c}-\color{darkblue}p_{i',c}\right)^2\\ & c \in \{x,y\} \end{aligned}

with:

• $$\color{darkred}d_{i,k}$$: squared distance between point $$i$$ and cluster centroid $$k$$ or zero if point $$i$$ is not part of cluster $$k$$
• $$\color{darkblue}p_{i,c}$$: coordinates of point $$i$$
• $$\color{darkred}\mu_{k,c}$$: cluster centroid coordinates
• $$\color{darkred}x_{i,k}$$: assignment of points $$i$$ to clusters $$k$$
• $$\color{darkred}n_k$$: number of points in cluster $$k$$
• $$\color{darkblue}\ell,\color{darkblue}u$$ are your constraints on the cluster size

Notes:

• This is a convex MIQCP.
• The objective is the "within-cluster sum-of-squares".
• As with many quadratic models, performance is not that good. 50 points is no problem, but after that things may slow down.
• I am not at all sure if this the best formulation.
• The big-M constraint can be replaced by indicator constraints. But unfortunately, most solvers don't support quadratic indicator constraints. I don't know why this is the case.
• The coordinates of the cluster centroids are not explicitly calculated as the mean of the coordinates of the points inside the cluster. The minimization will automatically take care of that. The centroid is the best location for $$\color{darkred}\mu_{k,c}$$. This can be proven by looking at the first-order conditions (after fixing the assignments). It is quite nice that we can drop completely from our model the nonlinear constraints that make up the centroid calculation. We are super lucky here.
• The model assumes the number of clusters is given. A standard way to find the "optimal" number of clusters is to solve with $$K=1,2,3,..$$ clusters and see when adding more clusters does not really improve the objective anymore. When plotting this we often see some kink in the graph that makes visually selecting this point not too difficult.
• In a sense, this is a multi-objective problem: minimize within-cluster SS and minimize the number of clusters.
• This model should always be feasible.
• The objective can be interpreted as $$\min \sum_{i,k,c}\left(\color{darkblue}p_{i,c}-\color{darkred}\mu_{k,c}\right)^2\cdot \color{darkred}x_{i,k}$$ Unfortunately, this makes the model a non-convex MINLP.
• This model is essentially k-means clustering. Of course, there are many alternative clustering approaches. Some obvious variants of this model are:
1. using L1 distances. That makes the model a linear MIP so easier to solve.
2. using distances instead of squared distances. That makes the model an MISOCP (after some reformulations) and even slower.
• There is symmetry with respect to how the clusters are numbered. Adding an ordering constraint, say $$\color{darkred}\mu_{k,x}\le \color{darkred}\mu_{k+1,x}$$ may help. It can reduce the solution time (sometimes substantially).
• thanks a lot for your answer. I will definitely try to run it ASAP. However, would you please comment on the usage of L1 distance instead of L2-norm. May 5, 2021 at 22:20
• The first constraint becomes: $\color{darkred}d_{i,k} \ge \sum_c \left|\color{darkblue}p_{i,c}-\color{darkred}\mu_{k,c}\right|-\color{darkblue}M_2(1-\color{darkred}x_{i,k})$. The absolute value can be reformulated using standard techniques (variable splitting, bounding). Thigs will stay convex so no additional binary variables needed. The big-M constant should change a bit. We can use indicator constraints to prevent the big-M constant. May 6, 2021 at 22:11