1 Introduction↩︎

Clustering is of fundamental importance for data analysis, with ubiquitous applications in various areas including pattern recognition [1], bioinformatics and biomedicine [2], and data management [3], [4]. Broadly, clustering a dataset requires grouping elements with high similarity, and separating dissimilar elements [5]. The clustering of a dataset is often viewed as an optimization problem, where the objective varies according to the desired properties of the optimal clustering. A large number of algorithms have been designed to cluster a dataset according to various objectives [6]. In particular, recent research for clustering methods has its focus on the design of rigorous and scalable approaches to cluster massive datasets [7]–[12].

The vast number of clustering objectives and methods, introduce the need for rigorous methods to identify “good quality” clusterings of a dataset [13], [14]. Clustering evaluation (or simply validation) assesses the quality of a given clustering relying on an evaluation measure. An evaluation measure is either external (i.e., supervised) or internal (i.e., unsupervised). Supervised measures rely on external knowledge about the data, e.g., ground truth labels—rarely available [15]. In contrast, internal measures evaluate the quality of a clustering relying uniquely on the dataset. Due to their broad applicability, internal measures are widely used in practice to identify high-quality clusterings [16].

The silhouette is one of the most commonly used internal measures for clustering evaluation [16], [17]. The silhouette \(s(e)\) of an element \(e\) assigned to a cluster \(C\) is defined as the ratio \((b(e)-a(e))/\max\{a(e),b(e)\}\) where \(a(e)\) is the average distance of \(e\) from the other elements of \(C\), and \(b(e)\) is the minimum average distance of \(e\) from the elements of a cluster \(C'\) different from \(C\). In other words, \(s(e)\) provides a simple to interpret score between -1 and 1, which captures the quality of the assignment of \(e\) to its group \(C\). That is, if \(e\) is closer (on average) to the elements in its cluster \(C\) than to the elements in the “closest” cluster \(C'\neq C\), then \(s(e)\) will be close to 1. In contrast, \(s(e)\) becomes negative when the element \(e\) is closer to the elements in \(C'\). As an example, consider elements \(e_{13}\) and \(e_{14}\) in [subfig:dataset] (where the distance is the Euclidean distance): \(i\)) the silhouette \(s(e_{13})\) is positive since \(e_{13}\) is closer to the elements of its cluster (\(e_{11}\) and \(e_{12}\)) compared to its average distance from the elements of other clusters; \(ii\)) the silhouette \(s(e_{14})\) is close to -1 since \(e_{14}\) is very close to all elements of cluster \(C_3\) and distant from the other elements of its cluster \(C_1\). Hence, the silhouette \(s(e)\) of an element \(e\) (or local silhouette) provides an interpretable and reliable score to assess the quality of the assignment of each element \(e\) to its cluster [18], [19].

Furthermore, the average silhouette of all elements, known as the silhouette of a clustering (or global silhouette), has a plethora of applications in many important areas [20]–[23], and is described in widely used data mining textbooks [15], [24]. As an example, the silhouette of a clustering can be used to identify a good value for the number \(k\) of clusters computed by popular clustering algorithms, e.g., \(k\)-means, or \(k\)-medoids [19], [25].

For general distance metrics, the exact computation of the silhouette for a \(k\)-clustering of \(n\) elements requires \(\Theta(n^2)\) distance calculations, irrespective of \(k\). The quadratic barrier constitutes a key challenge that prevents silhouette-based analyses on modern massive datasets. Surprisingly, for such a general setting, there are no methods to efficiently approximate the silhouette with provably high accuracy guarantees and breaking the \(\Theta(n^2)\) barrier. In contrast, various scalable techniques have been proposed to efficiently cluster large datasets with rigorous guarantees, especially in distributed environments [8], [26]–[29].

In this work, we introduce the first randomized approximation algorithms with rigorous (probabilistic) guarantees on the quality of their estimates for two key problems: \(i\)) estimating the silhouette of each element \(e\) of a given clustering; \(ii\)) estimating the global silhouette of a clustering. Our algorithms feature a memory-efficient distributed implementation in popular frameworks such as MapReduce [30]–[32], and the Massively Parallel Computing (MPC) model [33], in analogy with state-of-the-art scalable clustering algorithms. Similar to other randomized algorithms, the guarantees of our methods are controlled by two easy to interpret parameters \(\varepsilon\) and \(\delta\), controlling the trade-off between efficiency and accuracy. In particular, \(\varepsilon\) controls the (absolute) error between the output of our methods and the actual values being estimated, and \(\delta\) controls the error probability (all our results hold with probability at least \(1-\delta\)).

For the estimation of the silhouette of each element \(e\), we rely on the Probability Proportional to Sampling (pps) technique by [34], which can be used to obtain highly accurate estimates of the average distance between \(e\) and the elements of a given cluster. For a \(k\)-clustering, our method performs \(\mathcal{O}(nk\varepsilon^{-2} \log(nk/\delta))\) distance computations, producing estimates \(\widehat{s}\,(e)\) such that \(|s(e)-\widehat{s}\,(e)| \le f(\varepsilon)\) with \(f(\varepsilon) = \mathcal{O}\!\left(\varepsilon\right)\) with probability at least \(1-\delta\), simultaneously for all elements \(e\). Note the substantial improvement over the \(\Theta(n^2)\) distance computations of an exact method, at the expense of a small estimation error.

For the global silhouette estimation, we design three estimators, including a simple algorithm that computes the average silhouette using a small random subset of elements, and two slightly more complex algorithms relying on the pps approach. We show in our analysis that all our methods perform \(o(n^2)\) distance computations, and report an estimate with bounded error from the actual global silhouette, with controlled error probability. Finally, our experiments show that all our methods are extremely efficient and scalable in practice, especially on large datasets, where exact approaches are infeasible and existing techniques perform poorly.

Our contributions. To summarize our contributions are as follows.

• We present the first rigorous sampling-based approximation algorithms for two key problems involving the silhouette: \(i\)) estimating the silhouette of all elements of a clustered dataset; \(ii\)) estimating the global silhouette of a clustering. Our methods require the distance between elements to be defined through a metric.

• We prove that all our methods require a significantly smaller number of distance computations than \(\Theta(n^2)\), reporting high-quality estimates with (absolute) error controlled by parameter \(\varepsilon\), with probability controlled by parameter \(\delta\).

• We provide a distributed MapReduce design of our methods and discuss how our methods can be ported to the popular MPC framework. Our distributed methods use sublinear local memory (at each worker) and a constant number of rounds, two features that are highly desirable when processing massive datasets.

• We perform an extensive experimental evaluation on medium- and large-size datasets with multiple goals: \(i\)) assess the trade-off between the quality of the estimates and the number of distance computations of our methods for estimating the global silhouette of the entire clustering; \(ii\)) evaluate the accuracy of our method for the estimation of the local silhouette for all elements; \(iii\)) assess the performance of our distributed methods; \(iv\)) illustrate applications of our new methods, leveraging both the global and local silhouette.

Contributions over the conference version. A preliminary version of this work was presented at the SIAM International Conference on Data Mining, 2021, authored by a subset of the authors [35]. This work feartures substantial modifications and additions to that preliminary work, as summarized below.

• We propose two novel estimators for the global silhouette coefficient, in addition to the estimator originally devised in [35]. Our new methods require considerably fewer samples compared with the method appeared in [35], yielding a reduction in the time complexity for global estimation.

• We rigorously study the problem of local silhouette estimation, providing:

• We include all proofs omitted from [35] for space constraints.

• We show how the distributed MapReduce implementation of our methods can be ported to the popular MPC framework, which is widely adopted for large-scale data analysis.

• We perform a new and extensive experimental evaluation, over medium- and large-size real datasets available online, showing the effectiveness of our approaches in various settings.

The paper is structured as follows. In 2, we introduce the necessary preliminary notions. In 3, we discuss relevant previous work, with a focus on the silhouette approximation. In 4, we introduce our methods and their analysis. In 5, we present the results of our extensive experimental evaluation. Lastly, 6 summarizes our concluding remarks.

The code and necessary material to reproduce our results is available at the following repository https://github.com/iliesarpe/ScalableSilhouetteComputation.

2 Preliminaries↩︎

Let \(U\) be a metric space equipped with a distance function \(d(\cdot,\cdot)\), and let \(V = \{e_1, \dots e_n\} \subseteq U\) be a dataset of \(n\) distinct elements from \(U\). Let \(\mathcal{C}= \{C_1, \dots C_k\}\) be a \(k\)-clustering¹ of \(V\), that is \(\mathcal{C}\) consist of a partition of \(V\) into \(k\) disjoint and non-empty subsets called clusters. Formally, the silhouette \(s(e)\) of an element \(e \in V\) belonging to a cluster \(C\in \mathcal{C}\), as introduced by [17] is defined as, \[\label{eq:localSilh} s(e) \doteq \frac{b(e)-a(e)}{\max\{a(e), b(e)\}}\tag{1}\] where \[\label{eq:AandBterms} a(e) = \frac{\sum_{e' \in C} d(e,e')}{|C|-1} \enspace, \qquad b(e) = \min_{\substack{C_j \in \mathcal{C}, \\ C_j \neq C}}\frac{\sum_{e' \in C_j}d(e,e')}{|C_j|}\enspace.\tag{2}\] That is, \(a(e)\) is the average distance of element \(e\in V\) from all other elements of its cluster \(C\); while \(b(e)\) is the minimum average distance of \(e\in C\) from the elements in some cluster \(C_j\in \mathcal{C}\) different from \(C\). Note that by 1 it holds that \(s(e) \in [-1,1]\subseteq \mathbb{R}\). If element \(e\) is the only element in its cluster \(C\) such that \(C=\{e\}\), we define \(s(e)=0\), as proposed by [17]. As an example, consider \(e_{15}\) in 1, then \(s(e_{15})=0\).

