Neural Capacitated Clustering

1 Introduction

In recent years much progress has been achieved in applying deep learning methods to solve classical clustering problems. Due to its ability to leverage prior knowledge and information to guide the partitioning of the data, constrained clustering in particular has recently gained increasing traction. It is often used to incorporate existing domain knowledge in the form of pairwise constraints expressed in terms of must-link and cannot-link relations [Wagstaff et al., 2001]. However, another type of constraints has been largely ignored so far: cluster level constraints. This type of constraint can for example restrict each assignment group in terms of the total sum of weights which are associated with its members. The simplest case of such a constraint is the maximum size of the cluster, where each point exhibits a weight of one. In the more general case, weights and cluster capacities are real valued and can model a plenitude of practical applications.

Machine learning approaches are particularly well suited for cases in which many similar problems, i.e. problems from the same distribution, have to be solved. In general, most capacitated mobile facility location problems (CMFLP) [Raghavan et al., 2019] represent this setting when treating every relocation as a new problem. This is e.g. the case when considering to plan the location of disaster relief services or mobile vaccination centers for several days, where the relocation cost can be considered to be zero since the team has to return to restock at the end of the day. Other applications are for example the planning of the layout of factory floors which change for different projects or logistics problems like staff collection and dispatching [Negreiros et al., 2022]. Thus, we can utilize supervised learning to learn from existing data how to solve new unseen instances.

The corresponding formulation gives rise to well-known problems from combinatorial optimization, the capacitated p-median problem (CPMP) [Ross and Soland, 1977] where each center has to be an existing point of the data and the capacitated centered clustering problem (CCCP) [Negreiros and Palhano, 2006] where cluster centers correspond to the geometric center of their members. The general objective is to select a number K of cluster centers and find an assignment of points such that the total distance between the points and their corresponding centers is minimized while respecting the cluster capacity. Both problems are known to be NP-hard and have been extensively studied [Negreiros and Palhano, 2006].

1.1 Contributions

• We propose the first approach to solve general Capacitated Clustering Problems based on deep learning.
• Our problem formulation includes well-known problem variants like the CPMP and CCCP as well as more simple constraints on the cluster size.
• The presented approach achieves competitive performance on several artificial and real world datasets, compared to methods based on mathematical solvers while reducing run time by up to one order of magnitude. 1
• We present a cluster-first-route-second extension of our method as effective construction heuristic for the CVRP.

1 Additional ablation and generalization studies can be found in the supplementary material provided with the extended paper at https://arxiv.org/abs/2302.05134

2 Background

A typical clustering task is concerned with the grouping of elements in the given data and is normally done in an unsupervised fashion. This grouping can be achieved in different ways where we usually distinguish between partitioning and hierarchical approaches [Jain et al., 1999]. In this work we are mainly concerned with partitioning methods, i.e. methods that partition the data into different disjoint sub-sets without any hierarchical structure. Although clustering methods can be applied to many varying data modalities like user profiles or documents, in this work we consider the specific case of spatial clustering [Grubesic et al., 2014], that normally assumes points to be located in a metric space of dimension D, a setting often encountered in practical applications like the facility location problem [Ross and Soland, 1977].

2.1 Capacitated Clustering Problems (CCPs)

Let there be a set of n points N = {1, 2, . . . , n} with corresponding feature vectors xi ∈ RD of coordinates and a respective weight qi ∈ R associated with each point i ∈ N. Further, we assume that we can compute a distance measure d(xi, xj) for all pairs of points i, j ∈ N. Then we are concerned with finding a set of K capacitated disjoint clusters ck ∈ C, ck ⊂ N, ck ∩ cl = ∅ ∀k, l ∈ {1, . . . ,K} with capacities Qk > 0. The assignment of points to these clusters is given by the set of binary decision variables yik, which are 1 if point i is a member of cluster k and 0 otherwise.

Capacitated p-Median Problem

For the CPMP the set of possible cluster medoids is given by the set of all data points N and the objective is to minimize the distance (or dissimilarity) d(xi, xmk ) between all i and their cluster medoid mk:

