A distributed algorithm for solving large-scale p-median problems using expectation maximization

p-median problem formulation

Linear integer programming (LIP) can be used to obtain optimal solutions for small scale p-median problems but these solutions quickly become time and memory prohibitive as the values of $n$ and $p$ increase. The LIP formulation (Daskin, 2013; Rosing, Revelle & Rosing-Vogelaar, 1979; Senne, Lorena & Pereira, 2005; Resende & Werneck, 2004) for the discrete p-median problem is described below.

LIP formulation

Decision variables:

1. $X=\{X_j\}$, $\mathrm{\forall }j\in \{1,2,3\dots n\}$

2. $Y=\{Y_{ij}\}$, $\mathrm{\forall }i\in \{1,2,3\dots n\}$ and $\mathrm{\forall }j\in \{1,2,3\dots n\}$

Minimize $Z={\mathrm{\Sigma }}_{i=1}^n{\mathrm{\Sigma }}_{j=1}^n{w}_i{Y}_{ij}{d}_{ij}$ such that:

Constraints:

1. ${\mathrm{\Sigma }}_{j=1}^n{Y}_{ij}=1$, $\mathrm{\forall }i\in \{1,2,3\dots n\}$

2. ${\mathrm{\Sigma }}_{j=1}^n{X}_j=p$, $\mathrm{\forall }j\in \{1,2,3\dots n\}$

3. $Y_{ij}-X_j\le 0$, $\mathrm{\forall }i\in \{1,2,3\dots n\}$, $\mathrm{\forall }j\in \{1,2,3\dots n\}$

4. $X_j=\{0,1\}$, $\mathrm{\forall }j\in \{1,2,3\dots n\}$,

5. $Y_{ij}=\{0,1\}$, $\mathrm{\forall }i\{1,2,3\dots n\}$ and $\mathrm{\forall }j\in \{1,2,3\dots n\}$

$X_j$ is a binary variable that, if set, represents the selection of destination $j$ as one of the facilities. $Y_{ij}$, also a binary variable, if set, represents the assignment of destination $i$ to facility $j$. Z, the objective function is the demand ( $w_i$) weighted sum of distances between facility $j$ and destination $i$. This distance is denoted by $d_{ij}$. The distance value $d_{ij}$ contributes to the sum only if $Y_{ij}$ is set. Constraint 1 ensures that each destination is assigned only to a single facility. The selection of exactly $p$ facilities is ensured by Constraint 2. The assignment variable $Y_{ij}$ cannot be set if a facility is not located at destination $j$, this condition is enforced by Constraint 3. This integer programming problem can be solved by using branch and cut algorithms on the continuous linear programming solution. The size of the search space for a p-median problem solution is equal to the Stirling number of the second kind (Resource, 2020) or the number of ways to partition a set of $n$ items into $p$ parts. The branch and bound tree in the LIP solution may degenerate into an exhaustive search leading to exponential time complexity in the worst case scenario (Papadimitriou & Steiglitz, 1982).

Literature review

Since the linear integer programming solution is not scalable, many heuristic algorithms have been proposed to yield approximate solutions. Interchange-based algorithms have been shown to yield the best results out of all heuristic algorithms. Teitz & Bart (1968) proposed the original exchange or interchange algorithm in 1968. Modifications to the original algorithm (FI) to reduce the number of exchanges were proposed in Goodchild & Noronha (1983), Whitaker (1983), and Hansen & Mladenović (1997). Facilities currently in the solution are interchanged with facilities not in the solution one at a time, and the solution after this interchange is evaluated. The interchange that results in the maximum decrease in the cost function is selected in each iteration. The interchange algorithm has been shown to produce close to optimal results, but even the FI algorithm is not feasible for larger datasets because of the large number of pairwise exchanges. The total number of swaps for this algorithm is equal to $(p)(n-p)$ per iteration, where $n$ is the total number of destinations and $p$ is the required number of facilities. For each swap, the facility being introduced is tested as a candidate for each destination, therefore each iteration is $O(p{n}^2)$. Overall complexity of the algorithm is $O(tp{n}^2)$ where $t$ is the number of iterations required for convergence/termination. The value of $t$ depends on the initial solution. Other common heuristic algorithms include the Global/Regional Interchange Algorithm (GRIA) (Densham & Rushton, 1992), the greedy addition algorithm (Kuehn & Hamburger, 1963) and the alternate selection and allocation algorithm (Maranzana, 1964). These algorithms are discussed and evaluated in detail in Gwalani, Tiwari & Mikler (2021).

