Solution strategy based on Gaussian mixture models and dispersion reduction for the capacitated centered clustering problem

Gaussian mixture models

Gaussian Mixture Models (GMMs) have been widely studied for data modeling within the field of pattern classification based on Bayes decision theory (Theodoridis & Koutroumbas, 2010). For more accurate modeling of multi-dimensional patterns a mixture of Gaussian distributions can be used. Each mixture component is defined by two main parameters: a mean vector and a covariance matrix (Forsyth, 2012; Theodoridis & Koutroumbas, 1979; Theodoridis & Koutroumbas, 2010). Within this context, “clusters” can be characterized by each Gaussian component (mixture) which can model sub-sets of the whole set of patterns with maximum likelihood.

There are important differences between clustering that can be obtained with GMMs and centroid-based clustering techniques (such as K-Means or K-Nearest Neighbors). The following can be mentioned:

• Over-fitting (e.g., the model “overreacts” to minor fluctuations in the training data for prediction purposes) can be avoided with Gaussian distribution-based clustering.

• Clusters defined with Gaussian distributions can have different shapes and sizes. By contrast, centroid-based clustering algorithms tend to find clusters of comparable size of (more or less) symmetrical shape (Mohammed, Badruddin & Mohammed, 2016).

• At each iteration, Gaussian distribution-based clustering performs, for a given point, a “soft-assignment” to a particular cluster (there is a degree of uncertainty regarding the assignment). The centroid-based clustering performs a hard-assignment (or direct assignment) where a given point is assigned to a particular cluster and there is no uncertainty.

Due to these differences, the GMM-based clustering was considered as an alternative to generate feasible solutions for the CCCP. In terms of the CCCP formulation described in Section “Introduction” a cluster can be modeled by a single Gaussian probability density function (PDF). Hence, the location “patterns” of a set of clients X can be modeled by a mixture of K Gaussian PDFs where each PDF models a single cluster. If the set contains N clients, then $X=\left[x_{\mathrm{i}=1},x_{\mathrm{i}=2},...,x_{\mathrm{i}=\mathrm{N}}\right]$ and the mixture can be expressed as: (7) $$p(X)=\sum _{\mathrm{k}=1}^{\mathrm{K}}{P}_{\mathrm{k}}p(X|k)=\sum _{\mathrm{k}=1}^{\mathrm{K}}{P}_{\mathrm{k}}\mathcal{N}(X|m_{\mathrm{k}},S_{\mathrm{k}}),$$where k = 1,…, K and |K| = p is the number of Gaussian PDFs, p(X | k) represents the probabilities of each Gaussian PDF describing the set of clients X (Theodoridis & Koutroumbas, 2010) and P_k is the weight associated to each Gaussian PDF (hence, ${\sum }_{\mathrm{k}=1}^{\mathrm{K}}{P}_{\mathrm{k}}$ = 1.0). Each Gaussian component can be expressed as: (8) $$p(X|k)=\mathcal{N}(X|m_{\mathrm{k}},S_{\mathrm{k}})={\displaystyle \frac{1}{\sqrt{2\pi |S_{\mathrm{k}}|}}{e}^{-{\displaystyle \frac{1}{2}{(\mathrm{X}-{\mathrm{m}}_{\mathrm{k}})}^{\mathrm{T}}{\mathrm{S}}_{\mathrm{k}}^{-1}(\mathrm{X}-{\mathrm{m}}_{\mathrm{k}})}}\mathrm{\forall }k\in K,}$$where m_k is the mean vector and S_k is the covariance matrix of the k-th Gaussian PDF or k-th cluster. Note that m_k and each element of X (i.e., any x_i) must have the same size or dimension which is defined by l (in this case, l = 2 because each point consists of a x-coordinate and a y-coordinate). Finally, S_k is a matrix of size l × l.

For clustering purposes, m_k can model the mean vector of a sub-set of points in X which is more likely to be described by the Gaussian PDF k as estimated by Eq. (8). If the points in this sub-set of X are clustered, then m_k can represent the “centroid” of the cluster k and S_k can model the size and shape of the cluster k (Bishop, 2006; Theodoridis & Koutroumbas, 2010). Observe that X − m_k defines a distance or difference between each point in X and the centroid (located at m_k) of each cluster k. Thus, the estimation of probabilities considers the distance between each point x_i and each cluster k.

The parameters of the Gaussian PDFs (P_k, m_k and S_k) that best model (describe) each sub-set of the whole pattern X can be estimated by the iterative algorithm of Expectation-Maximization (EM) (Bishop, 2006; Theodoridis & Koutroumbas, 2010).