s.t.

where each point is assigned to only one cluster (2), all clusters respect the capacity constraint (3), medoids are selected minimizing the dissimilarity of cluster ck (4) and y being a binary decision variable (5).

Capacitated Centered Clustering Problem

In the CCCP formulation, instead of medoids selected among the data points, centroids µk are considered, which correspond to the geometric center of the points assigned to each cluster ck, replacing (4) with

which in the case of the Euclidean space for spatial clustering considered in this paper has a closed form formulation:

with |ck| as cardinality of cluster ck. This leads to the new minimization objective of

3 RelatedWork

ClusteringAlgorithms Traditional partitioningmethods to solve clustering problems, like the well-known k-means algorithm [MacQueen, 1967] have been researched for more than half a century. Meanwhile, there exists a plethora of different methods including Gaussian Mixture Models [McLachlan and Basford, 1988], density based models like DBSCAN [Ester et al., 1996] and graph theoretic approaches [Ozawa, 1985]. Many of these algorithms have been extended to solve other problem variants like the CCP. In particular [Mulvey and Beck, 1984] introduce a k-medoids algorithm utilizing a regret heuristic for the assignment step combined with additional re-locations during an integrated local search procedure while [Geetha et al., 2009] propose an adapted version of kmeans which instead uses a priority heuristic to assign points to capacitated clusters.

Meta-Heuristics Apart from direct clustering approaches there are also methods from the operations research community which tackle CCPs or similar formulations like the facility location problem. Different algorithms were proposed modeling and solving the CCP as General Assignment Problem (GAP) [Ross and Soland, 1977], via simulated annealing and tabu search [Osman and Christofides, 1994], with genetic algorithms [Lorena and Furtado, 2001] or using a scatter search heuristic [Scheuerer and Wendolsky, 2006].

Math-Heuristics In contrast to meta-heuristics, mathheuristics combine heuristic methods with powerful mathematical programming solvers like Gurobi [Gurobi Optimization, LLC, 2023], which are able to solve small scale instance to optimality and have shown superior performance to traditional meta-heuristics in recent studies. [Stefanello et al., 2015] combine the mathematical solution of the CPMP with a heuristic post-optimization routine in case no optimality was achieved. The math-heuristic proposed in [Gn¨agi and Baumann, 2021] comprises two phases: First, a global optimization phase is executed. This phase alternates between an assignment step, which solves a special case of the GAP as a binary linear program (BLP) for fixed medoids, and a median update step, selecting new medoids mk minimizing the total distance to all cluster members under the current assignment. This is followed by a local optimization phase relocating points by solving a second BLP for the sub-set of clusters which comprises the largest unused capacity. Finally, the PACK algorithm introduced in [L¨ahderanta et al., 2021] employs a block coordinate descent similar to the method of [Gn¨agi and Baumann, 2021] where the assignment step is solved with Gurobi and the step updating the centers is computed following eq. 7 according to the current assignment.

Deep Clustering Since the dawn of deep learning, an increasing number of approaches in related fields is employing deep neural networks. Most approaches in the clustering area are mainly concerned with learning better representations for downstream clustering algorithms, e.g. by employing auto-encoders to different data modalities [Tian et al., 2014; Xie et al., 2016; Guo et al., 2017; Yang et al., 2017], often trained with enhanced objective functions, which, apart from the representation loss, also include a component approximating the clustering objective and additional regularization to prevent the embedding space from collapsing. A comprehensive survey on the latest methods is given in [Ren et al., 2022]. More recently, deep approaches for constrained clustering have been proposed: [Genevay et al., 2019] reformulate the clustering problem in terms of optimal transport to enforce constraints on the size of the clusters. [Zhang et al., 2021] present a framework describing different loss components to include pairwise, triplet, cardinality and instance level constraints into auto-encoder based deep embedded clustering approaches. Finally, [Manduchi et al., 2021] propose a new deep conditional Gaussian Mixture Model (GMM), which can include pairwise and instance level constraints. Usually, the described deep approaches are evaluated on very large, high dimensional datasets like MNIST [LeCun and Cortes, 2010] or Reuters [Xie et al., 2016], on which classical algorithms are not competitive. This is in strong contrast to spatial clustering with additional capacity constraints, for which we propose the first deep learning based method.