Several meta-heuristic algorithms have also been proposed in the literature to improve the cost function values over those obtained from the interchange based algorithms. The most common meta-heuristics like genetic algorithms (Osman, Erhan & Zvi, 2003), tabu search (Rolland, Schilling & Current, 1997; Loranca, Velázquez & Analco, 2015) and simulated annealing (Chiyoshi & Galvão, 2000; Murray & Church, 1996) have been formulated and implemented for the p-median problem. Most of these meta-heuristic algorithms focus on improving the solutions obtained from local search or exchange algorithms in terms of the cost function value. The gains over the cost function value when compared to the classic exchange heuristics, however, are not substantial. The time required to compute the solutions on the other hand, is much higher. These methods do not scale well with respect to time as the numbers of destinations and facilities increase.

Taillard presents two decomposition methods to reduce run time, local optimization (LOPT) and decomposition (DEC) in Taillard (2003). The LOPT algorithm selects $p$ destinations as facilities randomly as the current solution. These facilities are added to a candidate list for subproblem creation. A subproblem is created in each iteration by selecting a facility $i$ at random, and its $r-1$ closest facilities from the candidate list. These $r$ facilities and demand points assigned to them form a subproblem that is solved using an interchange algorithm. The facilities selected from solving the subproblem are added to the candidate list if the solution (for the subproblem) has improved, else the original facility $i$ is removed from the candidate list. The algorithm terminates when the candidate list is empty. The DEC algorithm breaks the original problem into $t$ subproblems by selecting $t$ facilities in the original problem using an interchange algorithm. The service areas created for each of these $t$ facilities serve as the regions for the subproblems. The number of facilities, to be selected in each subproblem is equal to $int(p/t)$. The value of $t$ or the number of subproblems is set to $\lfloor \sqrt{p}\rfloor$ in this algorithm. It is shown in the study that the LOPT algorithm works better than DEC, however, it is sensitive to the initial solution and requires a considerable amount of time for termination. The results for large-scale problems ( $n=\text{85,900}$) in the study are shown relative to the best solution obtained from multiple runs of DEC, LOPT, and the ALT procedure. The ALT procedure is a variant of Maranzana’s alternate selection and allocation algorithm. The performance of these algorithms for large-scale problems was not compared with algorithms that have been known to produce close to optimal results. The decomposition into subproblems in the LOPT algorithm does not take into account the spatial distribution of the demand points. Therefore, the solutions to the subproblem may not solve the original problem optimally particularly for smaller values of $r$, and it becomes time prohibitive for larger values of $r$. The LOPT algorithm was shown to produce near-optimal solutions for small scale problems in their study. These decomposition-based techniques have the following limitations: (1) the number of subproblems is calculated using a formula that does not take into account the spatial configuration of the region, (2) the number of facilities to be selected in each subproblem is not demand weighted, and (3) the subproblems are generated by solving a p-median problem which may not successfully identify clusters in the region. The decomposition methods described in this study do not utilize the inherent clustering in the region, hence they can be improved upon both with respect to the cost function value and execution time.

Density-based decomposition

