Clustering points based on a distance matrix

Question

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?

Cross-posted:

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',
                                  linkage='average')
agg_results = agg.fit_predict(TT_np_matrix) # The solution is an array

Edit:

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).

Answer 1

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.

Answer 2

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.

Answer 3

If you are interested in a model-and-run approach, maybe you can have a look at LocalSolver. An example model for basic k-means is given here. LocalSolver finds near-optimal solutions in seconds until 10,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 LocalSolver allows delivering quality solutions quickly while scaling to large-scale instances.

Below is written the model using the LocalSolver modeling language, namely LSP. You can also find it written for Python, Java, C#, or C++ at the link mentioned above.

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 LocalSolver is a commercial solver. Nevertheless, academic and trial licenses are available for free.