4 ProposedMethod

4.1 Capacitated K-Means

The capacitated k-means method proposed by [Geetha et al., 2009] changes the assignment step in Lloyd's algorithm [Lloyd, 1982], which is usually used to compute kmeans clustering. To adapt the procedure to the CCP, the authors first select the K points with the highest weights q as initial centers, instead of selecting them randomly. Moreover, they introduce priorities ωik for each point i by dividing its weight qi by its distance to the center of cluster k:

Then the list of priorities is sorted and nodes are sequentially assigned to the centers according to their priority. The idea is to first assign points with large weights to close centers and only then points with smaller weight, which can be more easily assigned to other clusters. Next, the centroids are recomputed via the arithmetic mean of the group members (eq. 7). The corresponding pseudo code is given in Alg. 1.

While this heuristic works, it can easily lead to sub-optimal allocations and situations in which no feasible solution can be found, e.g. in cases where many nodes with high weight are located very far from the cluster centers. To solve these problems we propose several modifications to the algorithm.

4.2 Neural Scoring Functions

The first proposed adaption is to learn a neural scoring function fθ with parameters θ, which predicts the probability of

input : K, n, coordinates x, weights q, cluster capacity Q, convergence condition δ output : binary assignment matrix Y 1 M ←inittopk weights (x, q,K) 2 while not δ(x,M, Y ) do 3 Y←allzero (n,K) // reset assignment 4 Q←repeat (Q,K) // reset capacities 5 foreach i ∈ N do 6 compute priorities for all clusters (eq. 9) 7 sort priorities, insert into queue S 8 while S not empty do 9 get next i, k from S 10 if i unassigned and Qk ≥ qi then 11 Yik ← 1 // cluster assignment 12 Qk ←Qk - qi // update capacity 13 foreach k ∈ {1, . . . , K} do 14 update centroids via eq. 7 15 return Y

each node i to belong to cluster k:

For that purpose we first create a graph representation G = (V, E) of the points by connecting each point with its K nearest neighbors, producing edges eij ∈ E with edge weights d(xi, xj). Nodes vi ∈ V are created by concatenating the respective coordinates and weights [xi; qi]. This graph allows us to define a structure on which the relative spatial information of the different points can be efficiently propagated. We encode G with the Graph Neural Network (GNN) introduced in [Morris et al., 2019], which is able to directly work with edge weights by employing the graph operator defined as

where h(l-1) i ∈ R1×demb represents the embedding of node i at the previous layer l - 1, H(i) is the 1-hop graph neighborhood of node i, eji is the directed edge connecting nodes j and i, MLP1 and MLP2 are Multi-Layer Perceptrons MLP : Rdemb → Rdemb and σ() is a suitable activation function. Furthermore, we add residual connections and regularization to each layer. In our case we choose GELU [Hendrycks and Gimpel, 2016] and layer normalization [Ba et al., 2016] which outperformed ReLU and BatchNorm in preliminary experiments. The input layer projects the node features vi = [xi; qi] ∈ RD+1 to the embedding dimension demb using a feed forward layer, which then is followed by L GNN layers of the form given in eq. 11. In order to create embeddings hµk for the centers µk ∈ M we find the node j closest to µk (corresponding to the cluster medoid mk) and select its embedding h(L) j as hk. This embedding is concatenated with a globally pooled graph embedding hG ∈ Rdemb :

with MLPG : R2demb → Rdemb . Next, in order to model the interdependence of different centers, their resulting vectors are fed into a self attention layer (SA) [Vaswani et al., 2017] followed by another MLPµ : R2demb → Rdemb :