Clusters in a region can be identified by detecting changes in the spatial distribution of demand (population) in the region. The central place theory dictates that consumer preferences in any geographic area will lead to centers of various sizes to emerge in a region (Daniels, 2017). Therefore, any region is likely to have population centers with different population densities. Behavioral patterns also often lead to formation of clusters in a spatial region (Rushton, 1969). The goal of the density based decomposition is to identify these centers and assign the demand points to its corresponding center. We model any region as being composed of different sub-regions with corresponding mean and standard deviation of population points distributed normally around the mean and standard deviation. This Gaussian mixture model (GMM) captures the inter and intra subregion density variation adequately as it is simpler and faster than other density based clustering algorithms while performing well in identifying the dense regions. The expectation-maximization (EM) clustering algorithm (Dempster, Laird & Rubin, 1977) has been known to yield good results for detecting Gaussian clusters with different parameters in a dataset. It is an iterative method that attempts to predict parameters for statistical models in two steps: (1) expectation, which calculates the log-likelihood corresponding to the current estimates for the parameters, and (2) maximization, in which the model changes the parameters to increase the log-likelihood value. Since the number of clusters is an input to the algorithm, EM can be tried to fit models with a varying number of clusters and each fit can then be evaluated using the Bayesian information criterion (BIC). BIC can be used to identify the minimum descriptor length that can describe the data (Wikipedia.org, 2019). It penalizes the addition of new parameters to the model to avoid overfitting. We use expectation-maximization to fit the data in a GMM with $c$ components, starting from $c=2$ to $c=c_{max}$. Each demand point, $d_i$, represents a vector of size two, $\overrightarrow{x_i}=[x_{di},y_{di}]$. The objective of the Expectation-Maximization algorithm is to obtain a GMM for the data that represents the normally distributed sub-populations in the overall population (Brilliant.org, 2019). The population refers to the set of all demand points in the region, while each sub-population is a set of demand points that belong to the same bivariate normal distribution according to the model. Each component of the GMM has a mean vector, $\overrightarrow{{\mu }_i}$ and a covariance matrix ${\mathrm{\Sigma }}_i$. The probability distribution for the model is given by Brilliant.org (2019):

(1) $$p(\overrightarrow{x})={\mathrm{\Sigma }}_{i=1}^c{\theta }_i\mathcal{N}(\overrightarrow{x}|\overrightarrow{{\mu }_i},{\mathrm{\Sigma }}_i)$$where $\mathcal{N}(\overrightarrow{x}|\overrightarrow{{\mu }_i},{\mathrm{\Sigma }}_i)$ represents a bi-variate Gaussian distribution (Eq. (2)) corresponding to the component $i$ with parameters, $(\overrightarrow{{\mu }_i},{\mathrm{\Sigma }}_i)$. The sum of the component weights corresponding to each component, ${\theta }_i$ is equal to one, ${\mathrm{\Sigma }}_{i=1}^c{\theta }_i=1$.

(2) $$\mathcal{N}(\overrightarrow{x}|\overrightarrow{{\mu }_i},{\mathrm{\Sigma }}_i)=\frac{1}{\sqrt{{(2\pi )}^c|{\mathrm{\Sigma }}_i|}}{e}^{-\frac{1}{2}{(\overrightarrow{x}-\overrightarrow{{\mu }_i})}^T{\mathrm{\Sigma }}_i^{-1}(\overrightarrow{x}-\overrightarrow{{\mu }_i})}$$

Given a starting value for parameters, $({\theta }_i,{\mu }_i,{\mathrm{\Sigma }}_i)$ for each component in the model, the expectation step calculates the probabilities/expectation that a data point belongs to that component. The maximization step updates the parameters to maximize the overall expectation or likelihood value. Equation (3) shows the calculation of BIC for a model (Wikipedia.org, 2019; Neath & Cavanaugh, 2012).

(3) $$BIC=ln(n)c_m-2ln((\hat{L}))$$

$(\hat{L})$ is the maximized value of the likelihood function after the EM algorithm converges, $n$ is the number of data points, $c_m$ is the total number of model parameters. To detect the high demand point density sub-regions in the area, each demand point coordinate is input to the EM algorithm. The EM algorithm is executed for $c=2$ to $c=c_{max}$, and the model corresponding to each $M(c)$ is evaluated using BIC. It was observed that BIC analysis becomes unreliable for higher values of $c$. The gain in the maximum likelihood optimization function is higher than the penalty corresponding to the increase in the number of parameters. This results in decomposition of regions into independent components that maybe served by a common facility which is detrimental to the p-median objective function. To avoid this over decomposition, the $c_{max}$ value was set at 20 in these experiments.