The silhouette \(s(e)\) of an element \(e\in E\) is a local measure that evaluates the quality of a clustering \(\mathcal{C}\) with respect to the assignment of \(e\) to its cluster \(C\). Under a good assignment, the average distance of \(e\) to the points in \(C\) is much smaller than the average distance to the points of a different cluster, thus from 1 it holds \(s(e) \approx 1\). Hence, high-values of \(s(e)\) denote a “good” assignment of element \(e\) to its cluster \(C\). Conversely, if \(e\) is, on average, much closer to the elements of a different cluster \(C_j\neq C\) then by 1 it holds \(s(e)\approx -1\), i.e., a “bad” cluster assignment for \(e\). The individual \(s(e)\)’s are often represented through the silhouette plot [17], where the abscissae range in \([-1,1]\) (possible silhouette values) and the ordinate values correspond to the distinct elements, ordered by clusters.

The silhouette of a clustering \(\mathcal{C}\) evaluates the quality of a clustering \(\mathcal{C}\) as the average silhouette of all the elements \(e\in V\). Formally, the silhouette \(s(\mathcal{C})\) of a clustering \(\mathcal{C}\) is \[\label{eq:avgShil} s(\mathcal{C}) \doteq \frac{1}{n} \sum_{i = 1}^{n}s(e_i)\enspace.\tag{3}\]

Therefore, \(s(\mathcal{C})\) provides a global measure of the quality of the whole clustering \(\mathcal{C}\), where a value close to 1 denotes a high quality of the clustering \(\mathcal{C}\).

In practice, the score \(s(\mathcal{C})\) is widely used to identify a good value for \(k\) (e.g., [36]). That is, analysts often compare various clusterings \(\{\mathcal{C}_{k_i}\}_{i\ge 1}\) obtained with different values of \(k_i\), selecting the best value of \(k\) as the one maximizing \(s(\mathcal{C}_{k_i})\), for \(i\ge 1\). In fact, using the silhouette to select a good value of \(k\) for a \(k\)-clustering is highly effective and reliable, in contrast to other popular methods, such as the elbow method, which often lead to poor results [13].

It is clear from 1 that the straightforward computation of \(s(e)\) for each element \(e\in V\) and the subsequent computation of \(s(\mathcal{C})\) require \(\Theta(n^2)\) distance calculations. We remark that such a quadratic complexity is extremely resource demanding even when dealing with medium-sized datasets, i.e., impractical for modern massive data volumes.

Finally, we wish to recall that for an algorithm employing random sampling, its sample complexity refers to the total number of samples generated throughout its execution [37].

3 Related work↩︎

Since the introduction of the silhouette metric by [17], its use has been ubiquitous for clustering applications [38]–[41], including the use of the silhouette to obtain new clustering objectives for both \(k\)-medoids [25], [42] and \(k\)-means [43]. However, as observed before, obtaining the exact silhouette is impractical for large datasets, due to the \(\Theta(n^2)\) distance computations required. Therefore, various methods have been proposed to reduce the time complexity, or to simplify the silhouette calculation under certain assumptions. In this section, we review to the most relevant methods for silhouette approximation, which mostly apply to center-based clusterings. That is, clusterings where each cluster is associated with a distinguished center from the metric space (e.g., the centroids in the case of \(k\)-means clusterings).

[44] introduced a variant of the silhouette \(s(e)\), where the term \(b(e)\) is computed as the distance from the closest center, other than the cluster \(e\) belongs to (while \(a(e)\) is computed according to 2 ). Following such idea, [45] introduced the simplified silhouette. The simplified silhouette for an element \(e\) assigned to cluster \(C\) evaluates \(a(e)\) as the distance between \(e\) and the center of \(C\), and \(b(e)\) as the distance between \(e\) and the closest center of a cluster different from \(C\). The complexity of computing the simplified silhouette reduces to \(\mathcal{O}(nk)\). [45] and [46] show empirically that the simplified silhouette can be an effective evaluation measure for clusterings returned by Lloyd’s algorithm [47]. However, there is no evidence of the effectiveness of the simplified silhouette for other clusterings (e.g, based on arbitrary distance functions). In addition, as shown in our experimental evaluation in 5, the difference between the values of the silhouette (1 ) and the simplified silhouette, can be arbitrary large, making the simplified silhouette unreliable for rigorous silhouette-based analyses.

A different heuristic for improving the computation of the exact silhouette of center-based clusterings under Euclidean distances was introduced by [48]. For each element \(e\) of a cluster \(C\), the term \(a(e)\) is computed according to its definition (2 ), while for \(b(e)\) the heuristic first determines the average distance \(\bar{d}_{e,C'}\) between \(e\) and the elements of the cluster \(C' \neq C\), whose centroid is closest to \(e\). If, for any other cluster \(C'' \not\in \{C,C'\}\), it holds that the distance between \(e\) and the centroid of \(C''\) is greater than or equal to \(\bar{d}_{e,C'}\), then it sets \(b(e)=\bar{d}_{e,C'}\), otherwise it computes \(b(e)\) according to 2 . For Euclidean distances, it can be shown that in the former case \(\bar{d}_{e,C'}\) coincides exactly with the value of \(b(e)\) defined in 2 . The worst case complexity remains quadratic, but the heuristic likely performs less than \(\Theta\!\left(n^2\right)\) distance computations in practice.

The Apache Spark programming framework provides optimized methods for computing the silhouette of a clustering under \(d\)-dimensional squared Euclidean distances and under one formulation of cosine distance.³ For such specific distances, simple algebra suffices to show that the pre-computation of a limited number of values dependent on the coordinates of the points in each of the \(k\) clusters yields a parallelizable algorithm performing \(\mathcal{O}(nkd)\) distance computations. However, the above optimization does not apply to arbitrary distance metrics, as of interest in our work.

The methods presented in this paper apply to any metric distance and provide provably accurate silhouette estimations. Our methods rely on the use a Probability Proportional to Size (pps) sampling scheme, where each element is sampled with a probability proportional to a certain “size” measure. The use of pps sampling has been pioneered in the context of center-based clustering algorithms, distance query processing and centrality estimation in graphs [27], [29], [34], [49]. To the best of our knowledge, prior to the conference version of the current work [35] the use of pps for efficient clustering evaluation had not been explored.

4 Methods↩︎

In this section we present our methods. In 4.1 we focus on the estimation of the silhouette of all individual elements, while in 4.2 we focus on the estimation of the silhouette of the whole clustering. Finally, 4.3 describes scalable distributed implementations of our methods on the MapReduce and MPC frameworks.

4.1 Estimating the silhouette of each element↩︎

In this section we address the following problem,

Our main result, addressing 1 is the following.

1 shows that we can obtain a very accurate probabilistic approximation of each silhouette value \(s(e), e\in V\). For practical values of \(k\ll n\), and constant \(\varepsilon\) and \(\delta\), the result of 1 breaks the quadratic barrier of \(\Theta(n^2)\) distance computations, enabling the efficient and practical computation of silhouette values over large and massive datasets. Note that the \(\mathcal{O}\!\left(\varepsilon^{-2}\right)\) factor in the sample complexity and in the number of distance computations from 1, might become large for small \(\varepsilon\). While a quadratic dependence on \(\varepsilon\) cannot be avoided in a worst-case analysis [37], in practice, even large values of \(\varepsilon\) (e.g., close to 1) yield highly accurate estimates for all values in \(\widehat{s}\,(V)\).

In this section, we also provide evidence that pps is crucial to achieve our result. Namely, we argue that, under the same sample complexity of 1, a more intuitive uniform sampling approach performs arbitrarily poorly.

2 indicates that uniform sampling cannot guarantee an (arbitrarily small) estimation error \(\mathcal{O}\!\left(\varepsilon\right)\) using a sublinear sample size, in contrast with 1.

In the rest of this section, we describe our main algorithm (4.1.1) and prove the above two key results (in [sec:analysis,subsubsec:uniFails], respectively).

4.1.1 Algorithm↩︎

Consider the setting of 1. For each element \(e \in V\) and cluster \(C_j \in \mathcal{C}\), let \[\label{eq:weights} W_{C_{j}}(e) = \sum_{e' \in C_j}d(e, e')\enspace.\tag{4}\] Given the definition of \(W_{C_{j}}(e)\) in 4 , for an element \(e\) of a cluster \(C\) we can express the quantities \(a(e)\) and \(b(e)\) in the definition of the silhouette \(s(e)\) (from Equation 2 ) as \[a(e) = \frac{W_{C_{}}(e)}{|C|-1} \qquad \text{and} \qquad b(e) = \min_{C_j \neq C}\frac{W_{C_{j}}(e)}{|C_j|} \enspace.\]

Building on the above formulation, we rely on the accurate estimation of the terms \(W_{C_{j}}(e)\) to approximate the silhouette \(s(e)\) for each element \(e\in V\).

A first attempt to obtain an estimate of \(W_{C_{j}}(e)\) would be to sample uniformly at random a set of \(t\ge 1\) points from each cluster \(C_j\in \mathcal{C}\) to be used to estimate \(W_{C_{j}}(e)\). Unfortunately, in 4.1.3 we show an instance of \(V\) and \(\mathcal{C}\) where such an approach fails with large probability for the sample complexity guarantees of 1 over a large fraction of points (see 2). That is, with a uniform sampling approach we cannot obtain the accurate approximation of the silhouette values as captured in 1 with the guarantees provided by 1.

Therefore, to achieve an accurate estimation of the weights \(W_{C_{j}}(e)\), we exploit the more complex but effective Probability Proportional to Size (pps) sampling strategy proposed by [34]. The pps sampling strategy allows us to collect a suitable small sample \(F_{C_j}\) from each cluster \(C_j\). The sample \(F_{C_j}\) can be then used to approximate the value of \(W_{C_{j}}(e)\) of every element \(e \in V\), within user-defined error bounds and with high probability, using only weighted distances between \(e\) and elements in \(F_{C_j}\). The selection of the elements of \(F_{C_j}\) is not uniformly at random. Instead, it involves carefully designed sampling probabilities that favor the selection of distant elements from a suitable “central” element of each cluster \(C_j\) for \(j\in[k]\).

At high level, our algorithm silh-pps-all(2), consists of two phases: 1) For each cluster \(C_j \in\mathcal{C}\), the algorithm computes a sample \(F_{C_j}\) of expected size \(t= \ceil{(c\varepsilon^{-2})\ln{(5nk/\delta)}}\) for a suitably small constant \(c\); 2) each sample \(F_{C_j}\) is then used to estimate the weights \(W_{C_{j}}(e)\) for \(j\in [k]\), and such weights are combined to obtain \(\widehat{s}\,(e)\), i.e., the estimates of \(s(e)\) for each \(e\in V\). More in detail.