Since we feed a batch of all points N and all currently available centers µ to our model, they are encoded as a batch of context vectors hµk , one for each cluster k ∈ {1, . . . ,K}. The SA module then models the importance of all cluster center configurations hµ for each cluster k which allows the model to effectively decide about the cluster assignment probability of every node i to each cluster ck conditionally on the full information encoded in the latent context embeddings hµ. Finally, we do conditional decoding by concatenating every center embedding with each node and applying a final stack of (element-wise) MLPs. The architecture of our neural scoring function is shown in Figure 1.

Training the model We create the required training data and labels by running the math-heuristic solver of [Gn¨agi and Baumann, 2021] on some generated datasets to create a good (although not necessarily optimal) partitioning. Then we do supervised training using binary cross entropy 2 (BCE) with pairwise prediction of the assignment of nodes i to clusters k.

4.3 Neural Capacitated Clustering

To fully leverage our score estimator we propose several adaptions and improvements to the original capacitated kmeans algorithm (Alg. 1). The new method, which we dub Neural Capacitated Clustering (NCC) is described in Alg. 2.

Order of Assignment

Instead of sorting all center-node pairs by their priority and then sequentially assigning them according to that list, we fix an order given by permutation π for the centers and cycle through each of them, assigning one node at a time. Since the output of fθ is the log probability of point i belonging to cluster k and its magnitude does not directly inform the order of assignments of different nodes i and j in the iterative cluster procedure, we found it helpful to scale the output of fθ by the heuristic weights introduced in eq. 9. Thus, we assign that node i to cluster k, which has the highest scaled conditional priority and still can be accommodated considering the remaining capacity Qk. In case there remain any

unassigned points j at the end of an iteration, which cannot be assigned to any cluster since qj > Qk ∀k ∈ K, we assign them to a dummy cluster K + 1 located at the origin of the coordinate system. We observe in our experiments that already after a small number of iterations usually no nodes are assigned to the dummy cluster anymore, meaning a feasible allocation has been established. Moreover, since the neighborhood graph G does not change between iterations, we can speed up the calculation of priorities by pre-computing and buffering the node embeddings hi and graph embedding hG in the first iteration.

Re-Prioritization of Last Assignments

This is motivated by the observation that the last few assignments are the most difficult, since they have to cope with highly constrained center capacities. Thus, relying on the predefined cyclic order of the centers (which until this point has ensured that approx. the same number of nodes was assigned to each cluster) can lead to sub-optimal assignments in case some clusters have many nodes with very large or very small weights. To circumvent this problem we propose two different assignment strategies:

1. Greedy: We treat the maximum (unscaled) priority over all clusters as an absolute priority ¯ωi for all remaining unassigned points i:

Then the points are ordered by that priority and sequentially assigned to the closest cluster which can still accommodate them.

1. Sampling: We normalize the absolute priorities ¯ωi of all remaining unassigned points i via the softmax 3 function and treat them as probabilities according to which they are sequentially sampled and assigned to the closest cluster which can still accommodate them. This procedure can be further improved by sampling several assignment rollouts and selecting the configuration, which leads to the smallest resulting inertia.

The fraction α of final nodes for which the re-prioritization is applied we treat as a hyperparameter.

Algorithm 2: Neural Capacitated Clustering (NCC)

Weight-Adapted Kmeans++ Initialization

As found in the study of [Celebi et al., 2013] standard methods for the selection of seed points for centers during the initialization of k-means algorithms perform quite poorly. This is why the k-means++ [Arthur and Vassilvitskii, 2006] initialization routine was developed, which aims to maximally spread out the cluster centers over the data domain, by sampling a first center uniformly from the data and then sequentially sampling the next center from the remaining data points with a probability equal to the normalized squared distance to the closest already existing center. We propose a small modification to the k-means++ procedure (called ckm++), that includes the weight information into the sampling procedure by simply multiplying the squared distance to the closest existing cluster center by the weight of the data point to sample.

5 Experiments

