Although I asked this question on stackoverflow to possibly reach a broader audience, I wonder your inputs about this problem. Without giving much research into this, I thought p-center problem ($x_{ij}=1$ if point is clustered around point $j$; 0, otherwise) can be useful to construct a simple model to partition the points into equally-sized clusters. With a proper constraint, I could also let $x_{ij}\in\mathbb{R}_{\geq 0}$, I guess. I am just afraid that the problem size might be too large for the exact solution, hence I tried to find the cure in heuristic/algorithmic solutions world. How would you approach this problem?


I am trying to cluster a dataset with 3,663 points. I have (x,y) coordinates and point-to-point weights available. Since I would like to use the weights instead of Euclidean distances, I implemented the AgglomerativeClustering with Python's sklearn thanks to Vivek Kumar's answer here. I created 20 clusters and had the following number of points in each cluster:

Num_points_in_cluster = [673,54,111,88,312,464,332,366,141,128,33,86,426,39,45,1,63,105,2,194]

As seen, some clusters have only 1 or 2 point(s). I would instead like to have number of points in clusters almost uniformly distributed. Is there a way to play around this within AgglomarativeClustering or is there any other method (preferably in sklearn world) that meets the criterion? I think, an alternative solution could be solving a p-center-like problem, but the time performance of off-the-shelf linear solvers can be a bottleneck as I also would like a quick solution.

For code-checking purposes, here is my code, where TT_np_matrix is my weight matrix in numpy matrix type, e.g., TT_np_matrix = array([[2,3,4],[5,2,3],[4,5,8]]):

agg = AgglomerativeClustering(n_clusters=20, affinity='precomputed',
agg_results = agg.fit_predict(TT_np_matrix) # The solution is an array


Below is my (though I say "my," it most likely exists somewhere readily) formulation to solve the problem with binary/linear programming.

\begin{align} \min&\quad \sum_{i,j\in I}c_{ij}x_{ij}&\\ \text{s.t.}&\quad \sum_{j\in I} y_{j} =k & \\ &\quad \sum_{i\in I}x_{ij}\geq \lfloor |I|/k\rfloor y_j &\forall j\in I\\ &\quad \sum_{j\in I}x_{ij}=1 &\forall i\in I\\ &\quad x_{ij}\leq y_j &\forall i,j\in I\\ &\quad x,y\in\{0,1\}, \end{align}

where $x_{ij}=1$ if point $i\in I$ is centered around point $j\in I$; $0$, otherwise, and $y_j=1$ if point $j\in I$ is selected as a center; $0$, otherwise. Here, $k$ defines the number of clusters, and $I$ denotes the set of points to be clustered/partitioned.

Although I write this, I do not want to implement it as I believe it will be computationally more expensive (regardless of exactness) than an off-the-shelf so-called machine learning clustering algorithm. I am instead seeking an open-source implementable function or its example coded in Python (too much constraint, sorry).

  • $\begingroup$ The model you added in your edit is infeasible. The constraint that requires that every cluster has size at least $floor(|I|/k)$ is incorrect: this constraint should only hold for those $j\in I$ for which $y_j=1$. In any feasible solution, exactly $|I|$ $x_{ij}$ variables will have the value $1$, the others will have value $0$. In your second constraint you require that $floor(|I|/k)\times |I|$ variables take the value $1$. $\endgroup$ Commented Feb 23, 2021 at 18:52
  • $\begingroup$ I think, multiplying the r.h.s. with $y_j$ fixed the problem there. Do you see any other issues leading to infeasibility? $\endgroup$
    – tcokyasar
    Commented Feb 23, 2021 at 20:24
  • $\begingroup$ yes this looks fine. Note that your model requires that all clusters have the exact same size. The version of the model I posted below is slightly more flexible in the sense that you can control the size difference between clusters. If you indeed desire that all clusters have the same size, then I would recommend using your model. $\endgroup$ Commented Feb 24, 2021 at 8:36
  • $\begingroup$ Well, I think, adjusting the rhs will create a similar flexibility. Since we don’t know how much big M (though I choose a strict bound) impacts the time-to-solution, I will implement both and report the time results once I can break through my laziness barrier :) $\endgroup$
    – tcokyasar
    Commented Feb 24, 2021 at 10:19

3 Answers 3


To solve this using a p-center formulation, you could use this base model:

\begin{align} P: \min&\quad \sum_{i,j\in I}c_{ij}x_{ij}&\\ \text{s.t.}&\quad \sum_{j\in I} x_{ij} =1 & \forall i\in I\\ &\quad \sum_{i\in I}x_{ii}\leq maxClusters &\\ &\quad x_{ij}\leq x_{jj} &\forall i,j\in I\\ &\quad u\geq \sum_{i\in I}x_{ij} &\forall j\in I\\ &\quad l\leq \sum_{i\in I}x_{ij}+M(1-x_{jj}) &\forall j\in I\\ &\quad x,u,l\geq 0 &\\ &\quad x_{ii} \;binary & \forall i\in I \end{align} where $I$ is the set of items, $x_{ij}$ are the assignment variables and variables $u$, $l$ record the cardinality of respectively the largest and smallest cluster. You can then either add an additional constraint like: \begin{align} &\quad u-l \leq maxClusterSizeDifference & \end{align} to limit the maximum tolerated difference in cluster size, or, you could add an additional weighted objective term to the objective: \begin{align} &\quad w(u-l) & \end{align} ,where $w$ is some weight constant. Given that the above model is simple and quick to implement, you might give its implementation a try and see whether a MIP model finds a good enough solution within a very small amount of time. (try setting the MIP emphasis parameter of the solver to focus on finding good feasible solutions fast).