Phase 1. Each cluster \(C_j, j\in [k]\) is processed independently. If \(|C_j| \leq t\), then \(F_{C_j}\) is set to \(C_j\). Otherwise, we perform Poisson sampling over \(C_j\) where each element \(e \in C_j\) is included in \(F_{C_j}\) independently with a suitable probability \(p_e\). For \(e \in C_j\), probability \(p_e\) is determined as follows:

Phase 2. We now have access to a sample \(F_{C_j}\) obtained through Poisson sampling with individual probabilities \(p_e\) for each \(e\in V\) (computed over Phase 1), over each cluster \(C_j\). Hence, we use each \(F_{C_j}\) to obtain the estimates \[\widehat{W}_{C_{j}}(e) = \sum_{e' \in F_{C_j}}\frac{d(e, e')}{p_{e'}} \enspace,\] which are accurate estimators of \(W_{C_j}(e)\), as we prove in our analysis (Section 4.1.2). We then use values \(\widehat{W}_{C_{j}}(e), j\in[k], e\in V\) to compute \(\widehat{s}\,(e), e\in V\) as follows. Let \(e \in V\) with \(e\in C, C\in \mathcal{C}\), we compute estimates \(\widehat{a}(e) = \tfrac{\widehat{W}_{C_{}}(e)}{(|C|-1)}\) and \(\widehat{b}(e) = \min_{C_j \neq C} \left\{ \tfrac{\widehat{W}_{C_{j}}(e)}{|C_j|}\right\}\). Finally we obtain the estimates \(\widehat{s}\,(e)\) of the silhouette of each element \(e\in V\) as \[\label{eq:localSilhEst} \widehat{s}\,(e) = \frac{\widehat{b}(e)-\widehat{a}(e)}{\max\left\{\widehat{a}(e), \widehat{b}(e)\right\}} \enspace.\tag{5}\]

2 (silh-pps-all) reports all details of the resulting procedure that combines the phases described above.

4.1.2 Proof of 1↩︎

We prove 1 in two steps. First, we show that, with probability at least \(1 - 3\delta/5\), for each element \(e\in V\) the value \(\widehat{s}\,(e)\) computed by silh-pps-all approximates the actual value \(s(e)\) within a small absolute error expressed as a function \(f(\varepsilon)=\Theta(\varepsilon)\), with \(\varepsilon\in (0,1)\). Second, we upper bound the sample complexity and the number of distance computations with probability at least \(1 - 2\delta/5\).

Our key ingredient is the next lemma that provides a probabilistic upper bound on the relative approximation error of the estimates \(\widehat{W}_{C_{j}}(e)\) with respect to the actual values \(W_{C_{j}}(e)\) for an arbitrary element \(e\in V\), and \(j\in[k]\).

Let \(E\) be the event “the relative error guarantees stated in 1 hold for every element \(e \in V\) and every cluster \(C_j\)”. From now on our analysis conditions on the realization of event \(E\), even when not explicitly stated. Consider an arbitrary element \(e\in V\), and let \(\widehat{s}\,(e)\) be the estimate of the silhouette \(s(e)\) computed by silh-pps-all. The following key technical lemma, establishes a bound on the absolute error of the estimates for each element \(e\in V\). That is, we show \(|\widehat{s}\,(e) - s(e)| \le f(\varepsilon)\) as follows.

We observe that 2 provides rigorous probabilistic guarantees (holding with probability at least \(1-3\delta/5\)) on the approximation error of each element in the set \(\widehat{s}\,(V)\) returned by silh-pps-all, holding simultaneously for all the elements in \(V\). Note that, for example, when \(\varepsilon\le 1/2\) then \(4\varepsilon/(1-\varepsilon) \le 8 \varepsilon = \mathcal{O}\!\left(\varepsilon\right)\).

To complete the proof of 1, we now bound the sample complexity and the number of distance computations performed by 2, as follows.

Observe that the number of distance computations performed by silh-pps-all is bounded by \[\label{eq:distanceComp} \useshortskip \sum_{j=1}^k |C_j| (|F^0_{C_j}|+|F_{C_j}|)\enspace,\tag{7}\] where for ease of notation \(|F^0_{C_j}| = 0\) when \(|C_j| \leq t\). By combining 7 , together with 3, we obtain the bound on the sample size and the bound on the number of distance computations claimed in 1, which hold with probability at least \(1-2\delta/5\)—yielding 1 by [lem:mainerror,lem:singlesilh].

4.1.3 Proof of 2↩︎

We now show that simple Poisson sampling over each cluster \(C_j, j\in[k]\), with an expected sample size equal to silh-pps-all, is not sufficient to obtain an accurate estimate of weights \(W_{C_{j}}(e), j\in [k], e\in V\). That is, we will provide an instance of 1 where a uniform sampling technique with a sample complexity of \(\Theta(tk)\), i.e., roughly the same sample size of silh-pps-all, outputs estimates of \(s(e), e\in V\) with constant error \(\Omega(1)\) for a large fraction of elements in \(V\), proving 2. Our result shows that non-uniform sampling (e.g., pps sampling) is necessary to get the rigorous approximation guarantees provided by our algorithm silh-pps-all.

Building the instance. Let \(V=\{e_1,\dots,e_{m+1},e_1',\dots,e_{m}',e_{m+1}' \}\), \(|V|=2m+2\) for fixed and very large \(m\) and \(\mathcal{C}=\{C_1,C_2\}\) such that \(C_1= \{e_1,\dots, e_{m+1}\}, C_2= \{e_1',\dots, e_{m+1}'\}\) and fix \(d\) to be the Euclidean distance over \(\mathbb{R}\supset V\).⁵ For ease of presentation, and without loss of generality, we assume that the elements of each cluster are ordered according to their indexing, that is \(e_{i}< e_j\) (resp. \(e_i'> e_j'\)) if it holds \(1\le i< j\le m+1\). The elements of \(C_1\) are placed as follows: (\(i\)) the \(m\) elements \(e_i, 1\le i \le m\), span an interval of length \(2r\), \(r=o(1/m), r >0\) centered at the origin (i.e., \(d(e_i,0) \le r\)); (\(ii\)) \(e_{m_1+1} \ge 0\) is such that \(d(e_{m_1+1},0)\ge R_1\), for some large value \(R_1\) (that we will specify later). The structure of \(C_2\) is similar to the structure of \(C_1\): the \(m\) elements \(e_i', 1\le i \le m\) are within an interval of length \(r\) centered at \(-1/m\), and \(e_{m+1}'<0\) at distance \(R_2>0, R_2=\Theta(m)\) from \(0\). Furthermore, let the points \(e_i\) (resp. \(e_i'\)) for \(i\in [m]\) be equally spaced over the intervals containing them, that is \(d(e_i,e_{i+1})\) (resp. \(d(e_i',e_{i+1}')\)) equals \(\beta r\) with \(\beta \doteq 2/(m-1)\), for \(i\in [m-1]\). We display the geometry of the resulting dataset in 3. Finally let \(R_1\ge L'(R_2+1)\) for arbitrary large \(L'\).

Comparing uniform sampling with our algorithm silh-pps-all. Consider now a uniform Poisson sampling strategy (denoted as uniform sampling for short), that estimates \(W_{C_{j}}(e), j=1,2\) as follows:

It is easy to show that the estimators \(\widehat{W}_{C_{j}}^{\mathtt{U}}(e)\) are unbiased. We can also use the estimators \(\widehat{W}_{C_{j}}^{\mathtt{U}}(e)\) to approximate the silhouette \(s(e)\) in an analogous way as done in 5 . That is for each \(e\in C_1\), we consider the estimate \[\widehat{s}_{\mathtt{U}}(e) = \frac{\widehat{W}_{C_{2}}^{\mathtt{U}}(e)/(m+1) - \widehat{W}_{C_{1}}^{\mathtt{U}}(e)/m}{\max\left\{\widehat{W}_{C_{2}}^{\mathtt{U}}(e)/(m+1),\widehat{W}_{C_{1}}^{\mathtt{U}}(e)/m\right\}} \enspace.\] In the next lemma, we show that for a fixed sample size \(t=o(m)\) and with arbitrary large probability close to 1, it holds that \(|\widehat{s}_{\mathtt{U}}(e)-s(e)|\approx 2\). That is, the estimation error of the uniform sampling strategy corresponds to the maximum error attainable, for all \(e\in C_1\setminus\{e_{m+1}\}\). In contrast, the estimates provided by silh-pps-all, using the same value of \(t\) achieve, for \(e\in C_1\), \(|\widehat{s}\,(e)-s(e)|\approx 0\) with large probability, that can be made arbitrarily close to 1. A similar argument can be used to show a large constant estimation error for \(\widehat{s}_{\mathtt{U}}(e)\) over all elements in \(e\in C_2\) .

Proving our result. The final claim of 2 follows as a corollary of 2, and noting that \(|V|=n=2m+2\) and \(|C_1| = m+1\). We remark that our result provides evidence of the superior accuracy achieved by silh-pps-all over a naive uniform Poisson sampling approach. That is, our silh-pps-all provably reports much more accurate estimates of the values of \(s(e), e\in V\) that cannot be matched by a uniform sampling strategy, used to approximate the values of \(W_{C_{j}}(e),j\in[k]\), when both approaches retain the same sample complexity. We further complement our theoretical insights with our extensive experimental evaluation in 5.

4.2 Estimating the silhouette of a clustering↩︎

We now study the problem of obtaining a tight and accurate estimate of the silhouette of a clustering \(\mathcal{C}\), as captured in our next statement.

Recall that \(s(\mathcal{C})\) is largely used in practice for evaluating the quality of a clustering \(\mathcal{C}\)—where \(\mathcal{C}\) is often obtained optimizing popular clustering metrics such as the \(k\)-means or \(k\)-medoids objective functions [50]. Hence, obtaining an accurate approximation of \(s(\mathcal{C})\) is of high practical utility, e.g., for identifying a good value for the parameter \(k\) on large datasets where \(s(\mathcal{C})\) cannot be computed exactly.

In the next sections we will obtain several estimators for the value \(s(\mathcal{C})\), solving 2, where, in analogy with 1, we want \(f(\varepsilon) = \Theta(\varepsilon)\). In addition, we want to output our highly accurate estimate \(s'(\mathcal{C})\) performing significantly less distance computations required by the exact algorithm, i.e., \(\Theta(n^2)\).

4.2.1 Estimator based on the silhouette definition↩︎

Our first estimator (global-sampler, gl-s) is based on the definition of \(s(\mathcal{C})\) itself (see Equation 3 ). That is, \(s(\mathcal{C})\) corresponds to the average of the silhouette of all the elements \(e\in V\), where each term \(s(e)\) contributing to the average takes a value ranging from \(-1\) to \(1\). Therefore, our first estimator is based on sampling a sufficiently small number of elements from \(V\), and using the average of their exact silhouette value (i.e., \(s(e)\)) as proxy for the overall silhouette \(s(\mathcal{C})\).

The estimator is computed as follows: we sample \(\left\{ e_1,\dots,e_m \right\}\) random elements from \(V\) where each element \(e_\ell, \ell\in[m]\) is sampled uniformly at random from the set \(V\). The resulting estimator can be expressed as \[\label{eq:est1} \widehat{s}_1(\mathcal{C}) = \frac{1}{m}\sum_{\ell=1}^{m} s(e_{\ell}) \enspace.\tag{8}\]

The next result quantifies the guarantees offered by estimator \(\widehat{s}_1(\mathcal{C})\).

The time complexity (i.e., number of distance computations) required to compute \(\widehat{s}_1(\mathcal{C})\) is \(\mathcal{O}(mn)\), for fixed \(m\).

4.2.2 Estimator using silh-pps-all↩︎

Recall that in Section 4.1 we developed our algorithm silh-pps-all that given a dataset \(V\), a clustering \(\mathcal{C}\) of \(V\) and two parameters \(\varepsilon, \delta\in (0,1)\) computes a set \(\widehat{s}\,(V)\) such that for each \(e\in V\) it holds \(|\widehat{s}\,(e)-s(e)|\le (4\varepsilon)/(1-\varepsilon)\) with probability at least \(1-\delta\). Our second estimator (global-pps-full-sampler, gl-pps-f) approximates the silhouette \(s(\mathcal{C})\) using the estimates \(\widehat{s}\,(e)\) for each element \(e\in V\) as obtained by silh-pps-all, that is

\[\label{eq:est2} \widehat{s}_2(\mathcal{C}) = \frac{1}{n}\sum_{e \in V} \widehat{s}\,(e)\enspace,\tag{9}\] where the values \(\widehat{s}\,(e), e \in V\) are computed by silh-pps-all. The next theorem establishes an upper bound to the absolute error of \(\widehat{s}_2(\mathcal{C})\) with respect to \(s(\mathcal{C})\).

The proof is an immediate consequence of the definition of \(s(\mathcal{C})\) (3 ) and \(\widehat{s}_2(\mathcal{C})\) (9 ), and of the guarantees provided by [lem:mainerror,lem:singlesilh].

For a given (expected) sample size \(t\),⁶ the time complexity to compute \(\widehat{s}_2(\mathcal{C})\) is the same required by silh-pps-all, that is \(\mathcal{O}(nkt)\) with probability at least \(1-\delta\) (see 3, and the related discussion).

4.2.3 Estimator using two-phase sampling↩︎

Let us consider, for ease of notation, the set \(\widehat{s}\,(V) = \{\widehat{s}\,(e_1),\dots,\widehat{s}\,(e_n)\}\) of estimates computed with silh-pps-all. Note that the values in the set \(\widehat{s}\,(V)\) are fixed once the pps sampling phase of Algorithm 2 (i.e, Phase 1) has been performed, even if the values in \(\widehat{s}\,(V)\) have not explicitly been computed yet.⁷ We therefore approximate the estimator in Equation 9 as follows:

1. Execute Phase 1 of Algorithm 2 and obtain all samples \(F_{C_{j}}, j\in[k]\), so that the values of \(\widehat{s}\,(e_i), i\in[n]\) are determined (i.e., the set \(\widehat{s}\,(V)\)), even if not computed yet;

2. Select uniformly at random \(m\in\mathbb{N}, m\geq2\) indices \(i_1,\dots,i_m\) in the set \(\{1,\dots,n\}\);

3. For each \(i_j, j \in \left\{1,\dots,m\right\}\) evaluate \(\widehat{s}\,(e_{i_j})\);

4. Return the average \(\bar{s}(\mathcal{C})\) over the \(m\) estimates \(\widehat{s}\,(e_{i_j})\) obtained above.

The resulting estimator (global-pps-sub-sampler, gl-pps-s) can be expressed as follows, \[\label{eq:est3} \bar{s}(\mathcal{C})=\frac{1}{m}\sum_{j=1}^m \widehat{s}\,(e_{i_j})\enspace.\tag{10}\]

The following result shows that \(\bar{s}(\mathcal{C})\) is an unbiased estimator of \(\widehat{s}_2(\mathcal{C})\).

The following result provides the guarantees on the approximation of \(\widehat{s}_2(\mathcal{C})\) provided by its estimator \(\bar{s}(\mathcal{C})\).

We now prove the resulting error obtained from approximating the silhouette \(s(\mathcal{C})\) with \(\bar{s}(\mathcal{C})\).

It now remains to control \(|s(\mathcal{C}) - \bar{s}(\mathcal{C})|\leq \varepsilon_1+\varepsilon_2\) to hold with probability at least \(1-\delta\) within user defined accuracy and confidence \(\varepsilon, \delta\in(0,1)\). Since the final error is bounded by \(\varepsilon_1+\varepsilon_2\) we can set \(\varepsilon_2 = \varepsilon - \varepsilon_1\). Noting that \(\varepsilon_1 \le 4\varepsilon'/(1-\varepsilon')\) it is sufficient to set \(\varepsilon'=\varepsilon/q, q> 4\). Finally to control the confidence we can set \(\delta_1=\delta/2\) and \(\delta_2 = \delta/(2-\delta)\), yielding to the desired confidence to be at least \(1-\delta\)—solving 2 as desired.

When performing Phase 1 of silh-pps-all for a fixed \(t\), the time complexity to compute \(\bar{s}(\mathcal{C})\) is at most \(\mathcal{O}(n\log(nk/\delta)+mkt)\). That is the sum of the time complexity of Phase 1 (i.e., \(\mathcal{O}(n\log(nk/\delta)\)) and the time complexity of running Phase 2 of silh-pps-all only for the \(m\) elements \(e_{i_j}, j\in\left\{ 1,\dots,m\right\}\) (i.e., \(\mathcal{O}(mkt)\)).

We provide a summary of our sampling techniques designed to solve 2 in 1.

4.3 Distributed algorithms↩︎

In this section, we present a distributed design of our silhouette approximation algorithms (i.e., both for local and global estimation) using the MapReduce framework [30]–[32]. At the end of the section, we will briefly discuss how to simply adapt our design to the popular MPC model [33], [51], [52].

4.3.1 Computational model↩︎

MapReduce adopts a functional programming model where an algorithm is specified as a sequence of rounds. In each round:

Throughout the description of our algorithms, a key-value pair will be denoted as \((\it key \; | \; \it value_1, \ldots, \it value_r)\), where \(\langle \it value_1, \dots, \it value_r \rangle\) corresponds to an \(r\)-dimensional value, for (a small) constant \(r\). As customary in the MapReduce literature [31], we use the term reducer to refer to the application of a reduce function to a set of pairs with the same key, i.e., see step ([item:MRstep2]) in the above discussion.

When a MapReduce algorithm is executed in a distributed environment, the instances of the map (resp., reduce) functions of a round are assigned by the underlying system to the available computing nodes. That is, a MapReduce algorithm is agnostic to the actual number of computing nodes. Therefore, important performance indicators to be minimized for a MapReduce algorithm are:

In our setting, consider a \(k\)-clustering \(\mathcal{C} = \{C_1, \dots, C_k\}\) of a dataset \(V = \{e_0, \dots e_{n-1}\}\) of \(n\) elements, and assume that each element \(e\in V\) has additional information \(C(e)\) denoting the cluster to which \(e\in V\) is assigned to over \(\mathcal{C}\). Initially, we represent the input clustering by the following set of key-value pairs: \[\{ (i | e_i) : 0 \leq i < n \}\enspace.\] Reasonably, we assume that the values \(n\), \(k\), and \(t\) used in the MapReduce algorithm, and the size \(|C_j|\) for each \(j\in[k]\), are known globally. A key ingredient of our design is that in each round, the elements of \(V\), and the newly created intermediate data, will be partitioned into \(w\) subsets of size \(n/w\) each, where \(w \in [0,n-1]\) is a suitable design parameter. To simplify the description and the analysis of our methods, we assume that \(w \in \mathcal{O}\!\left(n^\alpha\right)\), for some (small) \(\alpha <1/2\). While our methods can be generalized to larger values of \(w\) (i.e., those corresponding to \(\alpha\geq 1/2\)) with minimal modifications, we highlight that \(w\) constitutes an upper bound to the maximum parallelism exploitable by the algorithm. Thus, assuming that \(w\in \mathcal{O}\!\left(\sqrt{n}\right)\) does not constitute a limitation in practice, especially for modern datasets of massive size \(n\).

4.3.2 Approximating the silhouette of each element↩︎

We now present how to adapt our algorithm silh-pps-all to the MapReduce framework. Our method requires 3 rounds: the first two rounds obtain the sample \(F_{C_j}\), for every cluster \(C_j\) and \(j\in[k]\), while the last round computes the estimates \(\widehat{s}\,(e)\) of \(s(e)\) of each element \(e\in V\).

Round 1.

Round 2.

Round 3.

Recall that 3 states that with probability at least \(1 -2\delta/5\), for every cluster \(C_j\) with \(|C_j| > t\), \(|F^0_{C_j}| = \mathcal{O}\!\left(\ln(k/\delta)\right)\) and \(|F_{C_j}| = \mathcal{O}\!\left(\varepsilon^{-2}\ln(nk/\delta)\right)\). Then, since \(t= \mathcal{O}\!\left(\varepsilon^{-2}\ln(nk/\delta)\right)\), it immediately follows that, with probability at least \(1-2\delta/5\) the maximum local memory required by any map or reduce phase in any of the above rounds is bounded by \(M_L = \mathcal{O}\!\left(n/w + w + k\varepsilon^{-2}\ln(nk/\delta)\right) = \mathcal{O}\!\left(n/w + k\varepsilon^{-2}\ln(nk/\delta)\right)\) and the aggregate memory is bounded by \(M_A = \mathcal{O}\!\left(wM_L\right)\).

4.3.3 Approximating the silhouette of a clustering↩︎

In Section 4.2, we introduced three methods for computing an approximation \(\widehat{s}(\mathcal{C})\) to the silhouette of the entire clustering \(\mathcal{C}\), see 1. The two estimators \(\widehat{s}_2(\mathcal{C})\) (see 9 ) and \(\bar{s}(\mathcal{C})\) (see 10 ), can be computed through simple aggregations of all, or an \(m\)-sample, of the individual silhouettes estimates \(\widehat{s}\,(e)\) for \(e\in V\). These two estimators can be computed in MapReduce through a straightforward adaptation of the 3-round algorithm described above, combining, in the reduce phase of the third round, the individual silhouette estimates (or the sample ones) within each of the \(w\) partitions, and then combining the resulting values in an extra round. Overall, this 4-round approach would use the same local and aggregate space bounds obtained in 4.3.2. Finally, concerning the sampling method that computes \(\widehat{s}_1(\mathcal{C})\) based on 8 , it can be simply implemented through the MapReduce framework, we omit the details for brevity, since the resulting approach is analogous to the techniques already introduced. Such final approach uses \(\mathcal{O}\!\left(1+(\log w / \log(n/w))\right) = \mathcal{O}\!\left(1\right)\) rounds, local memory \(M_L = \mathcal{O}\!\left(n/w + k\varepsilon^{-2}\ln(1/\delta)\right)\) and aggregate memory \(M_A=\mathcal{O}\!\left(wM_L\right)\),

4.3.4 Porting MapReduce algorithms to the MPC framework↩︎

The MPC framework models a parallel system with \(p\) machines, each equipped with a local memory of size \(s\). Similar to MapReduce, the computation is structured as a sequence of synchronous rounds. In a round, each machine

The MPC framework allows arbitrary communication patterns, subject to the only constraint that each machine may send or receive at most \(s\) point-to-point messages in each round, with messages sent in a round processed by their recipients in the next round. Similarly to MapReduce, the ultimate goal is to design MPC algorithms that run in as few rounds as possible, using local memory (heavily) sublinear in the aggregate memory size. However, MPC and MapReduce differ in the way they express parallelism:

It is simple to observe that the MapReduce algorithms described above (for both local and global silhouette estimation) do not rely on the dynamic parallelism offered by the MapReduce model. That is, in every round, the algorithms always partition the data into \(w\) disjoint subsets, and the computation essentially involves:

As a consequence, our algorithms can be immediately adapted to the MPC framework by associating one machine per partition, setting \(p=w\) and \(s=M_L\). Finally, since in our setting \(M_L=\Omega\!\left(n^{1/2}\right)\), all the required broadcast and aggregation operations can be performed in constant rounds using standard MPC techniques [53]. Thus, the overall number of rounds remains constant.

The following theorem, whose proof is a direct consequence of all the arguments developed in this section, provides a quantitative characterization of the performance of our distributed strategies.

5 Experimental evaluation↩︎

In this section, we present the results of our extensive experimental evaluation. We start by describing the setup of our experiments (5.1), and then we investigate the following issues.

5.1 Setup↩︎

5.1.1 Methods and baselines↩︎

We now introduce all methods considered in our experiments. We distinguish between methods for global (2) and local estimation (1) of silhouette values.

Global estimates. We compared our algorithms with two state-of-the-art baseline algorithms. In addition, we also considered a variant of our methods coupled with uniform sampling (see 4.1.3). The list of methods is as follows.

• Algorithms gl-s, gl-pps-f, and gl-pps-s(from 1), which compute the estimators in [eq:est1,eq:est2,eq:est3], respectively, using pps sampling to estimate individual silhouette values.

• Two algorithms, dubbed gl-uni-f and gl-uni-s, which compute the estimators in [eq:est2,eq:est3], respectively, using uniform sampling (rather than pps sampling) to estimate all individual silhouettes (see 4.1.3).

• “FS” by [48], which, for each element \(e\in V\) computes the term \(a(e)\) according to its definition (2 ), while estimates the term \(b(e)\) as \(1/|C_{\mu_e}| \sum_{e'\in C_{\mu_e}} d(e',e)\), where \(C_{\mu_e}\) corresponds to the cluster with closest center (or centroid) to \(e\in V\), different from the cluster \(C\) of \(e\). Observe that the estimate of \(b(e)\) is simplified with respect to the one proposed in the original paper. That is, we do not verify if the computed estimate is smaller than all the distances between \(e\) and all centers (or centroids) of the clusters in \(\mathcal{C}\setminus\{ C_{\mu_e} \}\)—to avoid performing an exact computation of \(b(e)\) when the check fails.

• “SIMPL” by [45], the simplified silhouette approach, which evaluates the silhouette \(s(e)\) of each element \(e\in V\) using the distances of \(e\) to the closest centers (or centroids) as proxies of the average distances. Formally, let \(\mu_1,\dots,\mu_k\) the centers of clusters \(C_1,\dots,C_k \in \mathcal{C}\) then for an element \(e\) in some cluster \(C\) the method SIMPL evaluates \(a(e) = d(e, \mu)\) where \(\mu\) is the center of \(C\), and \(b(e) = \min_{C_j\neq C} d(e,\mu_j)\).

Note that both the methods FS and SIMPL are deterministic, and do not offer quantitative guarantees on their estimation error, in contrast with our approaches from 1, where the estimation error can be controlled as a function of the input parameters.

Local estimates. When solving 1 we will compare our method silh-pps-all with the uniform-sampling based approach detailed in 4.1.3, that we denote by uni. That is, we do not consider FS that scales poorly on large data (our main focus for 1), and SIMPL, which has very high estimation errors.

5.1.2 Datasets↩︎

A summary of the datasets used in our experiments, and their statistics is reported in 2. We provide more details of the processing of each dataset in 10.1. We use small-size datasets (those marked with S in 2) and CR to showcase applications of [probl:GlobEst,probl:LocalEst] (see 5.5). As for medium-size datasets, they are mainly used in 5.2 and 5.3, and are processed as follows. We cluster each dataset using both the \(k\)-means and the \(k\)-medoids objectives. For the \(k\)-means objective we consider the Euclidean distance \(d_{\mathtt{eucl}}(x,y) = \sqrt{\sum_{i} (x_i-y_i)^2}\). For the \(k\)-medoids objective, in addition to \(d_{\mathtt{eucl}}\), we also consider several widely used distance functions:

• \(d_{\mathtt{cos}}\), the cosine distance \(d_{\mathtt{cos}}(x,y) = 1- \sum_{i} (x_i-y_i)/(\sqrt{\sum_i x_i^2} \sqrt{\sum_i y_i^2})\);

• \(d_{\mathtt{manh}}\), the Manhattan distance \(d_{\mathtt{manh}}(x,y) = \sum_{i} |x_i-y_i|\);

• \(d_{\mathtt{can}}\), the Canberra distance \(d_{\mathtt{can}}(x,y) = \sum_i |x_i-y_i|/(|x_i| + |y_i|)\enspace .\)

For each clustering objective (either \(k\)-means or \(k\)-medoids) and each distance, we use different values of the parameter \(k\in \{2,5,10,15,20\}\) to obtain different clusterings \(\mathcal{C}\). Hence, each distinct configuration (dataset, value of \(k\), distance, objective) leads to a different clustering \(\mathcal{C}\), which we use as an input to either 2 or 1.

Finally, for the large datasets from 2 we proceeded similarly to the medium-sized datasets, but we only clustered such data using the \(k\)-means objective to obtain different clusterings \(\mathcal{C}\) varying \(k\), with \(k\in\{2,5,10,15,20\}\). We do not use \(k\)-medoids given its well-known poor scalability on large data.

5.1.3 Environment and parameters↩︎

We provide extensive details on the setup of the experiments in 10.2. The code to reproduce our results is available online.⁸

For our experiments, we set a time limit of one-hour for the execution of each method on a given configuration. Unless otherwise specified, each algorithm was executed for 10 independent runs and the results show the average over the runs. For our approximate algorithms, we directly set the parameter \(t\) controlling the expected sample size, rather than obtaining \(t\) as a function of the other parameters (see 2). When comparing the algorithms from 1 for 2, and their uni variant, we will set their parameters to have a comparable runtime: we fix \(t\) for gl-pps-f and gl-uni-f and we set the sample size \(m\) of gl-s to \(kt\). Unless otherwise stated we set \(\delta=0.01\).

5.1.4 Exact silhouette scores↩︎

For the small and medium sized datasets from 2 we were able to compute the exact silhouette scores on all the various configurations described in 5.1.2 within one hour of computation. In contrast on large data, we were not able to obtain the exact silhouette scores of any of the multiple configurations within one-hour time limit. To obtain a good proxy for the exact values of \(s(e)\) for each element \(e\in V\), we then proceed as follows. We perform five independent executions of silh-pps-all and its uniform variant (uni) from 4.1.3, obtaining respectively \(\widehat{s}\,(e)_i\) and \(\widehat{s}_{\mathtt{U}}(e)_{i}\) with \(i\in\{1,\dots, 5\}\) with a high sample size of \(t=800\). We then estimate \(s(e) = (1/10) \sum_i {(\widehat{s}\,(e)_i + \widehat{s}_{\mathtt{U}}(e)_i)}\) to reduce the bias in the final result.

5.2 Global silhouette approximation↩︎

In this section we evaluate the accuracy of the various methods and baselines described in 5.1.1. We compare all methods considering the average error \(|\widehat{s}\,(\mathcal{C}) - s(\mathcal{C})| \in [0,2]\) over ten independent runs (see 5.1.3 for details), where \(\widehat{s}\,(\mathcal{C})\) corresponds to the approximation provided by the method. In our discussion we also mention the relative approximation error that corresponds to \(|\widehat{s}\,(\mathcal{C}) - s(\mathcal{C})|/s(\mathcal{C})\). For medium-sized datasets we set the parameters of the various methods (see 3) as: \(t=16\) for both gl-pps-f and gl-uni-f, \(m=\ceil{n/20}, t=16\) for gl-pps-s and gl-uni-s, and the sample size of gl-s to \(m'=tk\) depending on \(\mathcal{C}\) (see 5.1.3). Our choice of a small value of \(t\) is motivated by the size of our datasets, while the choice of \(m\) enforces the computation of the estimator in 10 on a strictly sublinear number of elements in \(n\).

5.2.1 Medium datasets↩︎

We show some representative results in 4, where we report the trade-off between the estimation error, and the speed-up achieved over the exact computation of the silhouette coefficient by the various methods over ten runs.

We start by observing that the FS method has a very high runtime almost matching the exact algorithm. The only configuration where the speedup is more than \(3\times\), with respect to the exact approach, is the one relative to the CR dataset with \(d_{\mathtt{can}}\) and \(k=10\). Such behavior is not surprising since, for small \(k\), it is likely that some clusters have roughly linear size, hence FS will need a quadratic number of distance computations to compute the exact \(a(e)\)’s. This makes the FS method not practical on large or medium sized data, since it scales poorly. Interestingly, we note that even if the FS method is a heuristic approach (not offering guarantees), it often produces estimates with the smallest error among all methods evaluated. This is due to the fact that, for our datasets and configurations, the cluster minimizing the term \(b(e), e\in V\) often corresponds to the cluster with the closest center in \(\mathcal{C}\) to \(e\).

Next, we note that the SIMPL method is often the fastest approach among all methods we compared (except on the SHU data clustered with \(k\)-means and \(k=2\)). Unfortunately, SIMPL is also the method with the largest approximation errors, which cannot be controlled by any parameter. For example, on the SHU data clustered with \(k\)-medoids and distance \(d_{\mathtt{eucl}}\) the exact silhouette is \(s(\mathcal{C})=0.24\), while SIMPL outputs an estimate with an absolute error of 0.2, implying a relative error of more than 80%. Similarly, for the CR data clustered with \(k\)-medoids, \(d_{\mathtt{can}}\) and \(k=10\), where the relative approximation error of SIMPL is more than 105%. Such large errors hinder the use of the SIMPL method in practice.

Concerning the randomized sampling algorithms we observe that gl-s often achieves the smallest estimation error compared to gl-pps-f and gl-uni-f–under our setting of parameters all three algorithms perform the same number of distance computations in expectation. Interestingly, for some inputs such as CR(with \(d_{\mathtt{can}}\) and \(k=2\)), SHU(with \(d_{\mathtt{cos}}\) and \(k=2\)), or IOT(with \(d_{\mathtt{cos}}\), \(k=10\)), our new gl-pps-f and gl-pps-s methods achieve small approximation errors, better than gl-s, while being slightly less efficient. Comparing gl-pps-f and gl-uni-f, the former often achieves smaller approximation errors with respect to the latter (e.g., on the SHU data with \(k\)-means and \(k=2\), or on the IOT data on most configurations). Such a higher accuracy, comes at the expense of a higher runtime. Interestingly, we observe significant differences in runtime only on medium-sized data. This is probably caused by the fact that Phase 1 of gl-pps-f has a non-negligible runtime compared to Phase 2 on medium-sized data such as IOT. Finally, we observe that both methods gl-pps-s and gl-uni-s retain the approximation guarantees offered by gl-pps-f and gl-uni-f, respectively, while often improving slightly (for gl-pps-s) or significantly (for gl-uni-s) the runtime required to evaluate the estimates.

5.2.2 Large datasets↩︎

For the experiments on large datasets, we clustered each dataset with \(k\)-means and distance \(d_{\mathtt{eucl}}\). We did not consider baseline FS for its poor scalability and the SIMPL method for its poor quality estimates (see 5.2.1). In addition, we used highly accurate estimates of \(s(\mathcal{C})\) as ground truth, given the poor scalability of the exact computation (see 5.1.4).

The results are reported in 5. First, as expected, all algorithms exhibit increasing accuracy as the sample size \(t\) grows. For example, on dataset GOW with \(k=5\), the error with \(t=128\) for gl-pps-f is less than half the error obtained with \(t=32\). Next, regarding the trade-off between runtime and accuracy of gl-s, gl-pps-f and gl-uni-f, which perform the same expected number of distance computations, we observe several interesting trends. In almost all cases gl-s exhibited the highest accuracy and smallest variance, while gl-pps-f provided more accurate estimates with (often significantly) smaller variance compared to gl-uni-f.

For what concerns gl-pps-s and gl-uni-s, which compute the estimators of gl-pps-f and gl-uni-f respectively, using only a sample of size \(m\) of the total elements \(n\), we observe that a small sample (\(m=n/20\)) is often sufficient to retain the accuracy of gl-pps-f and gl-uni-f. In addition, both gl-pps-s and gl-uni-s achieve remarkable speedups of up to \(10\times\) over gl-pps-f, gl-uni-f and gl-s. As in the case of medium datasets, the runtime of gl-pps-s is slightly higher than the one of  gl-uni-s given the overhead introduced by Phase 1 of 2.

Summary for issue I1. Our experiments show that neither FS nor SIMPL are suitable for estimating the global silhouette of a clustering, especially for large datasets. While FS features good accuracy, it does not yield significant speedups over the exact algorithm. On the other hand, SIMPL is very fast but it incurs high approximation errors, which cannot be controlled by any parameter.

As for randomized methods, gl-s often achieves the best trade-off between accuracy and runtime compared with gl-pps-f and gl-uni-f. gl-pps-f features better accuracy compared to gl-uni-f, at the expense of a higher runtime on most datasets. gl-pps-s(resp., gl-uni-s) attains comparable accuracy to gl-pps-f(resp., gl-uni-f) even when processing a small number of elements \(m\) compared to \(n\)—yielding a significant speedup on large datasets.

Therefore, for large datasets, where the exact computation of \(s(\mathcal{C})\) is impractical even for a single and fixed value of \(k\), our methods as well as gl-s provide accurate estimates efficiently. For example, on the GOW data, which has more than 6 million elements, gl-pps-f or gl-pps-s with \(m=n/10\) and \(t=128\) output very accurate silhouette estimates for all values of \(k\in \{2,5,10,20\}\) in less than five minutes and one minute, respectively.

5.3 Local silhouette estimation↩︎

We now validate our algorithm silh-pps-all that estimates all values of \(s(e), e\in V\), so to deal with issue I2. For ease of notation we denote silh-pps-all with pps. As baseline we consider the uniform Poisson-sampling method, that we denote with uni, which computes \(\widehat{s}_{\mathtt{U}}(e), e\in V\) using uniform Poisson sampling over each cluster \(C\in \mathcal{C}\) (see 4.1.3) . For each dataset and value of \(k\) considered in 4.2, we evaluate the accuracy of the local estimates in two settings. First, we consider a bucketing of the values of silhouette \(s(e), e\in V\). That is, we assign each element \(e\in V\) to a bucket based on its silhouette value \(s(e)\) where each bucket \(B_i\), with \(i=1,\dots,40\), corresponds to elements \(e\) with \(s(e) \in [L_i, U_i)\), where \(L_i = -1 + (i-1)/20\) and \(U_i = -1 + i/20\).⁹ We use the buckets \(B_i\) to evaluate the accuracy of the local estimates over each individual bucket. Second, we also analyze the accuracy of the pps and uni approaches over each individual cluster \(C_1,\dots,C_k\) in \(\mathcal{C}\).

5.3.1 Accuracy over buckets↩︎

For each bucket \(B_i\) let \(E_i\subseteq V\) be the elements of \(V\) with \(s(e)\in B_i\). We computed, for each independent run, the maximum error \(\max_{e\in E_i}\{|\widehat{s}\,(e)-s(e)|\}\) (where \(\widehat{s}\,(e)\) is obtained either with pps or uni), and averaged the maximum error over ten independent runs, to obtain the average maximum error (\(\bar{\eta}_{B_i}\)) and its standard deviation, shown in 6. Note that \(\bar{\eta}_{B_i}\) is a very precise metric to capture the local estimation performances, given by the fact that it considers the maximum error over each bucket \(B_i\), which can be very large, i.e., at most 2.

From 6 we first observe that both pps and uni approaches achieve smaller maximum error (\(\bar{\eta}_{B_i}\)) for buckets corresponding to elements with higher \(s(e)\). Such behavior is likely caused by the fact that for elements with high value of \(s(e)\), the difference between \(b(e)-a(e)\) (see 2 ) is high, hence even rough estimates of terms \(W_{C_{j}}(e)\) yield a good approximation of \(s(e)\). The previous observation also explains the larger error over elements with silhouette \(s(e)\) close to 0. Comparing the average error \(\bar{\eta}_{B_i}\) we observe that, pps significantly reduces the error (and its standard deviation) over all buckets with respect to uni on most inputs. For example, on the IOT data with \(d_{\mathtt{can}}\) and \(k=5\) the estimates reported by pps have half of the maximum error achieved by uni. Similarly, on the GOW data for \(k=2, t=32\) the standard deviation of the error of pps is significantly smaller on each bucket compared to uni. We observe that on some configurations pps achieves similar performances to the uni approach (e.g., GOW for \(k=10\)). Nevertheless, we never observed pps reporting estimates with significantly smaller accuracy compared to uni, highlighting the high quality of the estimates of pps as captured by 1. We provide some additional results in 11.

5.3.2 Accuracy over clusters↩︎

Next we investigate the local approximation error of pps and uni across each cluster of \(\mathcal{C}\). That is, for each cluster we computed both the maximum error \(\max_{e\in C} \{|\widehat{s}\,(e)-s(e)|\}\) and the average error \((1/|C|) \sum_{e\in C} |\widehat{s}\,(e)-s(e)|\) (where \(\widehat{s}\,(e)\) is computed either with uni or pps).

In 7, we visualize the maximum and average errors for different input configurations, over ten independent runs. We observe the pps method reports estimates with comparable or significantly smaller average or maximum error compared to uni. As an example, on the MT dataset with \(k=2\), the average and maximum errors of pps for Cluster 1 are respectively bounded by 0.12 and 0.21, while, in contrast, those of uni are bounded by 0.17 and 0.28, respectively. Similarly, on the GOW dataset with \(k=10\), on all but cluster-index 7, our silh-pps-all always attains smaller average and maximum errors with respect to uni. We report additional results in 11.

Summary for issue I2. We investigated the accuracy for the extremely challenging task of computing the local silhouette value of all elements \(e\in V\) (i.e., 1). Our experiments show that our silh-pps-all algorithm provides highly-accurate estimates for all elements \(e\in V\), with significantly reduced variance compared to the uni baseline. That is, the estimates of pps are more accurate independently of the actual silhouette values (5.3.1), and independently of the cluster to which the elements belong (5.3.2), which are desirable properties for practical applications.

5.4 Distributed implementation↩︎

As discussed in 4.3, our silh-pps-all can be efficiently implemented in a distributed setting, such as the one modeled by the MapReduce/MPC frameworks, using a constant number of rounds. In this section, we assess the scalability of this distributed implementation of silh-pps-all, thus dealing with Issue I3. Since our experiments are performed on a single machine we rely on a simple parallel implementation retaining the steps in 4.3 (see 10.2 for details of the platform used, and 10.3 for more details on our parallelization strategy).

Our goal is to study the parallel speedup of silh-pps-all varying the number of threads in \(\{2,4,8,16\}\), varying the number of clusters \(k\in\{2,5,10,15,20\}\), and varying \(t\in \{512,1024\}\). We considered our largest datasets from 2. For each configuration we perform five runs, and recorded the speedup over the sequential implementation. That is, let \(T_i\) be the average time (over the five runs) to execute silh-pps-all with \(i\in \{1,2,4,8,16\}\) threads, i.e., \(T_1\) is associated with the sequential implementation. We will study the average speedup \(T_1/T_i\).

The results are reported in 8. We observe a remarkable speed-up over the sequential implementation of the parallel algorithm up to 8 threads on all \(k\) values (except for \(k=2\), see discussion below). More specifically, silh-pps-all achieves an almost linear speedup using up to 8 threads. The speedup becomes slightly smaller than linear for a larger number of threads, which is expected due to the increased overhead introduced by synchronization. Next, we observe that the parallel speedup is significantly greater for larger values of \(k\) and \(t\), since our implementation parallelizes the work over different clusters in \(\mathcal{C}\) (see 10.3). That is, for smaller values of \(k\) and \(t\) the degree of parallelism in the computation is not sufficient to compensate the overheads incurred by the multi-threaded implementation, i.e., synchronization. Finally, we note that a larger value of \(t\) leads to a higher speedup (e.g., on MT with \(t=1024\) the speedup with 16 threads is about 14 for \(k=20\)). Remarkably, when using \(t=1024\) on PH with \(k\ge 5\) and 16 threads, we obtain a speedup of at least \(10\times\) compared to the sequential execution of silh-pps-all. Notably, over such configuration, our parallel implementation of silh-pps-all computes extremely accurate estimates in less than 30 seconds (see 4 in 11, and recall that the PH dataset cannot be processed by the exact algorithm even within hours).

Summary for issue I3. A simple distributed implementation of silh-pps-all achieves remarkable speedups over its sequential version. The parallel speedup and the high accuracy achieved by silh-pps-all enable the estimation of the silhouette of all elements of massive datasets. That is, our method enables an extremely efficient processing of datasets that cannot be analyzed through current (sequential) methods.

5.5 Applications↩︎

5.5.1 Silhouette plot construction↩︎

We now show that the local estimates in output to silh-pps-all can be used to build a silhouette plot of a clustering \(\mathcal{C}\) that closely matches the exact one (e.g., see an example of silhouette plot in [subfig:silhPlot]). In this setting, we focus on BR, our smallest dataset, for ease of visualization,¹⁰ and for the fact that we can compute the exact values of \(s(e)\), for all \(e\in V\).

9 reports an example of two clusterings \(\mathcal{C}\), obtained respectively with \(k\)-medoids(top) and with \(k\)-means(bottom). In the top panel, we observe that the estimates provided by silh-pps-all are generally more accurate than those provided by uni, whose estimates are rather far from their actual values, especially for Cluster 2. In the bottom panel, we first note that both silh-pps-all and uni provide similar values for the estimate of the global silhouette coefficient \(s(\mathcal{C})\) (dashed line). Interestingly, the two methods yield rather different silhouette plots. Namely, the plot obtained with the estimates of silh-pps-all is much closer to the one obtained with the exact algorithm. For example, the estimates of silh-pps-all on Cluster 1 never exceed 0.75, in contrast with uni, whose estimates can be larger than 0.8 (note that the exact values for Cluster 1 are all less than \(0.75\)). Similarly, for Cluster 9, silh-pps-all assigns negative values to very few elements as opposed to uni.

5.5.2 Selection of the best value of \(k\)↩︎

Let \(V\) be a dataset and let \(K=\{k_1,\dots, k_{|K|}\}\) be a set containing different values of \(k\) used to cluster \(V\). The global silhouette is often used to select the value of \(k\in K\) that yields the best clustering for \(V\). Specifically, for each \(k\in K\) we obtain a clustering \(\mathcal{C}_k\) (e.g., using \(k\)-means) and select \(\mathcal{C}_{k^*}\) as the best clustering for \(V\), where \(k^* =\arg\max_{k\in K}s(\mathcal{C}_k)\) [54]. In this subsection we compare the effectiveness of gl-pps-s versus gl-uni-s in detecting the best value of \(k\) to be used for clustering. We focus on these two estimators, since, as argued before, they are accurate and their computation is scalable (see summary for Issue I1 in 5.2). Thus, they can be used to compute high-quality approximations of \(s(\mathcal{C}_k)\) for many different candidate values of \(k\).

We considered two datasets for which we can compute \(s(\mathcal{C}_k)\) exactly and efficiently, namely CR and WI. We clustered each dataset with \(k\)-means, with different values of \(k \in K\), with \(K=\{2,\dots, 15\}\), and for each resulting clustering \(\mathcal{C}_k\) we computed the exact global silhouette \(s(\mathcal{C}_k)\). Then, for each of the two algorithms gl-pps-s and gl-uni-s, both run with \(m= n/10\) and \(t=10\), we performed \(N=200\) independent runs for each distinct value of \(k\). Let \(\widehat{s}\,(\mathcal{C}_k)_{j}\) be the estimate of either gl-pps-s or gl-uni-s on dataset \(V\) clustered with value \(k\in [2,15]\) at iteration \(j\in [1,N]\). We let \(\boldsymbol{k}^{\mathsf{exact}}\) be the vector of the 14 values of \(k \in K\), sorted by non-increasing values of the corresponding global silhouette, that is, \(s(\mathcal{C}_{\boldsymbol{k}^{\mathsf{exact}}_1}) \ge s(\mathcal{C}_{\boldsymbol{k}^{\mathsf{exact}}_{2}}) \ge \cdots \ge s(\mathcal{C}_{\boldsymbol{k}^{\mathsf{exact}}_{14}})\). Similarly, we let \(\boldsymbol{k}^{\mathrm{\small pps}}(j)\) (resp. \(\boldsymbol{k}^{\mathrm{\small uni}}(j)\)) denote the vector of the 14 values of \(k\in K\), sorted by non-increasing values of the estimates of the corresponding global silhouette provided by gl-pps-s(resp. gl-uni-s) in the \(j\)-th run, for \(j \in [1,N]\). Note that, given a run \(j\in[1,N]\) the best value of \(k\in K\) returned by an approximate method corresponds to \(\boldsymbol{k}^{\mathrm{\small pps}}_1(j)\), or \(\boldsymbol{k}^{\mathrm{\small uni}}_1(j)\). To assess the accuracy of gl-pps-s and gl-uni-s for selecting a reasonable value of \(k\in K\) with respect to the exact approach, we now define a cumulative accuracy function of each approximate method. Given \(\ell \in [1,|K|]\), the cumulative accuracy evaluated at \(\ell\) corresponds to the fraction of runs in which the approximate method selects as best value \(k\in K\) a value in the top-\(\ell\) positions of \(\boldsymbol{k}^{\mathsf{exact}}\) (i.e., \(k\in \{\boldsymbol{k}^{\mathsf{exact}}_1,\dots,\boldsymbol{k}^{\mathsf{exact}}_\ell\}\)). More formally given the set of \(N\) sorted vectors \(\widehat{\boldsymbol{k}}(j)\), where \(\widehat{\boldsymbol{k}}(j)\) is either \(\boldsymbol{k}^{\mathrm{\small pps}}(j)\) or \(\boldsymbol{k}^{\mathrm{\small uni}}(j)\), \(j\in[1,N]\), the cumulative accuracy corresponds to, \[\text{Cumulative accuracy}(\ell, \widehat{\boldsymbol{k}}) \doteq \frac{1}{N} \sum_{j=1}^{N} \sum_{i=1}^{\ell} \mathbf{1}\left[\widehat{\boldsymbol{k}}_1(j) = \boldsymbol{k}^{\mathsf{exact}}_i\right] \enspace .\] Where \(\boldsymbol{1}[\cdot]\) corresponds to the indicator function. Observe that, as a function of \(\ell\), the cumulative accuracy eventually hits the value of \(1\) for \(\ell = |K|\). Clearly, a function approaching \(1\) very quickly indicates that the associated estimation method features high precision and a small variance in selecting a good value of \(k\), according to \(\boldsymbol{k}^{\mathsf{exact}}\).

Results are reported in 10, where we observe that the gl-pps-s method has significantly higher accuracy than gl-uni-s. For example, more than 97% of the 200 runs of gl-pps-s yield the selection as best value of \(k\) one of the top-3 values in \(\boldsymbol{k}^{\mathsf{exact}}\) on dataset CR, and one of the top-2 values in \(\boldsymbol{k}^{\mathsf{exact}}\) on the WI dataset. In contrast, gl-uni-s achieves smaller accuracy: less than 95% of the 200 runs correctly yield the selection of one the top-3 values of \(k\) in \(\boldsymbol{k}^{\mathsf{exact}}\) on the CR and less than 92% of the runs detect one of the top-2 values in \(\boldsymbol{k}^{\mathsf{exact}}\) on the WI dataset. We observe that the global accuracy is affected by how the exact silhouette values vary across the top-\(\ell\) elements of \(\boldsymbol{k}^{\mathsf{exact}}\) (we show the ordering of \(\boldsymbol{k}^{\mathsf{exact}}\) in 5). For example, in the above experiments the silhouette values associated with the top-2 values of \(\boldsymbol{k}^{\mathsf{exact}}\) do not differ much, making the detection of the best \(k\) extremely challenging, especially for gl-uni-s.

Summary for issue I4. We studied two different and widely adopted applications of the silhouette coefficient. First, we showed that the local estimates of silh-pps-all can be used to approximate the silhouette plots of a clustering \(\mathcal{C}\), matching exact plots. Second, we showed that our algorithm gl-pps-s for global silhouette estimation can be safely used in lieu of the expensive exact algorithm to detect the best clustering granularity (i.e., the value \(k\) yielding higher \(s(\mathcal{C}_k)\)) for a given dataset \(V\), improving over existing methods. Our experiments showcase the large applicability of the methods designed in this work for both local and global estimation scenarios.

6 Conclusions↩︎

In this work, we introduced an efficient, sampling-based algorithm to estimate the silhouette of all elements of a clustered dataset. Our new method (silh-pps-all), dramatically reduces the quadratic complexity required by the exact computation, while providing quantitative guarantees on the approximation error. In addition, we introduced several efficient estimators for the global silhouette providing a tight characterization of their required sample size.

Finally , we designed a distributed variant of our techniques, which requires small (sublinear) local memory, linear aggregate memory, and a constant number of computation rounds.

We validated our techniques through extensive experiments that show the accuracy, efficiency, and parallel scalability of our new methods for local and global silhouette estimation. For the local silhouette values, we showed that silh-pps-all provides highly accurate estimates compared to existing (heuristic) methods. For the global silhouette value, we showed that gl-pps-s(based on silh-pps-all’s sampling schema) often achieves the best trade-off between number of samples processed and accuracy. That is, gl-pps-s reports highly accurate estimates with small variance while being extremely efficient. We performed further experiments that showcase the use of our methods for two widespread applications:

These are some open issues deserving further investigation:

7 Acknowledgments↩︎

This research is funded by the Ministry of University and Research of Italy within the Complementary National Plan PNC-I.1 “Research initiatives for innovative technologies and pathways in the health and welfare sector, D.D. 931 of 06/06/2022, PNC0000002 DARE - Digital Lifelong Prevention CUP: B53C22006440001”. This research is also funded by the ERC Advanced Grant REBOUND (834862), and the Wallenberg AI, Autonomous Systems and Software Program (WASP) funded by the Knut and Alice Wallenberg Foundation.

8 Technical tools↩︎

The next results is a standard tool for data mining and machine learning algorithms, to provide concentration for sums of independent and bounded random variables [37], and is crucial for the analysis of our estimators.

The next results (belonging to the large farmily of bounds known as Chernoff-Hoeffding bounds) are used in the proof of 1. The form of the bounds that we use is from the textbook by [55].

9 Estimating other clustering validation measures↩︎

Our strategy that estimates the terms \(W_{C_{j}}(e), j\in[k], e\in V\) through the pps approach can be used to approximate various other measures for internal clustering evaluation. In particular, all those measures that are based on sums of distances between elements of the dataset \(V\). This is the case of measures that consider the cohesion (i.e., the average intracluster distance) and the separation (i.e., the average intercluster distance) of a clustering \(\mathcal{C}\) [15].

More precisely, consider a \(k\)-clustering \(\mathcal{C} = \{C_1, \dots, C_k\}\) and define the cohesion and separation of \(\mathcal{C}\) as \[\begin{align} Coh(\mathcal{C}) & = & \frac{1}{2} \frac{\sum_{j=1}^k \sum_{e',e'' \in C_j} d(e',e'')}{\sum_{j=1}^k {|C_j| \choose 2}} \enspace,\\ Sep(\mathcal{C}) & = & \frac{\sum_{1 \leq j_1 < j_2 \leq k}\sum_{e' \in C_{j_1}}\sum_{e'' \in C_{j_2}} d(e',e'')}{\sum_{1 \leq j_1 < j_2 \leq k} (|C_{j-1}||C_{j_2}|)} \enspace, \end{align}\] respectively. These measures have also been used to evaluate the average cluster reliability on networks, where distances correspond to connection probabilities [10], [56]. We can rewrite the above measures in terms of the sums \(W_{C_{}}(e)\) defined in 4 , as follows \[\begin{align} Coh(\mathcal{C}) = \frac{1}{2} \frac{\sum_{j=1}^k \sum_{e \in C_j} W_{C_j}(e)}{\sum_{j=1}^k {|C_j| \choose 2}} \quad \text{ and } \quad Sep(\mathcal{C}) = \frac{\sum_{1 \leq j_1 < j_2 \leq k}\sum_{e \in C_{j_1}} W_{C_{j_2}}(e)}{\sum_{1 \leq j_1 < j_2 \leq k} (|C_{j_1}||C_{j_2}|)} \enspace . \end{align}\]

Clearly, highly accurate approximations of terms \(W_{C_{j}}(e), j\in[k], e\in V\) with low relative errors yield approximations with low relative error for \(Coh(\mathcal{C})\) and \(Sep(\mathcal{C})\). Specifically, let \(\widehat{\rm Coh}(\mathcal{C})\) and \(\widehat{\rm Sep}(\mathcal{C})\) as the approximation of \(Coh(\mathcal{C})\) and \(Sep(\mathcal{C})\) respectively. Where each approximation is obtained by substituting the term \(W_{C}(e)\) with the value \(\widehat{W}_{C_{}}(e)\) computed within silh-pps-all. The following result is an immediate consequence of 1.

10 Reproducibility↩︎

10.1 Datasets↩︎

We consider various datasets from the UCI and SNAP repositories, that are publicly available online.¹¹ We provide a script to download and process each dataset with our code.¹² For the \(k\)-means objective we used the implementation from scikit-learn,¹³ which uses Lloyd’s algorithm coupled with an initialization obtained through \(k\)-means++. For the \(k\)-medoids objective dataset we used the state-of-the-art approach by [57].¹⁴

10.2 Setup↩︎

We implemented our algorithms and the baseline methods in C++20, and compiled it under gcc 9.5.0, with all optimization flags set. Except otherwise stated, all the experiments were performed single-threaded on a 72-core machine Intel Xeon Gold, running Ubuntu 20.04. To cluster the various datasets, we used Python 3.12.3. For the \(k\)-means objective, we used the default implementation from the Scikit-learn library,¹⁵ while for the \(k\)-medoids objective we used the \(k\)-medoids library [58].¹⁶ Our distributed implementation of silh-pps-all is done using the OpenMP library for C++, given that our execution is on a single many-core machine.

10.3 Parallel implementation↩︎

To parallelize our algorithm silh-pps-all we rely on the OpenMP library. Our implementation adopts the following parallelization strategy:

• We parallelize the for loops in [line:initSampleW,line:forLoopPPS] of 2, inspired by Round 1 and Round 2 of our Mapreduce algorithm;

• We parallelize the for loop in [line:forLoopAllElements] of 2, inspired by Round 3 of our Mapreduce algorithm.

For each of the above parallel loop executions, we perform synchronization on the necessary data-structures, yielding some of the bottlenecks discussed in 5.4.

11 Additional results↩︎

In this section we complement our extensive experimental evaluation, by presenting additional results to those of 5.

• In 11 we show, for all the missing large dataset of 2. The results of the experiments in 5.2.2 that compare the global estimates of the various randomized algorithms.

• In 12 we show additional results for the accuracy of pps and uni approaches over various buckets—extending results from 5.3.1.

• In 13 we show additional results on the average and maximum error on each cluster \(\mathcal{C}\) of pps and uni for the setting of 5.3.2.

• In 4 we show the average sequential runtime required by silh-pps-all considering the setting of 5.4.

• In 5 we show the values of \(s(\mathcal{C}_k)\) for each \(k\in \boldsymbol{k}\) under the setting of 5.5.

References↩︎

---

1. We will drop \(k\) when clear from the context.↩︎

2. Recall for \(\boldsymbol{x},\boldsymbol{y}\in \mathbb{R}^z\) the Euclidean distance \(d(\boldsymbol{x},\boldsymbol{y})\) is \(\sqrt{\sum_i (\boldsymbol{x}_i-\boldsymbol{y}_i)^2}\) where \(\boldsymbol{x}_i, i\in[z]\) denotes the \(i\)-th coordinate of \(\boldsymbol{x}\).↩︎

3. https://spark.apache.org/↩︎

4. A central element, is an element close to a majority of the elements of \(C_j\), see the proof of 1 for a formal definition.↩︎

5. Recall for \(a,b\in \mathbb{R}\) the Euclidean distance corresponds to \(d(a,b) = |a-b|\), denoted also as length.↩︎

6. Recall that the number of samples in each cluster in silh-pps-all is \(\mathcal{O}(t)\) with probability at least \(1-2\delta/5\).↩︎

7. After Phase 1, the values \(\widehat{s}\,(e_1),\dots,\widehat{s}\,(e_n)\) are determined, i.e., there is no randomness once the samples \(F_{C_{j}}\) for \(j\in[k]\) are fixed.↩︎

8. https://github.com/iliesarpe/ScalableSilhouetteComputation.↩︎

9. The bucket corresponding to index \(i=40\) also contains the right endpoint \(L_{40}=1\).↩︎

10. To avoid cluttered plots, it is good practice to visualize only a small subset of the elements of \(V\) over silhouette plots; even on large datasets.↩︎

11. https://archive.ics.uci.edu/datasets/, https://snap.stanford.edu/data/↩︎

12. https://github.com/iliesarpe/ScalableSilhouetteComputation↩︎

13. https://scikit-learn.org/↩︎

14. https://pypi.org/project/kmedoids/↩︎

15. https://scikit-learn.org/stable/↩︎

16. https://pypi.org/project/kmedoids/↩︎