We implement our model and the simple baselines in PyTorch [Paszke et al., 2019] version 1.11 and use Gurobi version 9.1.2 for all methods that require it. All experiments are run on a i7-7700K CPU (4.20GHz). We use L = 4 GNN layers, an embedding dimension of demb = 256 and a dimension of dh = 256 for all hidden layers. More details on training our neural scoring function we report in the supplementary. 4

5.1 Capacitated Clustering

For the experiments we use the CCCP formulation of the CCP (eq. 8) which considers centroids instead of medians. While there are several possible ways to select a useful number K of clusters, like the Elbow method [Yuan and Yang, 2019], here we adopt a practical approach consisting of solving the problem with the random assignment baseline method for a number of seeds and choosing the minimal resulting number of clusters as K. For the n = 200 datasets we run all methods for 3 different seeds and report the mean cost with standard deviation and average run times. Since Gurobi requires a run time to be specified because it otherwise can take arbitrarily long for the computation to complete, we set reasonable total run times of 3min for n = 200 and 15min for the original sizes. If Gurobi times out, we return the last feasible assignment if available. Otherwise, we report the result for that instance as infeasible and set its cost to the average cost of the rnd-NN baseline. The validation accuracy of our neural scoring function on the different training sets is shown in Figure 2. We evaluate our method in a greedy and a sampling configuration, which we tune on a separate validation set: (g20-1) stands for 1 greedy rollout for a fraction of α = 0.2 and (s-25-32) for 32 samples for α = 0.25.

4 We open source our code at https://github.com/jokofa/NCC

Datasets We performexperiments on artificial data and two real world datasets. The artificial data with instances of size n = 200 we generate based on a GMM. As real world datasets for capacitated spatial clustering we select the wellknown Shanghai Telecom (ST) dataset [Wang et al., 2019] which contains the locations and user sessions for base stations in the Shanghai region. In order to use it for our CCP task we aggregate the user session lengths per base station as its corresponding weight and remove base stations with only one user or less than 5min of usage in an interval of 15 days as well as outliers far from the city center, leading to a remaining number of n = 2372 stations. We set the required number of centers to K = 40. The second dataset we assemble by matching the internet access sessions in the call record data of the Telecom Italia Milan (TIM) dataset [Barlacchi et al., 2015] with the Milan cell-tower grid retrieved from OpenCelliD [OCID, 2021]. After pre-processing it contains n = 2020 points to be assigned to K = 25 centers. We normalize all weights according to K with a maximum capacity normalization factor of 1.1. The experiments on the real world data are performed in two different settings: The first setting simply executes all methods on the full dataset, while the second setting sub-samples the data in a random local grid to produce 100 test instances of size n = 200, with weights multiplied by a factor drawn uniformly from the interval [1.5, 4.0) for ST and [2.0, 5.0) for TIM to produce more variation in the required K. The exact pre-processing steps and subsampling procedure are described in the supplementary.

Baseline Methods

• random: sequentially assigns random labels to points while respecting cluster capacities.
• rnd-NN: selects K random points as cluster centers and sequentially assigns nearest neighbors to these clusters, i.e. points ordered by increasing distance from the center, until no capacity is left.
• topk-NN: similar to random-NN, but instead selects the K points with the largest weight as cluster centers.
• CapKMeans: the capacitated k-means algorithm of [Geetha et al., 2009] with ckm++ initialization which outperformed the original topk weights initialization.
• PACK: the block coordinate descent math-heuristic introduced in [L¨ahderanta et al., 2021] (using Gurobi).
• GB21: the two phase math-heuristic proposed by [Gn¨agi and Baumann, 2021] also using the Gurobi solver.