Models with lower BIC values are better than models with higher BIC values, hence the value of $c$ for which the minimum value of BIC is reported is selected as the number of clusters, $k$, for the region. The region is decomposed into $k$ subregions or components using the EM algorithm. Each subregion has a different $\overrightarrow{\mu }$ and $\overrightarrow{\mathrm{\Sigma }}$ corresponding to the positions of demand points. Each of these subregions, however, may not represent a dense or clustered region as the EM detects changes in distribution. These distributions could be homogeneous or non-homogeneous. A homogeneous distribution represents a random scatter of demand points in the region. A heterogeneous distribution is a significant deviation from a homogeneous distribution of points in the region. Figure 4 shows an example of the creation of subproblems using the EM algorithm. The initial demand points, membership of demand points to their corresponding Gaussian component, and the corresponding subregions are shown in Figs. 4A–4C respectively. The destinations are distributed in 20 components based on the BIC analysis. The subregions are created by unionizing the geometries corresponding to the demand points.

Density-based decomposition: merge step

It can be seen in Fig. 4 that 10 out the 20 subregions do not represent clustered or dense areas in terms of the demand point distribution, and therefore may not be good candidates for an independent subproblem. We propose a merge step to merge the sparser and/or homogeneous regions into denser ones, thereby preventing them from being solved independently.

A subregion is labeled sparse if the demand density in the subregion is less than the demand density of the entire region. Demand density is defined as the total demand per unit area. The total demand in a subregion is calculated by aggregating the demand at all of its demand points. The heterogeneity (or homogeneity) of the spatial distribution of points in a region can be detected by comparing the existing distribution of the subregion with the expected binomial distribution using a ${\chi }^2$ test. The details of heterogeneity detection are explained in the appendix “Detecting Heterogeneity in a Region”. Once the non-clustered (homogeneous or sparse) subregions have been identified, the demand points corresponding to these subregions are reassigned to the remaining dense subregions. Each demand point in a non-clustered subregion is assigned to its closest dense subregion. The distance to a dense subregion is measured as the distance between the demand point to be reassigned and the centroid of the subregion. Thus, the non-clustered regions are dissolved in this way reducing the number of subregions and subproblems from $k$ to $q$. Figure 4D shows the subregions corresponding to Fig. 4C after the demand points corresponding to the non-clustered subregions are merged into the dense subregions. We propose to use these subregions as independent subproblems.

If the distribution of demand points in the original problem is completely homogeneous, the merge step will lead to a single subproblem corresponding to the original problem, hence the decomposition will not succeed. This limitation of the algorithm is discussed more in “Challenges and Distance-based Decomposition”.

Divide and conquer

The distribution or density-based decomposition of the region yields $q$ subregions after the merge step is employed. The number of facilities to be selected in these subregions is determined based on the proportion of the population in the subregion with respect to the entire population. The facilities are assigned to each subregion proportional to the total demand in the subregion (integer division). The remaining facilities are assigned to the subregions in a round-robin fashion in decreasing order of demand at subregions. These subproblems are then solved exactly using MILP² solver if possible. If the scale of the subproblem is too large to solve using MILP, then the subproblem is solved approximately using FI. Scale based memory profiling was conducted to determine the threshold values of $n$ and $p$ to employ FI instead of MILP. The subproblem solutions are combined to obtain a global solution in the reduce step. The union of facilities selected in each subproblem forms the current solution, X, for the original problem. The $n$ demand points in the original problem are assigned to their closest facility in the current solution. This solution then can be further improved by applying an exchange based search algorithm (FI/GRIA) to this solution. Since the solution is near the optima, the search takes much less time than when used on a randomly chosen initial solution.

Time complexity