The above model would serve as a good benchmark. If speed is really of the essence, I would implement some heuristic. E.g. first partition all your nodes $I$ into two clusters (say $I_1$ and $I_2$) using the above model. Then recursively, reapply this splitting procedure for each of the new clusters. Continue until the desired number of clusters is reached.

Finally, you can use some balanced k-means clustering algorithm.

  • $\begingroup$ Joris, thanks for your answer. After reading your model, I also wanted to write down what I was thinking. So, I edited my question. Assuming my model is also correct, I believe the complexity of yours (in terms of the number of variables) is lower ($I^2+2$ versus $I^2+I$) given $|I|>2$. I also think both models can be treated with continuous variables (for $x$ in yours and both $x$ and $y$ in mine). But, your model involves a big $M$ that can turn things a bit less efficient. It looks like I will give these a try and see how they perform on my instance. $\endgroup$
    – tcokyasar
    Commented Feb 23, 2021 at 11:14
  • $\begingroup$ @tcokyasar You can't make all variables continuous: this will lead to fractional solutions. In my model the $x_{ii}$ variables need to be binary variables; the remaining ones can be continuous. $\endgroup$ Commented Feb 23, 2021 at 18:56
  • $\begingroup$ You can get rid of the big-$M$ in the lower bound constraints by rewriting as $lx_{jj} \leq \sum_{i \in I} x_{ij} $ $\endgroup$
    – Sune
    Commented Feb 23, 2021 at 21:16
  • $\begingroup$ @Sune Then you end up with $lx_{jj}$ variables, which you then still need to tie to the $l$ variable: $l\leq lx_{jj} + M(1-x_{jj})$ for all $j\in I$. So this didn't really get rid of the big-M. $\endgroup$ Commented Feb 23, 2021 at 21:54
  • $\begingroup$ Ah, you are right, I was too quick to comment. I saw $l$ as a fixed lower bound on the cluster size. $\endgroup$
    – Sune
    Commented Feb 24, 2021 at 6:38

I'm not familiar with agglomarative clustering, but in general terms you are dealing with a bicriterion optimization problem. One criterion has to do with the "affinities" of points in a cluster. You may want to maximize (?) something to do with the pairwise weights. That might be the sum of the weights of all pairs sharing a cluster, or the minimum weight of any pair sharing a cluster, or something else. The second criterion has to do with cluster size. You may want to minimize the range of cluster sizes, or maximize the size of the smallest cluster, or minimize the standard deviation of the cluster sizes, or something else.

So your first step is to figure out exactly what your two criteria are. You second step is then to figure out how to manage the competition between them. You can optimize a weighted sum, or optimize one subject to a limit on how bad the other can get, or use some other multiobjective scheme (for instance, goal programming), or just settle for any Pareto efficient solution. Depending on your objectives, you may also need to constrain the number of clusters, to avoid getting each point as its own cluster or one cluster containing all the points. (It's not guaranteed you need to constrain the number of clusters or the size of a cluster.)

Once you have nailed down the two objectives and how to merge them, you can probably define a mixed integer programming model, or look for a heuristic or metaheuristic to get a "good" solution. The heuristic might be a clustering method or might just be a general purpose (meta)heuristic.


If you are interested in a model-and-run approach, maybe you can have a look at Hexaly. An example model for basic k-means is given here. Hexaly finds near-optimal solutions in seconds, with up to 100,000 points to be clustered.

Instead of the classical Boolean modeling approach, using so-called set variables makes the model natural to read and compact in memory. Then, the relaxation-free heuristics embedded inside Hexaly allow delivering quality solutions quickly while scaling to large-scale instances.

The model is written below using the Hexaly Modeler. You can also find it written for Python, Java, C#, or C++ at the above link.

function model() {
  // Set decisions: clusters[c] represents the points in cluster c
  clusters[1..k] <- set(nbObservations);

  // Each point must be in one cluster and one cluster only
  constraint partition[c in 1..k](clusters[c]);

  // Compute variances
  for [c in 1..k] {
    local cluster <- clusters[c];
    local size <- count(cluster);

    // Compute the centroid of the cluster
    centroid[d in 0..nbDimensions-1] <- size == 0 ? 0 :
            sum(cluster, o => coordinates[o][d]) / size;

    // Compute the variance of the cluster
    squares[d in 0..nbDimensions-1] <- sum(cluster,
        o => pow(coordinates[o][d] - centroid[d], 2));
    variances[c] <- sum[d in 0..nbDimensions-1](squares[d]);

  // Minimize the total variance
  obj <- sum[c in 1..k](variances[c]);
  minimize obj;

Disclaimer: Please note that Hexaly is a commercial solver. Nevertheless, academic and trial licenses are available for free.

  • $\begingroup$ If I am reading the code correctly, I see that LocalServer uses coordinates rather than a weight/distance matrix. Can you please clarify this? $\endgroup$
    – tcokyasar
    Commented Feb 23, 2021 at 10:31
  • $\begingroup$ @tcokyasar If centroids must be points of your dataset, then the Facility Location model will be closer to your need: localsolver.com/docs/last/exampletour/facilitylocation.html. In this case, you can use your "weights" to evaluate the variance. If centroids can be arbitrary points, your weights cannot be used directly, indeed. Maybe you know the function behind the weights, to be used instead of Euclidian distances. With LocalSolver, you plug any external function into your optimization model, whatever its nature (nonlinear, nonconvex, or even possibly noisy, then nonsmooth). $\endgroup$
    – Hexaly
    Commented Feb 23, 2021 at 14:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.