Results We evaluate all baselines in terms of inertia, which is defined as the total squared distance of all points to their assigned center. On the generated data (Table 1) our method outperforms all other baselines in terms of inertia while being much faster than all other methods with comparable performance. Furthermore, for three runs with different random seeds our method achieves results with significantly smaller standard deviations and is able to solve all instances within the given time. Results for the sub-sampled datasets are reported in Table 2. On ST our method beats close to all methods in terms of inertia and is only slightly outperformed by GB21 while being 5× faster. On TIM we achieve similar performance compared to GB21 while being much faster and outperform all other approaches. In particular, we are more than one order of magnitude faster than the next best performing baseline PACK. Furthermore, our method again achieves very small standard deviation for greedy as well as sampling based assignments, showing that it reliably converges to good (local) optima. As expected the random baseline leads to very bad results. The naive baselines rnd-NN and topk-NN are better but still significantly worse than the more advanced methods, achieving inertia which are 3-7 times higher than that of the best method. Compared to CapKMeans NCC leads to improvements of 41%, 37% and 17% respectively for the generated, ST and TIM data while achieving even faster run times, showing the impact of our proposed model adaptations. The inertia and run times for the full datasets are directly reported in the headings of Figures 3 and 4, which display the cluster assignments for Milan and Shanghai, explicitly drawing the convex hull of each cluster for better visualization. While both math-heuristics are able to outperform NCC on these large-scale instances in terms of inertia, our method is close to one order of magnitude faster. Furthermore, the results show that our method is able to find useful cluster structures which especially for TIM are more homogeneous than those of GB21 and show vast improvements compared to CapKMeans.

5.2 Capacitated Vehicle Routing

To show the efficacy of our approach we extend it to a clusterfirst-route-second (C1R2) construction method for Capacitated Vehicle Routing Problems (CVRP). The CVRP is an extension of the traveling salesman problem (TSP) in which K capacitated vehicles have to serve the demand of n customers from a fixed depot node [Toth and Vigo, 2014].

AlgorithmModifications To adapt our method for the different problem we include an additional MLP in our scoring function, which encodes the depot node, and concatenate the depot embedding with hG and hk in eq. 13. In our algorithm we add the depot node to each cluster k during the center update (Algorithm 2, line 22) and add the distance from each node to the depot to the priority weights. After the clustering we use the fast TSP solver provided by VeRyPy [Rasku et al., 2019] to route the nodes in each assigned group.

Dataset To evaluate our algorithm we choose the benchmark dataset of [Uchoa et al., 2017] which consists of 100

instances of sizes between 100 and 1000 points sampled according to varying distributions (uniform, clustered, etc.) and with different depot positions and weight distributions. We split the benchmark into three sets of problems with size N1 (100 ≤ n < 250), N2 (250 ≤ n < 500) and N3 (500 ≤ n).

Baselines In our experiments we compare against several classical C1R2 approaches: First, the sweep algorithm of [Gillett and Miller, 1974], which starts a beam at a random point and adds nodes in turn by moving the beam around the depot. We restart the beam at each possible point and run it clock and counter-clock wise. Next, sweep+, which instead of routing nodes in the order in which they were passed by the beam, routes them by solving a TSP with Gurobi. The petal algorithm introduced in [Foster and Ryan, 1976] creates "petal" sets by running the sweep algorithm from different starting nodes and then solves a set covering problem with Gurobi to join them. Finally, for comparison (although not C1R2) the powerful auto-regressive neural construction method POMO of [Kwon et al., 2020] which is trained with deep reinforcement learning and uses additional instance augmentation techniques. It is evaluated either greedily (g) or with sampling (s) and a beam width of n (size of the instance).

Results As shown in Table 3, our extended approach performs very competitive on the benchmark, beating all C1R2 approaches from the classical literature and being close to POMO on the small and medium sized instances (N1 and N2) while significantly outperforming it on the large instances (N3). Moreover, our method achieves the smallest fleet size of all methods, very close to the optimal fleet size Koptimal.

6 Conclusion

This paper presents the first deep learning based approach for the CCP. In experiments on artificial and real world data our method NCC shows competitive performance and fast and robust inference. Moreover, we demonstrate its usefulness as constructive method for the CVRP, achieving promising results on the well-known Uchoa benchmark.

Ethical Statement

There are no ethical issues.