The time complexity for a serial version of FI is $O(p^2{n}^2)$, assuming $t$, the number of iterations, $t=cp$, for some constant $c,$ where $c<1$. This assumption is based on experimental analysis on the change in number of iterations with change in initial solution in Gwalani, Tiwari & Mikler (2021). If the problem is decomposed into $q$ parts, then on average each subproblem is of scale ( $n/q,p/q)$, hence the complexity reduces to $p^2{n}^2/q^4$ for all problems which are solved concurrently. The expectation maximization step is $O(n(2+3+\dots c_{max}))=O(n)$, for $c_{max}<<n$ (the number of iterations needed for convergence are much smaller than $n$, hence that factor is neglected). The reassignment step is another $O(n)$ operation. The improvement step at the end is again $O(tp{n}^2)$ but the number of iterations required for convergence is much smaller, hence the improvement term reduces to $O(p{n}^2)$. We show in the article, that a reasonably good solution can be obtained even without the improvement step by using distributed computing.

Synthetic datasets

We compare the performance of the discussed algorithm on synthetic datasets with varying scale and spatial distributions. In the synthetic datasets, the number of demand points is varied by changing the total demand in the region. All demand units are scattered in the region grid depending on the required spatial distribution. These demand units are then combined to create destinations. Each destination is assumed to have a demand selected from a normal distribution, $\mathcal{N}\text{(1,500, 400)}$. These values were chosen to replicate the creation of census block groups, a census geographic unit in USA. The centroid of the merged demand units for a destination defines the coordinates of the demand point. Because of the stochastic nature of the this process, not all points on the fine grid get assigned to a demand unit and in turn to a destination, this results in the blank unoccupied islands within the region. The process is explained in detail in Gwalani, Tiwari & Mikler (2021). The value of $p$ is chosen to be approximately equal to $n/4$ to create large-scale problems.

We used three synthetic datasets to evaluate the EM-FI algorithm for proof of concept. Figure 5 shows the demand distribution in these datasets, the cost function and execution time results obtained for these datasets using the three heuristic approaches. It can be seen that the cost function values are reasonably good (within 2% of the best-known solution) for EM-FI even when a global improvement is not employed in the reduce step for the datasets with multiple clusters. Additionally, the EM-FI strategy is at least five-eight times faster than FI. The EM-FI:Reassignment does not perform as well for Dataset III because it contains only one large cluster and the decomposition breaks it into multiple parts even when they belong to the same GMM component which results in the poorer results. Using EM-FI:Improvement however, yields results that are within 0.5% of the best-known solution with a much lower execution time than the stochastic version of FI. In fact, statistical analyses using the Man-Whitney U test³ showed that the cost function values obtained from EM-FI:Improvement are statistically significantly lower than the cost function values obtained using FI.

Real datasets

We use response planning in preparation for bio-emergencies as an application of the p-median location model to evaluate the algorithm on real datasets. Public health emergencies can be catastrophic to the health and well-being of local communities. These emergencies require a timely and well-planned response from the local public health authorities to minimize casualties. The Centers for Disease Control and Prevention recommends setting up ad-hoc clinics or dispensing locations to distribute supplies and medical countermeasures (MCMs) to the affected population within a stipulated time frame (Centers for Disease Control and Prevention, 2006). We formulate this problem as the p-median problem where census block group centroids in the region represent population centers or demand points in the region, and the dispensing locations are the facilities that need to be located. The results are demonstrated for three regions in the United States with varying spatial distribution of the census block group centroids, 1) Dallas County, Texas, 2) Texas Department of State Health Services Region 6, Texas and 3) Los Angeles County, California. Figure 6 shows the demand point distributions, problem scales and results for these datasets. It can be seen that since there are more obvious distribution changes and the clusters are more spread out and distinct for Dallas and Los Angeles counties, EM-FI:Reassignment is able to produce reasonably good results (within 5% of the best-known solutions) with respect to the cost function values for these datasets. EM-FI:Improvement reduces the cost function values further for these datasets with relatively low execution time gain. The block group centroids in Region 6 are mostly centered around the city of Houston (Fig. 6B). The minimum BIC value was obtained for 18 components for this region. Figure 7A shows the decomposition of block groups in Region 6 into 18 parts. The components after the merge step is performed on the sparser and homogeneous regions are shown in Fig. 7B. EM-FI:Reassignment produces poorer results because of the decomposition of the single large component in the center into subregions. Consequently, this causes the EM-FI:Improvement algorithm to take more time as the initial solution in the reduce step is not that close to the optimal solutions (local or global). Despite this limitation, the EM-FI:Improvement algorithm outperforms FI with respect to both execution time and cost function values. Table 2 shows the comparison of results with respect to both execution time and cost function values between FI and EM-FI:Reassignment for these datasets. It can be seen that the execution time is almost two orders of magnitude less for the EM-FI:Reassignemnt when compared to FI even on a personal quad core machine. The execution time column in the table shows the average time across 30 runs and the first and third quartile values. The EM-FI:reassignment produces reasonably good solutions for clustered datasets even with respect to the cost function values (within 4% of the best solution). The cost column in the table shows the minimum relative cost and the first and third quartiles values for relative costs.