The advantage of this Gaussian approach for clustering is that faster inference about the points x_i that belong to a specific cluster k may be obtained considering all points. In this context, it is important to mention that due to the probabilistic nature of the inference process, a single point x_i is not directly assigned to a specific cluster (as it is required by the CCCP) because all points have probabilities associated to all clusters (i.e., “soft-assignment”). Also, this approach does not consider the capacity of each cluster and the demand associated to each client point. Thus, restrictions about the quantity of points x_i associated to each cluster are not integrated in this algorithm. In order to consider these requirements and restrictions a modification on the EM algorithm was performed. This is described in the following sections.

Dispersion reduction in GMMs performance for the CCCP

Capacitated Centered Clustering Problem data, which consists of x and y coordinates, represents a particular challenge for GMMs. This is because the values of coordinates can affect the computation of the exponential element of Eq. (8). Also, dispersion of data may affect convergence of the EM algorithm for P_k, m_k and S_k. Thus, specific re-scaling or coding of CCCP data was required to reduce the effect of dispersion and coordinates’ values on the computations of GMMs.

In order to reduce dispersion of the CCCP data the compression algorithm presented in Fig. 1 was performed. It is important to mention that compression is only performed on the coordinates’ values, and not on the number of data points (thus, the size of the instance remains the same).

As presented, the original x–y coordinates were replaced by their unique indexes. This led to elimination of empty spaces between data points. Coding of the compressed data was performed within the interval $\left[0,1\right]$ as presented in Fig. 1. This coding made the compressed data more suitable for fast computation of Eq. (8).

An important issue regarding the use of GMMs for the CCCP is the restriction on the number of Gaussian components. In this regard, extensive research has been performed to determine the most suitable number of Gaussian components for different sets of data (McLachlan & Rathnayake, 2014). However, available research does not consider the capacity aspect of clustering which is the main feature of CCCP data. Nevertheless, a restriction on the number of Gaussian components must be considered because, as discussed in Fraley & Raftery (1998), the EM algorithm may not be suitable for cases with very large numbers of components.

Due to this situation a maximum number of 30 Gaussian components was considered. This restriction was based on the DONI database, which has some of the CCCP instances with the largest number of client points (Fernandes-Muritiba et al., 2012) (i.e., N = 13,000 client points with p = 30 facilities). A discussion on future extensions for the Gaussian approach and the EM algorithm to address larger number of components is presented in the final section.

Standard EM algorithm

Figure 2 presents the structure of the standard EM algorithm (Bishop, 2006; Theodoridis & Koutroumbas, 1979; Theodoridis & Koutroumbas, 2010) considering the variables defined by Eqs. (7) and (8). In this case, X represents the array of compressed and coded x − y coordinates of all points x_i (the array structures presented in Figs. 2 and 3 follow the EM formulation described in Theodoridis & Koutroumbas (1979)).

As presented in Fig. 2, the locations of the N clients are stored into the array X which consists of a matrix of dimension l × N where: (a) N is the number of client points, and (b) l = 2 as each column vector of X consists of the values cx_i and cy_i that identify the compressed and coded x − y coordinates of the client point x_i where i = 1,…, N.

The EM algorithm starts with initial values for m_k, S_k and P_k. Values for m_k and S_k were randomly generated as follows: (9) $$m_{\mathrm{k}}={\left[\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{m}(\mathrm{c}{\mathrm{x}}_{\mathrm{m}\mathrm{i}\mathrm{n}},\mathrm{c}{\mathrm{x}}_{\mathrm{m}\mathrm{a}\mathrm{x}}),\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{m}(\mathrm{c}{\mathrm{y}}_{\mathrm{m}\mathrm{i}\mathrm{n}},\mathrm{c}{\mathrm{y}}_{\mathrm{m}\mathrm{a}\mathrm{x}})\right]}^T\mathrm{\forall }k\in K,$$ (10) $$S_{\mathrm{k}}=\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{m}(0.0,0.1)\times I_{\mathrm{l}}\mathrm{\forall }k\in K,$$where (cx_min, cx_max) and (cy_min, cy_max) are the minimum and maximum values throughout all compressed and coded x and y coordinates respectively, and I_l is the identity matrix of size l × l.

For P_k a lower bound for K was obtained by considering the total demand of the points x_i and the capacity of the clusters C_k. Because all clusters have the same capacity, C_k = C. Then, K and P_k were obtained as follows: (11) $$K={\displaystyle \frac{\sum _{\mathrm{i}=1}^{\mathrm{N}}{d}_{\mathrm{i}}}{C},}$$ (12) $$P_{\mathrm{k}}={\displaystyle \frac{1}{|K|}\mathrm{\forall }k\in K.}$$

The stage of Expectation starts with these initial values for m_k, S_k and P_k. An initial computation of assignment or “responsibility” scores γ (z_ik) is performed to determine which x_i is more likely to be associated to a particular cluster (and thus, to belong to this cluster) with parameters m_k, S_k and P_k (Bishop, 2006). Observe that, as presented in Fig. 2 (Step 2), γ (z_ik) is computed by means of Eq. (8). These scores can lead to provide values for the decision variable y_ik of the CCCP objective function (see Eq. (1)). This process will be discussed in the following section.

Then, the stage of Maximization integrates the scores γ (z_ik) into the re-estimation of the cluster’s parameters which leads to $m_{\mathrm{k}}^{\mathrm{n}\mathrm{e}\mathrm{w}}$, $S_{\mathrm{k}}^{\mathrm{n}\mathrm{e}\mathrm{w}}$ and $P_{\mathrm{k}}^{\mathrm{n}\mathrm{e}\mathrm{w}}$. Convergence is achieved if the total error or absolute difference between the previous and new estimates is less than a given threshold e (in this case, e=0.5). If convergence is not achieved, then $m_{\mathrm{k}}\leftarrow m_{\mathrm{k}}^{\mathrm{n}\mathrm{e}\mathrm{w}}$, $S_{\mathrm{k}}\leftarrow S_{\mathrm{k}}^{\mathrm{n}\mathrm{e}\mathrm{w}}$ and $P_{\mathrm{k}}\leftarrow P_{\mathrm{k}}^{\mathrm{n}\mathrm{e}\mathrm{w}}$.

Modified EM algorithm with direct assignment and capacity restrictions

As discussed in Section “Gaussian Mixture Models” the advantage of the GMMs approach is that faster inference about the points x_i that belong to a specific cluster k can be obtained. This inference is performed based on the probabilities defined by Eq. (8) where the parameters of the Gaussian PDFs are estimated with the standard EM algorithm. However under this approach a single point x_i is not directly assigned to a specific cluster (as it is required by the CCCP) because all points have probabilities associated to all clusters (i.e., “soft-assignment”). Also, while the standard EM algorithm can determine the most likely sub-sets of client points in X to be covered by each cluster k, these sub-sets are not restricted by the capacity of the cluster and the demand of each assigned client point. Hence, the standard EM algorithm was modified in order to comply with the requirements and restrictions of the CCCP.

As previously mentioned, the score γ (z_ik) computed in Step 2 of the EM algorithm (see Fig. 2) represents the likelihood or responsibility of the cluster k on the generation of the point x_i (Bishop, 2006). As presented in Fig. 3 all scores are stored into a matrix γ of dimension K × N, where each column vector ${\left[\gamma (z_{\mathrm{i}1}),...,\gamma (z_{\mathrm{i}\mathrm{K}})\right]}^{\mathrm{T}}$ contains the assignment scores of the point x_i to all clusters (thus, ${\sum }_{\mathrm{k}=1}^{\mathrm{K}}\gamma (z_{\mathrm{i}\mathrm{k}})=1.0\mathrm{\forall }i\in N$). These scores represent the basis for defining the decision variable y_ik.

From Section “Introduction” it was explained that y_ik = 1 if the point x_i is assigned to cluster k and y_ik = 0 otherwise (each point x_i can be assigned to only one cluster). Based on the scores of γ (z_ik) it was determined that for each vector x_i the cluster assignment would be defined by the cluster k with maximum γ(z_ik) value. If two or more clusters share the same maximum likelihood, then one of them is randomly assigned. An example of this assignment process is presented in Fig. 3.

By determining the unique assignment of each point x_i to each cluster k at Step 2 of the EM algorithm (see Fig. 2), the number of points assigned to each cluster is obtained. This leads to determine the cumulative demand of the points assigned to each cluster. This information is stored into the vector: (13) $$\mathrm{D}\mathrm{e}\mathrm{m}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{s}=\left[D_1,D_2,...,D_{\mathrm{K}}\right],$$where D_k represents the cumulative demand of the points assigned to cluster k and it must satisfy D_k ≤ C_k. This vector is important to comply with the capacity restrictions because it was found that homogenization of the cumulative demands D_k contributes to this objective. Homogenization is achieved by minimizing the coefficient of variation between all cumulative demands: (14) $$min\mathrm{C}\mathrm{V}={\displaystyle \frac{\mathrm{S}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{r}\mathrm{d}\mathrm{D}\mathrm{e}\mathrm{v}\mathrm{i}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}(\mathrm{D}\mathrm{e}\mathrm{m}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{s})}{\mathrm{M}\mathrm{e}\mathrm{a}\mathrm{n}(\mathrm{D}\mathrm{e}\mathrm{m}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{s})}.}$$

The objective function defined by Eq. (14) is integrated within the evaluation step of the standard EM algorithm (see Step 4 of Fig. 2). This leads to the modified EM algorithm with capacity restrictions where convergence is based on two objectives:

• minimization of the error (less than a threshold of e = 0.5) between the estimates of m_k, P_k and S_k;

• minimization of the coefficient of variation of Demands.

An important issue regarding the minimization of CV is that there is the possibility that min CV may not lead to comply the condition D_k ≤ C_k for all clusters k. This was addressed by the following strategy: if min CV leads to a cluster k where the condition D_k ≤ C_k is not complied, then the farthest assigned points are re-assigned to other clusters with closer centroids and available capacity.

Figure 4 presents an example of the clustering achieved without and with the objective function defined by Eq. (14). As presented, the clusters are more accurately defined with the integration of Eq. (14). Also, clusters comply with the capacity restrictions of the CCCP.

It is important to remember that the proposed DRG meta-heuristic consists of the modified EM algorithm which provides the cluster assignation for each point x_i considering compressed and coded CCCP data. Thus, Fig. 4 presents the decoded and uncompressed points x_i (i.e., original cx_i and cy_i coordinates) assigned to each cluster based on the assignments of the DRG meta-heuristic with the modified EM algorithm on compressed and coded data.

Further improvement of the assignments and the centroids at m_k are achieved by a Greedy algorithm that verifies that all points x_i are assigned to the closest cluster. This leads to additional exchange and insertion/deletion of complying client points x_i between clusters. These operations must comply with the capacity restrictions of each cluster. The assignments are updated if the insertion/deletion is valid. This leads to re-estimation of the locations of the m_k centroids.

DRG vs. VNS-SA-CS-TS-GA

Table 2 presents the best results of the DRG meta-heuristic for the considered instances. Information regarding the runs performed by each method to report the best result is also presented when available. Also, information regarding the programing language and the hardware used by the authors of the reviewed methods were also included.

As reported in the literature, CS and TS are the most competitive methods to solve the CCCP with a mean error gap of 0.78% and 0.61% respectively (and thus, were considered to update the benchmark solutions). Then, the proposed DRG meta-heuristic stands as the next most competitive method with a mean error gap of 2.58%. The VNS, SA and GA methods show a more significant mean error gap with 4.69%, 13.79% and 7.03% respectively. Also, VNS, SA and GA show a higher variability in performance which is characterized by an increased standard deviation when compared with their mean error gap. The DRG shows a balanced mean and standard deviation, thus its performance is more robust and consistent through all different CCCP instances.

DRG vs. CKM-TCG

Tables 3 and 4 present the results of the DRG meta-heuristic for the instances where reference data of the CKM and TCG methods were available. As presented in Table 3, the DRG meta-heuristic is more competitive than the CKM method. Also, as previously observed, the DRG is more consistent.

When compared to the TCG method, this is more competitive than the DRG approach even though the error gaps are minimal (less than 1.5%).

Error and speed vs. size of the instance

Figure 5 presents the graphical review of the error gaps of all methods based on the size of the test instance. TS and CS are located on the x-axis as they are the benchmark methods. It can be observed that, as the size of the instance increases, the error gap of SA, GA and VNS significantly increases. TCG presents very small error gaps with small instances (less than 1,000 points) and CKM reports error gaps comparable to SA for small-medium size instances (less than 5,000 points). In contrast, DRG performs consistently through all instances, decreasing its error gap as the instance grows from medium to large size (up to 13,0221 points).

Regarding speed, Fig. 6 presents the computational (running) times reported by the reviewed methods. While TS and CS are the benchmark methods, these take over 25,000 s to reach the best-known solution for the largest instance. Note that for these methods, their computational times exponentially increase for instances larger than 6,000 points.

In contrast, SA is very consistent with a computational time of approximately 1,000 s through all instances. GA significantly increases for instances larger than 6,000 points (up to 7,000 seconds for the largest instance). However, these methods have the largest error gaps as reviewed in Fig. 5. The speed pattern of DRG is very similar to that of GA, however, as reviewed in Fig. 5, its error gap is the closest to the benchmark methods for instances larger than 6,000 points.

It is important to mention that this comparison may not be fair due to the differences in the programming language and the hardware used for implementation and testing of all the methods. In order to compare running speed all methods should be tested with the same hardware and be implemented with the same programming language by the same software developer (the developer’s expertise may also affect the speed performance of the software). Due to the difficulty of achieving this task, in practice the comparison is commonly performed on the best results obtained by other methods as performed in Chaves & Nogueira-Lorena (2011) and Fernandes-Muritiba et al. (2012). Running time is measured in order to determine if the proposed method can provide a solution within reasonable time considering standard resources of hardware and software. Due to this situation, the DRG can provide very suitable results (<2.6% error gap) within very reasonable time.