Challenges and distance-based decomposition

There are three drawbacks for using the EM-FI:Improvement algorithm for any general large-scale problem. (1) Over decomposition of a region with a single large cluster into subregions as a result of the EM-BIC analysis. The merge step may not combine the subregions belonging to the same dense cluster and lead to over-decomposition, as illustrated in Fig. 7. This over-decomposition yields to poorer results for this kind of distribution for EM-FI:Reassignment and increases the run time for EM-FI:Improvement. (2) Additionally, even for a multi-clustered, large problem, the final step in EM-FI improvement can prove to be a bottleneck in terms of execution time because of the size of the problem, even if the solution after EM-FI reassignment is close to the optimal solution. (3) For an entirely homogeneous distribution, the density based decomposition yields the original problem after the merge step, hence there are no time gains.

These three problems can be resolved by replacing the EM-FI improvement step with another decomposition step. The objective of this level of decomposition is to divide the current sources selected after the divide and conquer step (EM-FI reassignment) into $k$ parts using a distance based clustering algorithm like k-means. The idea behind this decomposition is that the destinations that are far away from a source are unlikely to be served by it given the distance based optimization function. The value of $k$ is determined by using the cluster quality criterion, Davies–Bouldin Index (Davies & Bouldin, 1979), on clusters created for $k=2$ to $k=k_{max}$. The Davies–Bouldin Index selects the clustering solution that minimizes the ratio of intra-cluster dissimilarity (scatter) and inter-cluster dissimilarity. The index evaluates the ratio, $R_{ij}$, for all pairs of clusters, $C_i$ and $C_j$, $i\ne j$, in the solution. The intra-cluster scatter or the numerator of $R_{ij}$ is the sum of the average distance between all members in cluster $C_i$ to the centroid of $C_i$ and the average distance between all members in cluster $C_j$ to the centroid of $C_j$. While the inter-cluster dissimilarity between $C_i$ and $C_j$ is the distance between the two cluster centroids. For each cluster $C_i$, the maximum value of $R_{ij}$ is selected as its index, the Davies–Bouldin measure for the solution is the average index value across all clusters. A lower value of this index indicates better clustering. Each part is treated as an independent subproblem with the sources in that part as the initial solution to the subproblem for serving the union of destinations assigned to them. These subproblems are solved using FI concurrently, and since the initial solutions are already close to the optima these subproblems converge sooner. This step of distance based decomposition is repeated for $k=k-1$ (until $k=1$) by clustering the sources obtained from solving the subproblems created in the previous iteration. If the objective cost function value remains unchanged across two consecutive distance based decomposition iterations, then the process is terminated before $k=1$. The final step in this process for $k=1$ (if the algorithm does not terminate before that) is identical to the EM-FI Improvement step, however, the initial solution is now even closer to the optimum because of the previous decompositions, and hence it is much faster. We conducted experiments for the real dataset problems (10 runs for each dataset) on a Google compute machine (n1-highmem-16, 16 vCPUs, 104 GB RAM) employing the iterative distance based decomposition after EM-FI:Reassignment instead of the improvement step. We executed six MILP problems and 10 FI subproblems concurrently on this machine because of higher memory and CPU capabilities. Figure 8 shows the changes in execution time and relative cost function values with each iteration of distance based decomposition, after the reassignment step in a single execution for Region 6. Table 3 shows the cost and execution time results for these experiments and compares them with the results obtained from EM-FI:Improvement on the i7 dual quad core machine. The time gains from both the use of a more powerful machine and distance-based decomposition are shown in the table. Additionally, we solved a p-median problem to select 2,500 facilities in the entire state of Texas on the Google machine. There are 15,811 census block groups in the state of Texas. A problem of this scale would have taken days to execute if solved using the serial implementation of FI. A solution within 4% of the final solution was obtained in about 1.5 h for this problem and the final solution was obtained in less than 9 h.

If the density based decomposition yields the original problem because of the absence of clusters, the region can be divided into $k$ subproblems using the distance-based decomposition strategy described above. The first iteration solution can be improved by dividing the sources obtained from the first iteration into $k=k-1$ parts iteratively. These decomposition strategies, however, may yield a single large subproblem that becomes the bottleneck to obtain a solution. This was observed in the case of the Texas dataset which caused the execution time to increase. The algorithm can be used recursively to overcome this problem.

Another challenge with the EM-FI decomposition algorithm is the unreliable nature of the BIC criterion. Higher numbers of clusters seemed to yield lower BIC values and lead to over-decomposition, especially in case of centered distributions. There are many cluster quality metrics and criteria that measure the intra-cluster similarity and the inter-cluster similarity in different ways. These metrics include cluster validation indexes such as Calinski and Harabasz score (Caliński & Harabasz, 1974), Dunn’s Index (Dunn, 1973), Davies–Bouldin index (Davies & Bouldin, 1979), and silhouette coefficient (Rousseeuw, 1987), if the ground truth is not known. It was observed that these indexes yield vastly different results with varying distributions. BIC seemed to yield the best results for identifying density-based clusters and the Davies–Bouldin Index was selected to evaluate distance-based clusters. For a given region, denser regions can be identified by using DBSCAN without the dependence on the number of clusters, but this algorithm requires two additional parameters, minimum number of neighbors and the radius of neighborhood, which are not generalizable across problems as density is a relative concept. Additionally, the EM algorithm is faster than DBSCAN, hence more suitable for efficiency gains. The algorithm can be made even more efficient by employing parallel versions of expectation-maximization (Lee, Leemaqz & McLachlan, 2016) if necessary, although the expectation maximization step in the algorithm is much faster than the p-median steps, therefore the time gains are expected to be minimal.

Detecting heterogeneity in a region

The heterogeneity (or homogeneity) of the spatial distribution of points in a region is identified by overlaying a fine grid over the region. Let the number of demand points in the subregion be $n_{sub}$. The expected distribution for the number of cells in the grid with $i$ ( $0\le i\le n_{sub}$) demand points is compared with the obtained frequency distribution using the ${\chi }^2$-statistic. The probability that a demand point is located in any cell in a homogeneous distribution is $p_{cell}=1/N$, N is the total number of equal-sized cells in the grid. The probability that a cell will be occupied by $i$ demand points, $p(i)$ is equal to $p(i)=\left(\begin{array}{c}{n}_{sub}\\ i\end{array}\right)p_{cell}^i(1-p_{cell}{)}^{n_{sub}-i}$. Therefore, the expected number of cells with $i$ demand points is $E(i)=Np(i)$. This distribution is compared with the actual number of cells with $i$ demand points, $F(i)$, in the subregion to check for deviation from a random distribution. The ${\chi }^2$ statistic can then be calculated as ${\chi }^2={\mathrm{\Sigma }}_{i=0}^{n_{sub}}\frac{{(E(i)-F(i))}^2}{E(i)}$. This ${\chi }^2$ probability distribution has $n_{sub}$ degrees of freedom. The homogeneity hypothesis is evaluated by calculating the probability that the ${\chi }^2$ value lies on the chi-square distribution function. We use the significance level of $5\mathrm{\%}$ for this analysis ( $\alpha =0.05$).