Scalable and Distributed Silhouette Approximation

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

Scalable and Distributed Silhouette Approximation

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

cloud distributedcomputing parallelcomputing

distributed computing, parallel computing, cloud

arXiv:2607.01993v1 [cs.DS] 2 Jul 2026

Scalable and Distributed Silhouette Approximation Ilie Sarpe [email protected] KTH Royal Institute of Technology

Federico Altieri [email protected] University of Padova

Andrea Pietracaprina [email protected] University of Padova

Geppino Pucci [email protected] University of Padova

Fabio Vandin [email protected] University of Padova Abstract The silhouette is one of the most widely used measures to assess the quality of a 𝑘clustering of a dataset of 𝑛 elements. Its popularity stems from the fact that its evaluation requires no information beyond the clustering assignment. In addition, the silhouette is extremely easy to interpret, with a plethora of applications in various domains. The silhouette provides a score to measure the quality of a clustering as a whole or for each individual element. However, the exact computation of the: (𝑖) silhouette of each element of a dataset; and (𝑖𝑖) the global silhouette of the entire clustering; require Θ(𝑛 2 ) distance calculations, under general metrics. The quadratic complexity Θ(𝑛 2 ) is extremely prohibitive, especially on massive modern datasets. Surprisingly, existing approximate methods using O (𝑛 2 ) distance calculations are heuristics and do not offer provable and controllable guarantees on the quality of their results. In this work, we introduce the first rigorous and efficient algorithms to estimate: (𝑖) the (local) silhouette of each element of a dataset; and (𝑖𝑖) the (global) silhouette; of any metric 𝑘-clustering. Our methods are based on sampling, performing O (𝑛𝑘𝜀 −2 ln(𝑛𝑘/𝛿)) distance computations, and providing estimates with additive error O (𝜀) with probability at least 1 − 𝛿. That is, the user-defined parameters 𝜀 and 𝛿 in (0, 1) control the trade-off between accuracy and efficiency. Furthermore, we introduce a scalable and distributed design of our methods for the MapReduce and Massively Parallel Computing (MPC) frameworks. Our distributed algorithms require a constant number of rounds and sublinear local memory, under practical parameter setting. Finally, we perform extensive experiments to compare our methods with state-of-the-art approaches. The results show that our new techniques yield the best trade-off between accuracy and efficiency for both local and global silhouette estimation. In addition, our methods scale efficiently to massive datasets for which an exact computation of the silhouette is not practical.

Introduction

Clustering is of fundamental importance for data analysis, with ubiquitous applications in various areas including pattern recognition [Xu and Tian, 2015], bioinformatics and biomedicine [Ikotun et al., 2023], and data management [Aggarwal and Reddy, 2014, Karypis and Kumar, 1998]. Broadly, clustering a dataset requires grouping elements with high similarity, and separating dissimilar elements [Kleinberg, 2002]. 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 [Hennig et al., 2015]. In particular, recent research for clustering methods has its focus on the design of rigorous and scalable approaches to cluster massive datasets [Barger and Feldman, 2016, Ceccarello et al., 2017, 2019, Malkomes et al., 2015, Mazzetto et al., 2019, Pranjal and Balcan, 2013]. The vast number of clustering objectives and methods, introduce the need for rigorous methods to identify “good quality” clusterings of a dataset [Liu et al., 2010, Schubert, 2023]. 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 [Tan et al., 2005]. 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 [Hassan et al., 2024]. The silhouette is one of the most commonly used internal measures for clustering evaluation [Hassan et al., 2024, Rousseeuw, 1987]. The silhouette 𝑠 (𝑒) of an element 𝑒 assigned to a cluster 𝐶 is defined as the ratio (𝑏 (𝑒) − 𝑎(𝑒))/max{𝑎(𝑒), 𝑏 (𝑒)} where 𝑎(𝑒) is the average distance of 𝑒 from the other elements of 𝐶, and 𝑏 (𝑒) is the minimum average distance of 𝑒 from the elements of a cluster 𝐶 ′ different from 𝐶. In other words, 𝑠 (𝑒) provides a simple to interpret score between -1 and 1, which captures the quality of the assignment of 𝑒 to its group 𝐶. That is, if 𝑒 is closer (on average) to the elements in its cluster 𝐶 than to the elements in the “closest” cluster 𝐶 ′ ≠ 𝐶, then 𝑠 (𝑒) will be close to 1. In contrast, 𝑠 (𝑒) becomes negative when the element 𝑒 is closer to the elements in 𝐶 ′. As an example, consider elements 𝑒 13 and 𝑒 14 in Figure 1a (where the distance is the Euclidean distance): 𝑖) the silhouette 𝑠 (𝑒 13 ) is positive since 𝑒 13 is closer to the elements of its cluster (𝑒 11 and 𝑒 12 ) compared to its average distance from the elements of other clusters; 𝑖𝑖) the silhouette 𝑠 (𝑒 14 ) is close to -1 since 𝑒 14 is very close to all elements of cluster 𝐶 3 and distant from the other elements of its cluster 𝐶 1 . Hence, the silhouette 𝑠 (𝑒) of an element 𝑒 (or local silhouette) provides an interpretable and reliable score to assess the quality of the assignment of each element 𝑒 to its cluster [Dinh et al., 2019, Lovmar et al., 2005]. 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 [Hossain and Angryk, 2007, Ng and Han, 1994, Sellam et al., 2016, Wiwie et al., 2015], and is described in widely used data mining textbooks [Han et al., 2011, Tan et al., 2005]. As an example, the silhouette of a clustering can be used to identify a good value for the number 𝑘 of clusters computed by popular clustering algorithms, e.g., 𝑘-means, or 𝑘-medoids [Dinh et al., 2019, Lenssen and Schubert, 2024]. For general distance metrics, the exact computation of the silhouette for a 𝑘-clustering of 𝑛 2

𝑒 15 𝑒 13 𝑒 12 𝑒 11 𝑒8 𝑒7 𝑒9 𝑒6 𝑒5 𝑒 10 𝑒4 𝑒1 𝑒3 𝑒2

𝑒 10 𝑒3

𝑒7 𝑒2

𝑒6 𝑒8 𝑒9

𝑒1 𝑒5

𝑒4

𝑒 11

𝑒 12

𝑒 13

𝑒 15

−1

𝑒 14

𝑒 14 0 𝑠 (C) 𝑠 (𝑒)

(b)

(a)

Figure 1: (a): dataset 𝑉 = {𝑒 1 . . . , 𝑒𝑛 } for 𝑛 = 15 with elements in the Euclidean plane R2 . Different shapes represent 𝑘 = 4 different clusters, that is C = {𝐶 1 = {𝑒 1, . . . , 𝑒 3, 𝑒 14 }, 𝐶 2 = {𝑒 4, . . . , 𝑒 10 }, 𝐶 3 = {𝑒 11, 𝑒 12, 𝑒 13 }, 𝐶 4 = {𝑒 15 }}. (b): silhouette 𝑠 (𝑒) of the elements 𝑒 ∈ 𝑉 (where the distance is the Euclidean distance): values are grouped by clusters and sorted. The dashed line represents the value of the silhouette of the clustering C, marked with 𝑠 (C). Such representation is known as the silhouette plot [Rousseeuw, 1987]. elements requires Θ(𝑛 2 ) distance calculations, irrespective of 𝑘. 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 Θ(𝑛 2 ) barrier. In contrast, various scalable techniques have been proposed to efficiently cluster large datasets with rigorous guarantees, especially in distributed environments [Bahmani et al., 2012, Balcan et al., 2013, Ene et al., 2011, Feldman and Langberg, 2011, Malkomes et al., 2015]. In this work, we introduce the first randomized approximation algorithms with rigorous (probabilistic) guarantees on the quality of their estimates for two key problems: 𝑖) estimating the silhouette of each element 𝑒 of a given clustering; 𝑖𝑖) estimating the global silhouette of a clustering. Our algorithms feature a memory-efficient distributed implementation in popular frameworks such as MapReduce [Dean and Ghemawat, 2008, Leskovec et al., 2014, Pietracaprina et al., 2012], and the Massively Parallel Computing (MPC) model [Im et al., 2023], 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 𝜀 and 𝛿, controlling the trade-off between efficiency and accuracy. In particular, 𝜀 controls the (absolute) error between the output of our methods and the actual values being estimated, and 𝛿 controls the error probability (all our results hold with probability at least 1 − 𝛿). For the estimation of the silhouette of each element 𝑒, we rely on the Probability Proportional to Sampling (pps) technique by Chechik et al. [2015], which can be used to obtain highly accurate estimates of the average distance between 𝑒 and the elements of a given cluster. For a 𝑘-clustering, our method performs O (𝑛𝑘𝜀 −2 log(𝑛𝑘/𝛿)) distance computations, producing estimates b 𝑠 (𝑒) such that |𝑠 (𝑒) − b 𝑠 (𝑒)| ≤ 𝑓 (𝜀) with 𝑓 (𝜀) = O (𝜀) with probability at least 1 − 𝛿, simultaneously for all 3

elements 𝑒. Note the substantial improvement over the Θ(𝑛 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 𝑜 (𝑛 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: 𝑖) estimating the silhouette of all elements of a clustered dataset; 𝑖𝑖) 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 Θ(𝑛 2 ), reporting high-quality estimates with (absolute) error controlled by parameter 𝜀, with probability controlled by parameter 𝛿. • 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: 𝑖) 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; 𝑖𝑖) evaluate the accuracy of our method for the estimation of the local silhouette for all elements; 𝑖𝑖𝑖) assess the performance of our distributed methods; 𝑖𝑣) 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 [Altieri et al., 2021]. 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 [Altieri et al., 2021]. Our new methods require considerably fewer samples compared with the method appeared in [Altieri et al., 2021], yielding a reduction in the time complexity for global estimation. • We rigorously study the problem of local silhouette estimation, providing: 𝑖) a sound algorithm to solve the local silhouette estimation problem, based on techniques appeared in our previous version; 𝑖𝑖) rigorous bounds on the sample complexity, fixing existing issues with previous proofs in the literature [Chechik et al., 2015]; 𝑖𝑖𝑖) we prove that the 4

guarantees of our new method based on pps sampling cannot be matched by a simpler and more intuitive uniform sampling-based approach—supporting the design of our method. • We include all proofs omitted from [Altieri et al., 2021] 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 Section 2, we introduce the necessary preliminary notions. In Section 3, we discuss relevant previous work, with a focus on the silhouette approximation. In Section 4, we introduce our methods and their analysis. In Section 5, we present the results of our extensive experimental evaluation. Lastly, Section 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.

Preliminaries

Let 𝑈 be a metric space equipped with a distance function 𝑑 (·, ·), and let 𝑉 = {𝑒 1, . . . 𝑒𝑛 } ⊆ 𝑈 be a dataset of 𝑛 distinct elements from 𝑈 . Let C = {𝐶 1, . . . 𝐶𝑘 } be a 𝑘-clustering 1 of 𝑉 , that is C consist of a partition of 𝑉 into 𝑘 disjoint and non-empty subsets called clusters. Formally, the silhouette 𝑠 (𝑒) of an element 𝑒 ∈ 𝑉 belonging to a cluster 𝐶 ∈ C, as introduced by Rousseeuw [1987] is defined as, 𝑏 (𝑒) − 𝑎(𝑒) 𝑠 (𝑒) (1) max{𝑎(𝑒), 𝑏 (𝑒)}

where

𝑎(𝑒) =

′ 𝑒 ′ ∈𝐶 𝑑 (𝑒, 𝑒 )

|𝐶 | − 1

𝑏 (𝑒) = min

𝐶 𝑗 ∈C, 𝐶 𝑗 ≠𝐶

′ 𝑒 ′ ∈𝐶 𝑗 𝑑 (𝑒, 𝑒 )

|𝐶 𝑗 |

(2)

That is, 𝑎(𝑒) is the average distance of element 𝑒 ∈ 𝑉 from all other elements of its cluster 𝐶; while 𝑏 (𝑒) is the minimum average distance of 𝑒 ∈ 𝐶 from the elements in some cluster 𝐶 𝑗 ∈ C different from 𝐶. Note that by Equation (1) it holds that 𝑠 (𝑒) ∈ [−1, 1] ⊆ R. If element 𝑒 is the only element in its cluster 𝐶 such that 𝐶 = {𝑒}, we define 𝑠 (𝑒) = 0, as proposed by Rousseeuw [1987]. As an example, consider 𝑒 15 in Figure 1, then 𝑠 (𝑒 15 ) = 0. The silhouette 𝑠 (𝑒) of an element 𝑒 ∈ 𝐸 is a local measure that evaluates the quality of a clustering C with respect to the assignment of 𝑒 to its cluster 𝐶. Under a good assignment, the average distance of 𝑒 to the points in 𝐶 is much smaller than the average distance to the points of a different cluster, thus from Equation (1) it holds 𝑠 (𝑒) ≈ 1. Hence, high-values of 𝑠 (𝑒) denote a “good” assignment of element 𝑒 to its cluster 𝐶. Conversely, if 𝑒 is, on average, much 1We will drop 𝑘 when clear from the context.

closer to the elements of a different cluster 𝐶 𝑗 ≠ 𝐶 then by Equation (1) it holds 𝑠 (𝑒) ≈ −1, i.e., a “bad” cluster assignment for 𝑒. The individual 𝑠 (𝑒)’s are often represented through the silhouette plot [Rousseeuw, 1987], where the abscissae range in [−1, 1] (possible silhouette values) and the ordinate values correspond to the distinct elements, ordered by clusters. Example 1. Consider the dataset from Figure 1a and 𝑑 (·, ·) being the Euclidean distance.2 Figure 1b shows its silhouette plot, with points grouped by cluster. We observe that, for 𝑒 11 ∈ 𝐶 3 it holds 𝑠 (𝑒 11 ) > 0 since 𝑒 11 is much closer to 𝑒 12 and 𝑒 13 than to the points of any other cluster different from 𝐶 3 . In contrast, consider 𝑒 14 , then 𝑏 (𝑒 14 ) = (𝑑 (𝑒 11, 𝑒 14 ) + 𝑑 (𝑒 12, 𝑒 14 ) + 𝑑 (𝑒 13, 𝑒 14 ))/3 < (𝑑 (𝑒 1, 𝑒 14 ) + 𝑑 (𝑒 2, 𝑒 14 ) + 𝑑 (𝑒 3, 𝑒 14 ))/3 = 𝑎(𝑒 14 ). Therefore 𝑠 (𝑒 14 ) < 0, capturing the intuition that assigning 𝑒 14 to 𝐶 3 instead of 𝐶 1 leads to a better clustering. The silhouette of a clustering C evaluates the quality of a clustering C as the average silhouette of all the elements 𝑒 ∈ 𝑉 . Formally, the silhouette 𝑠 (C) of a clustering C is 1 ∑︁ 𝑠 (C) 𝑠 (𝑒𝑖 ) . 𝑛 𝑖=1 𝑛

(3)

Therefore, 𝑠 (C) provides a global measure of the quality of the whole clustering C, where a value close to 1 denotes a high quality of the clustering C. Example 2. Consider Figure 1b under the setting of Example 1: we observe that for the global silhouette of the clustering C it holds 𝑠 (C) ≈ 0.36. That is, on average, most points are well assigned over C, as we can observe from Figure 1a. Clearly, a positive value for 𝑠 (C), does not imply that all points have 𝑠 (𝑒) > 0, e.g., 𝑠 (𝑒 14 ) < 0.

In practice, the score 𝑠 (C) is widely used to identify a good value for 𝑘 (e.g., [Lletı et al., 2004]). That is, analysts often compare various clusterings {C𝑘𝑖 }𝑖≥1 obtained with different values of 𝑘𝑖 , selecting the best value of 𝑘 as the one maximizing 𝑠 (C𝑘𝑖 ), for 𝑖 ≥ 1. In fact, using the silhouette to select a good value of 𝑘 for a 𝑘-clustering is highly effective and reliable, in contrast to other popular methods, such as the elbow method, which often lead to poor results [Schubert, 2023]. It is clear from Equation (1) that the straightforward computation of 𝑠 (𝑒) for each element 𝑒 ∈ 𝑉 and the subsequent computation of 𝑠 (C) require Θ(𝑛 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 [Mitzenmacher and Upfal, 2017].

Related work

Since the introduction of the silhouette metric by Rousseeuw [1987], its use has been ubiquitous for clustering applications [Arbelaitz et al., 2013, Moulavi et al., 2014, Tomasini et al., 2016, Xiong 2 Recall for 𝒙, 𝒚 ∈ R𝑧 the Euclidean distance 𝑑 (𝒙, 𝒚) is

of 𝒙.

√︁Í

𝑖 (𝒙𝑖 − 𝒚𝑖 )

2 where 𝒙 , 𝑖 ∈ [𝑧] denotes the 𝑖-th coordinate 𝑖

and Li, 2013], including the use of the silhouette to obtain new clustering objectives for both 𝑘-medoids [Kaufmann and Rousseeuw, 1987, Lenssen and Schubert, 2024] and 𝑘-means [Lai et al., 2024]. However, as observed before, obtaining the exact silhouette is impractical for large datasets, due to the Θ(𝑛 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 𝑘-means clusterings). Van der Laan et al. [2003] introduced a variant of the silhouette 𝑠 (𝑒), where the term 𝑏 (𝑒) is computed as the distance from the closest center, other than the cluster 𝑒 belongs to (while 𝑎(𝑒) is computed according to Equation (2)). Following such idea, Hruschka et al. introduced the simplified silhouette. The simplified silhouette for an element 𝑒 assigned to cluster 𝐶 evaluates 𝑎(𝑒) as the distance between 𝑒 and the center of 𝐶, and 𝑏 (𝑒) as the distance between 𝑒 and the closest center of a cluster different from 𝐶. The complexity of computing the simplified silhouette reduces to O (𝑛𝑘). Hruschka et al. and Wang et al. [2017] show empirically that the simplified silhouette can be an effective evaluation measure for clusterings returned by Lloyd’s algorithm [Lloyd, 1982]. 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 Section 5, the difference between the values of the silhouette (Equation (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 Frahling and Sohler [2008]. For each element 𝑒 of a cluster 𝐶, the term 𝑎(𝑒) is computed according to its definition (Equation (2)), while for 𝑏 (𝑒) the heuristic first determines the average distance 𝑑¯𝑒,𝐶 ′ between 𝑒 and the elements of the cluster 𝐶 ′ ≠ 𝐶, whose centroid is closest to 𝑒. If, for any other cluster 𝐶 ′′ ∉ {𝐶, 𝐶 ′ }, it holds that the distance between 𝑒 and the centroid of 𝐶 ′′ is greater than or equal to 𝑑¯𝑒,𝐶 ′ , then it sets 𝑏 (𝑒) = 𝑑¯𝑒,𝐶 ′ , otherwise it computes 𝑏 (𝑒) according to Equation (2). For Euclidean distances, it can be shown that in the former case 𝑑¯𝑒,𝐶 ′ coincides exactly with the value of 𝑏 (𝑒) defined in Equation (2). The worst case complexity remains quadratic, but the heuristic likely performs less than Θ 𝑛 2 distance computations in practice. The Apache Spark programming framework provides optimized methods for computing the silhouette of a clustering under 𝑑-dimensional squared Euclidean distances and under one formulation of cosine distance.3 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 𝑘 clusters yields a parallelizable algorithm performing O (𝑛𝑘𝑑) 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 3 https://spark.apache.org/

clustering algorithms, distance query processing and centrality estimation in graphs [Balcan et al., 2013, Chechik et al., 2015, Cohen et al., 2018, Feldman and Langberg, 2011]. To the best of our knowledge, prior to the conference version of the current work [Altieri et al., 2021] the use of pps for efficient clustering evaluation had not been explored.

Methods

In this section we present our methods. In Section 4.1 we focus on the estimation of the silhouette of all individual elements, while in Section 4.2 we focus on the estimation of the silhouette of the whole clustering. Finally, Section 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, Problem 1. Given a dataset 𝑉 from an arbitrary metric space, a 𝑘-clustering C, and two parameters 𝜀, 𝛿 ∈ (0, 1), obtain b 𝑠 (𝑉 ) = {(𝑒,b 𝑠 (𝑒) : 𝑒 ∈ 𝑉 )} such that, with probability at least 1 − 𝛿, |b 𝑠 (𝑒) − 𝑠 (𝑒)| = O (𝜀) simultaneously for each element 𝑒 ∈ 𝑉 , minimizing the number of distance computations performed. Our main result, addressing Problem 1 is the following. Theorem 1. Consider a dataset 𝑉 , a 𝑘-clustering C, and two parameters 𝜀, 𝛿 ∈ (0, 1). There exists a randomized algorithm based on pps sampling with sample complexity O (𝑘𝜀 −2 log(𝑛𝑘/𝛿)) that, with probability 1 − 𝛿: 1. returns a set b 𝑠 (𝑉 ) = {(𝑒,b 𝑠 (𝑒) : 𝑒 ∈ 𝑉 )} such that, for all 𝑒 ∈ 𝑉 , |b 𝑠 (𝑒) − 𝑠 (𝑒)| ≤ 4𝜀/(1 − 𝜀); 2. performs O (𝑛𝑘𝜀 −2 log(𝑛𝑘/𝛿)) distance computations.

Theorem 1 shows that we can obtain a very accurate probabilistic approximation of each silhouette value 𝑠 (𝑒), 𝑒 ∈ 𝑉 . For practical values of 𝑘 ≪ 𝑛, and constant 𝜀 and 𝛿, the result of Theorem 1 breaks the quadratic barrier of Θ(𝑛 2 ) distance computations, enabling the efficient and practical computation of silhouette values over large and massive datasets. Note that the −2 O 𝜀 factor in the sample complexity and in the number of distance computations from Theorem 1, might become large for small 𝜀. While a quadratic dependence on 𝜀 cannot be avoided in a worst-case analysis [Mitzenmacher and Upfal, 2017], in practice, even large values of 𝜀 (e.g., close to 1) yield highly accurate estimates for all values in b 𝑠 (𝑉 ). 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 Theorem 1, a more intuitive uniform sampling approach performs arbitrarily poorly.

Theorem 2. There exists a large dataset 𝑉 , and a 2-clustering C of 𝑉 , such that the same algorithm as the one of Theorem 1, modified replacing pps sampling with uniform sampling, returns a set of estimatesb 𝑠 (𝑉 ) = {(𝑒,b 𝑠 (𝑒) : 𝑒 ∈ 𝑉 )}, with arbitrarily large probability, such that: |b 𝑠 (𝑒 ′)−𝑠 (𝑒 ′)| ≥ 1 ′ for at least 𝑛/2 − 1 elements 𝑒 ∈ 𝑉 . 8

Theorem 2 indicates that uniform sampling cannot guarantee an (arbitrarily small) estimation error O (𝜀) using a sublinear sample size, in contrast with Theorem 1. In the rest of this section, we describe our main algorithm (Section 4.1.1) and prove the above two key results (in Sections 4.1.2 and 4.1.3, respectively). 4.1.1

Algorithm

Consider the setting of Problem 1. For each element 𝑒 ∈ 𝑉 and cluster 𝐶 𝑗 ∈ C, let ∑︁ 𝑊𝐶 𝑗 (𝑒) = 𝑑 (𝑒, 𝑒 ′) .

(4)

𝑒 ′ ∈𝐶 𝑗

Given the definition of 𝑊𝐶 𝑗 (𝑒) in Equation (4), for an element 𝑒 of a cluster 𝐶 we can express the quantities 𝑎(𝑒) and 𝑏 (𝑒) in the definition of the silhouette 𝑠 (𝑒) (from Equation (2)) as 𝑎(𝑒) =

𝑊𝐶 (𝑒) |𝐶 | − 1

and

𝑏 (𝑒) = min 𝐶 𝑗 ≠𝐶

𝑊𝐶 𝑗 (𝑒) |𝐶 𝑗 |

Building on the above formulation, we rely on the accurate estimation of the terms 𝑊𝐶 𝑗 (𝑒) to approximate the silhouette 𝑠 (𝑒) for each element 𝑒 ∈ 𝑉 . A first attempt to obtain an estimate of 𝑊𝐶 𝑗 (𝑒) would be to sample uniformly at random a set of 𝑡 ≥ 1 points from each cluster 𝐶 𝑗 ∈ C to be used to estimate 𝑊𝐶 𝑗 (𝑒). Unfortunately, in Section 4.1.3 we show an instance of 𝑉 and C where such an approach fails with large probability for the sample complexity guarantees of Theorem 1 over a large fraction of points (see Theorem 2). That is, with a uniform sampling approach we cannot obtain the accurate approximation of the silhouette values as captured in Problem 1 with the guarantees provided by Theorem 1. Therefore, to achieve an accurate estimation of the weights 𝑊𝐶 𝑗 (𝑒), we exploit the more complex but effective Probability Proportional to Size (pps) sampling strategy proposed by Chechik et al. [2015]. The pps sampling strategy allows us to collect a suitable small sample 𝐹𝐶 𝑗 from each cluster 𝐶 𝑗 . The sample 𝐹𝐶 𝑗 can be then used to approximate the value of 𝑊𝐶 𝑗 (𝑒) of every element 𝑒 ∈ 𝑉 , within user-defined error bounds and with high probability, using only weighted distances between 𝑒 and elements in 𝐹𝐶 𝑗 . The selection of the elements of 𝐹𝐶 𝑗 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 𝐶 𝑗 for 𝑗 ∈ [𝑘]. At high level, our algorithm silh-pps-all (Algorithm 1), consists of two phases: 1) For each cluster 𝐶 𝑗 ∈ C, the algorithm computes a sample 𝐹𝐶 𝑗 of expected size 𝑡 = ⌈(𝑐𝜀 −2 ) ln (5𝑛𝑘/𝛿)⌉ for a suitably small constant 𝑐; 2) each sample 𝐹𝐶 𝑗 is then used to estimate the weights 𝑊𝐶 𝑗 (𝑒) for 𝑗 ∈ [𝑘], and such weights are combined to obtain b 𝑠 (𝑒), i.e., the estimates of 𝑠 (𝑒) for each 𝑒 ∈ 𝑉 . More in detail. Phase 1. Each cluster 𝐶 𝑗 , 𝑗 ∈ [𝑘] is processed independently. If |𝐶 𝑗 | ≤ 𝑡, then 𝐹𝐶 𝑗 is set to 𝐶 𝑗 . Otherwise, we perform Poisson sampling over 𝐶 𝑗 where each element 𝑒 ∈ 𝐶 𝑗 is included in 𝐹𝐶 𝑗 independently with a suitable probability 𝑝𝑒 . For 𝑒 ∈ 𝐶 𝑗 , probability 𝑝𝑒 is determined as follows:

Algorithm 1: silh-pps-all Input: Clustering C = {𝐶 1, . . . , 𝐶𝑘 } of 𝑉 = {𝑒 1, . . . , 𝑒𝑛 } and 𝜀, 𝛿 ∈ (0, 1), 𝑡 ≥ 2. Output: 𝑠ˆ(𝑉 ) = {(𝑒 1,b 𝑠 (𝑒 1 ), . . . , (𝑒𝑛 ,b 𝑠 (𝑒𝑛 ))}. ⊲ Phase 1: PPS sampling −2 1 𝑛 ← |𝑉 |; 𝑡 ← ⌈(𝑐𝜀 ) ln (5𝑛𝑘/𝛿)⌉; 2 for each cluster 𝐶 𝑗 ∈ C do 3 if |𝐶 𝑗 | ≤ 𝑡 then 𝐹𝐶 𝑗 ← 𝐶 𝑗 ; 4 else 5 𝐹𝐶0 𝑗 ← Poisson sampling of 𝐶 𝑗 with probability (2/|𝐶 𝑗 |) ln(5𝑘/𝛿); Í ¯ ← 𝑒 ′ ∈𝐶 𝑗 𝑑 (𝑒, ¯ 𝑒 ′); 6 for each 𝑒¯ ∈ 𝐹𝐶0 𝑗 do 𝑊𝐶 𝑗 (𝑒) 7 for each 𝑒 ∈ 𝐶 𝑗 do ¯ ¯ : 𝑒¯ ∈ 𝐹𝐶0 𝑗 }; 8 𝛾𝑒 ← max{𝑑 (𝑒, 𝑒)/𝑊 𝐶 𝑗 (𝑒) 9 𝛾𝑒 ← max{1/|𝐶 𝑗 |, 𝛾𝑒 }; 10 𝑝𝑒 ← min{1, 𝑡𝛾𝑒 }; 11

𝐹𝐶 𝑗 ← Poisson sampling of 𝐶 𝑗 with probabilities {𝑝𝑒 : 𝑒 ∈ 𝐶 𝑗 };

⊲ Phase 2: Silhouette estimation 12 for each 𝑒 ∈ 𝑉 do 13 Let 𝑒 belong to cluster 𝐶; 14 for each cluster 𝐶 𝑗 do b𝐶 𝑗 (𝑒) = Í𝑒 ′ ∈𝐹 𝑑 (𝑒, 𝑒 ′)/𝑝𝑒 ′ ; 15 𝑊 𝐶𝑗 16

17 18 19

b𝐶 (𝑒)/(|𝐶 | − 1); 𝑎b(𝑒) ← 𝑊 b𝐶 𝑗 (𝑒)/|𝐶 𝑗 | : 𝐶 𝑗 ≠ 𝐶}; 𝑏b(𝑒) ← min{𝑊 b 𝑠 (𝑒) ← (𝑏b(𝑒) − 𝑎b(𝑒))/(max{b 𝑎 (𝑒), 𝑏b(𝑒)});

return b 𝑠 (𝑉 ) = {(𝑒,b 𝑠 (𝑒)) : 𝑒 ∈ 𝑉 };

• First, we select an initial sample 𝐹𝐶0 𝑗 , by Poisson sampling. That is, each 𝑒 ∈ 𝐶 𝑗 is included in 𝐹𝐶0 𝑗 independently with fixed probability (2/|𝐶 𝑗 |) ln(5𝑘/𝛿). Each sample 𝐹𝐶0 𝑗 will contain, with sufficiently high probability, an aforementioned “central” element of 𝐶 𝑗 .4 ¯ for each sampled element 𝑒¯ ∈ 𝐹𝐶0 𝑗 . Then, for each 𝑒 ∈ 𝐶 𝑗 • We then compute the value 𝑊𝐶 𝑗 (𝑒) we set 𝑝𝑒 = min{1, 𝑡𝛾𝑒 }, where 𝛾𝑒 is computed as                   1  ¯  𝑑 (𝑒, 𝑒) 𝛾𝑒 = max , max .  |𝐶 𝑗 | 𝑒∈𝐹 ¯  ¯ 𝐶0 𝑊𝐶 𝑗 (𝑒)   𝑗    |{z} | {z }          𝑇1 𝑇 2  

4A central element, is an element close to a majority of the elements of 𝐶 , see the proof of Lemma 1 for a formal 𝑗 definition.

For 𝑒 ∈ 𝐶 𝑗 , 𝛾𝑒 is the maximum between two terms 𝑇1 and 𝑇2 . Where, 𝑇1 corresponds to the probability of selecting 𝑒 uniformly among all elements of 𝐶 𝑗 , ensuring that all points have a sufficiently large probability of being sampled; while 𝑇2 corresponds to the maximum ¯ of elements 𝑒¯ ∈ 𝐹𝐶0 𝑗 . The intuition is that relative contribution of 𝑒 to the weights 𝑊𝐶 𝑗 (𝑒) ¯ implies that element 𝑒 ∈ 𝐶 𝑗 is a higher relative contribution of element 𝑒 ∈ 𝐶 𝑗 to 𝑊𝐶 𝑗 (𝑒) ¯ with 𝑒¯ ∈ 𝐹𝐶0 𝑗 —captured by the extremely important for the accurate estimation of 𝑊𝐶 𝑗 (𝑒) value 𝛾𝑒 through 𝑇2 , which yields a high sampling probability 𝑝𝑒 for 𝑒 ∈ 𝐶 𝑗 . Phase 2. We now have access to a sample 𝐹𝐶 𝑗 obtained through Poisson sampling with individual probabilities 𝑝𝑒 for each 𝑒 ∈ 𝑉 (computed over Phase 1), over each cluster 𝐶 𝑗 . Hence, we use each 𝐹𝐶 𝑗 to obtain the estimates b𝐶 𝑗 (𝑒) = 𝑊

∑︁ 𝑑 (𝑒, 𝑒 ′)

𝑒 ′ ∈𝐹𝐶 𝑗

𝑝𝑒 ′

which are accurate estimators of 𝑊𝐶 𝑗 (𝑒), as we prove in our analysis (Section 4.1.2). We then use b𝐶 𝑗 (𝑒), 𝑗 ∈ [𝑘], 𝑒 ∈ 𝑉 to compute b values 𝑊 𝑠 (𝑒), 𝑒 ∈ 𝑉 as follows. Let 𝑒 ∈ 𝑉 with 𝑒 ∈ 𝐶, 𝐶 ∈ C, we b𝐶 (𝑒) 𝑊 b𝐶 (𝑒) 𝑊 𝑗 b compute estimates 𝑎b(𝑒) = (|𝐶 |−1) and 𝑏 (𝑒) = min𝐶 𝑗 ≠𝐶 |𝐶 𝑗 | . Finally we obtain the estimates

b 𝑠 (𝑒) of the silhouette of each element 𝑒 ∈ 𝑉 as b 𝑠 (𝑒) =

𝑏b(𝑒) − 𝑎b(𝑒) n o . b max 𝑎b(𝑒), 𝑏 (𝑒)

(5)

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

Proof of Theorem 1

We prove Theorem 1 in two steps. First, we show that, with probability at least 1 − 3𝛿/5, for each element 𝑒 ∈ 𝑉 the value b 𝑠 (𝑒) computed by silh-pps-all approximates the actual value 𝑠 (𝑒) within a small absolute error expressed as a function 𝑓 (𝜀) = Θ(𝜀), with 𝜀 ∈ (0, 1). Second, we upper bound the sample complexity and the number of distance computations with probability at least 1 − 2𝛿/5. Our key ingredient is the next lemma that provides a probabilistic upper bound on the relative b𝐶 𝑗 (𝑒) with respect to the actual values 𝑊𝐶 𝑗 (𝑒) for an approximation error of the estimates 𝑊 arbitrary element 𝑒 ∈ 𝑉 , and 𝑗 ∈ [𝑘]. Lemma 1. Let 𝑐 = 18. If 𝑡 = ⌈(3𝑐/𝜀 2 ) ln (5𝑛𝑘/𝛿)⌉, then with probability at least 1 − 3𝛿/5, for every b𝐶 𝑗 (𝑒) computed by silh-pps-all is such that element 𝑒 ∈ 𝑉 and every cluster 𝐶 𝑗 , the estimate 𝑊 b𝐶 𝑗 (𝑒) − 𝑊𝐶 𝑗 (𝑒) 𝑊 𝑊𝐶 𝑗 (𝑒)

≤𝜀 .

Proof. Our proof expands the poorly sketched argument in the proof of [Chechik et al., 2015, Lemma 12]. Consider an arbitrary cluster 𝐶 𝑗 . If |𝐶 𝑗 | ≤ 𝑡 the statement follows trivially, since b𝐶 𝑗 (𝑒) = 𝑊𝐶 𝑗 (𝑒). Thus, we focus on the case |𝐶 𝑗 | > 𝑡. For an element 𝑒 ∈ 𝐶 𝑗 , let 𝑚(𝑒) denote 𝑊 the median of the distances from 𝑒 to all other elements of 𝐶 𝑗 , that is, the distance between 𝑒 and the ⌈|𝐶 𝑗 |/2 + 1⌉-th closest element to 𝑒 in 𝐶 𝑗 . Let 𝑒 min = arg min𝑒 ′ ∈𝐶 𝑗 𝑚(𝑒 ′). The element 𝑒 is called well positioned if 𝑚(𝑒) ≤ 2𝑚(𝑒 min ). Using the triangle inequality, it is easy to see that the ⌈|𝐶 𝑗 |/2⌉ elements of 𝐶 𝑗 closest to 𝑒 min are well positioned. Hence, the initial random sample 𝐹𝐶0 𝑗 will contain a well positioned element with probability at least

|𝐶2𝑗 | 𝛿 2 5𝑘 ≥1− 1− 1− ln , |𝐶 𝑗 | 𝛿 5𝑘

since (1 + 𝑞/𝑥)𝑥 ≤ exp(𝑞), for 𝑥 ≥ 1, |𝑞| ≤ 𝑥, and |𝐶 𝑗 |/2 > 𝑡/2 > ln(5𝑘/𝛿). As proved in Chechik et al. [2015, Lemma 9], if 𝐹𝐶0 𝑗 contains a well positioned element, then for 𝑐 = 18 we have that for all 𝑒 ∈ 𝐶 𝑗 , the values 𝛾𝑒 computed by Algorithm 1 satisfy 𝛾𝑒 ≥

𝑑 (𝑒, 𝑒 ′) 1 max . 𝑐 𝑒 ′ ∈𝐶 𝑗 𝑊𝐶 𝑗 (𝑒 ′)

(6)

To obtain a worst-case bound, we assume that for each 𝑒 ∈ 𝐶 𝑗 it holds that 𝑝𝑒 = 𝑡𝛾𝑒 where b𝐶 𝑗 (𝑒 ′) of 𝑝𝑒 ≤ 1. In fact, when it holds 𝑝𝑒 = 1, the contribution of element 𝑒 to the weights 𝑊 b𝐶 𝑗 (𝑒 ′). all elements 𝑒 ′ ∈ 𝑉 is always computed exactly, implying no error over the estimates 𝑊 From Equation (6), we have that 𝑝𝑒 ≥ (𝑡/𝑐) max𝑒 ′ ∈𝐶 𝑗 (𝑑 (𝑒, 𝑒 ′)/𝑊𝐶 𝑗 (𝑒 ′)). Fix an arbitrary element 𝑒 ∈ 𝐶 𝑗 , and for every 𝑒 ′ ∈ 𝐶 𝑗 let 𝑋𝑒 (𝑒 ′) be the random variable taking value 𝑑 (𝑒, 𝑒 ′)/𝑝𝑒 ′ , with probability 𝑝𝑒 ′ , and 0 otherwise. We have that ∑︁ b𝐶 𝑗 (𝑒) = 𝑊 𝑋𝑒 (𝑒 ′) . 𝑒 ′ ∈𝐶 𝑗

b𝐶 𝑗 (𝑒)] = 𝑊𝐶 𝑗 (𝑒) by the linearity of expectation and the definiIt is immediate to see that E[𝑊 tion of variables 𝑋𝑒 (𝑒 ′). Let 𝜏𝑒 = 𝑊𝐶 𝑗 (𝑒)/(𝑡/𝑐) and, for every 𝑒 ′ ∈ 𝐶 𝑗 , let 𝑌𝑒 (𝑒 ′) = 𝑋𝑒 (𝑒 ′)/𝜏𝑒 be a normalized variable associated to each 𝑋𝑒 (𝑒 ′). Note that 𝑌𝑒 (𝑒 ′) ≤ (𝑑 (𝑒, 𝑒 ′)/𝑊𝐶 𝑗 (𝑒))/ (max𝑒 ′′ ∈𝐶 𝑗 {𝑑 (𝑒 ′′, 𝑒 ′)/𝑊𝐶 𝑗 (𝑒 ′′)}) ≤ 1 and therefore 𝑌𝑒 (𝑒 ′) ∈ [0, 1], hence  ∑︁    𝑊𝐶 𝑗 (𝑒) 𝑡 ′   𝜇𝑒 E  𝑌𝑒 (𝑒 )  = = . 𝜏𝑒 𝑐 𝑒 ′ ∈𝐶 𝑗   

The above discussion implies that

∑︁ 𝑡ª © ∑︁ ª b𝐶 𝑗 (𝑒) ≤ (1 + 𝜀)𝑊𝐶 𝑗 (𝑒) = P © P 𝑊 𝑌𝑒 (𝑒 ′) ≤ (1 + 𝜀) ® = P 𝑌𝑒 (𝑒 ′) ≤ (1 + 𝜀)𝜇𝑒 ® . 𝑐 ′ ′ ¬ «𝑒 ∈𝐶 𝑗 ¬ «𝑒 ∈𝐶 𝑗

Since variables 𝑌𝑒 (𝑒 ′)’s are independent random variables, and each 𝑌𝑒 (𝑒 ′) takes either value (1/𝜏𝑒 )𝑑 (𝑒, 𝑒 ′)/𝑝𝑒 ′ ≤ 1 with probability 𝑝𝑒 ′ , or 0 otherwise, by the application of Theorem 7 (see 12

Section A) we conclude that

2 𝛿 𝜀 𝑡 © ∑︁ ª ′ P ≤ , 𝑌𝑒 (𝑒 ) ≥ (1 + 𝜀)𝜇𝑒 ® ≤ exp − 3 𝑐 5𝑛𝑘 ′ «𝑒 ∈𝐶 𝑗 ¬ where the last inequality follows from the choice of 𝑡 as from the statement. Therefore, we have that b𝐶 𝑗 (𝑒) ≤ (1 + 𝜀)𝑊𝐶 𝑗 (𝑒) ≥ 1 − 𝛿 . P𝑊 5𝑛𝑘 A symmetrical argument can be used to also show that b𝐶 𝑗 (𝑒) ≥ (1 − 𝜀)𝑊𝐶 𝑗 (𝑒) ≥ 1 − 𝛿 . P𝑊 5𝑛𝑘 By the union bound, the probability that there exists a cluster 𝐶 𝑗 such that the initial sample 𝐹𝐶0 𝑗 does not contain a well positioned element is bounded by 𝑘𝛿/(5𝑘) = 𝛿/5. Now, by conditioning on the event that for all clusters 𝐶 𝑗 for 𝑗 ∈ [𝑘] the initial sample 𝐹𝐶0 𝑗 contains a well positioned element, we proceed with a second application of the union bound to prove our claim. That is, the probability that there exists an element 𝑒 ∈ 𝑉 and a cluster 𝐶 𝑗 for which b𝐶 𝑗 (𝑒) − 𝑊𝐶 𝑗 (𝑒))/𝑊𝐶 𝑗 (𝑒)| > 𝜀 is at most 2𝑛𝑘𝛿/(5𝑛𝑘) = 2𝛿/5. Therefore, the probability that |(𝑊 b𝐶 𝑗 (𝑒) − 𝑊𝐶 𝑗 (𝑒))/𝑊𝐶 𝑗 (𝑒)| ≤ 𝜀 for every 𝑒 ∈ 𝑉 is at least 1 − 𝛿/5 − 2𝛿/5 = 1 − 3𝛿/5, yielding |(𝑊 the statement and concluding the proof. □ Let 𝐸 be the event “the relative error guarantees stated in Lemma 1 hold for every element 𝑒 ∈ 𝑉 and every cluster 𝐶 𝑗 ”. From now on our analysis conditions on the realization of event 𝐸, even when not explicitly stated. Consider an arbitrary element 𝑒 ∈ 𝑉 , and let b 𝑠 (𝑒) be the estimate of the silhouette 𝑠 (𝑒) computed by silh-pps-all. The following key technical lemma, establishes a bound on the absolute error of the estimates for each element 𝑒 ∈ 𝑉 . That is, we show |b 𝑠 (𝑒) − 𝑠 (𝑒)| ≤ 𝑓 (𝜀) as follows. Lemma 2. Under event 𝐸, silh-pps-all outputs estimates b 𝑠 (𝑒) for each element 𝑒 ∈ 𝑉 such that, |b 𝑠 (𝑒) − 𝑠 (𝑒)| ≤

4𝜀 . 1−𝜀

Proof. First note that if 𝑎(𝑒) = 0 (that is, {𝑒} is a singleton cluster), then it follows that 𝑎b(𝑒) = 0 and 𝑠 (𝑒) = b 𝑠 (𝑒) = 0, and the claim trivially holds. We thus consider the case 𝑎(𝑒) > 0 (also, observe that 𝑏 (𝑒) > 0 by definition). When event 𝐸 holds, the relative errors of terms: (𝐼 ) |(b 𝑎 (𝑒) − 𝑎(𝑒))/𝑎(𝑒)|; and (𝐼𝐼 ) |(𝑏b(𝑒) − 𝑏 (𝑒))/𝑏 (𝑒)|, are both upper bounded by 𝜀 > 0. That is, the bound on term (𝐼 ) follows immediately from the definition of 𝑎b(𝑒) and the relative error b𝐶 𝑗 (𝑒) provided by Lemma 1. Concerning the bound on the error of term (𝐼𝐼 ), bound for 𝑊 b𝐶 𝑗 (𝑒)/|𝐶 𝑗 |. Similar consider an arbitrary cluster 𝐶 𝑗 and let 𝑏𝐶 𝑗 (𝑒) = 𝑊𝐶 𝑗 (𝑒)/|𝐶 𝑗 | and 𝑏b𝐶 𝑗 (𝑒) = 𝑊 to the argument for term (𝐼 ), it is easy to show that |(𝑏b𝐶 𝑗 (𝑒) − 𝑏𝐶 𝑗 (𝑒))/𝑏𝐶 𝑗 (𝑒)| ≤ 𝜀. Recall that 𝑏 (𝑒) = min𝐶 𝑗 ≠𝐶 𝑏𝐶 𝑗 (𝑒) and that 𝑏b(𝑒) = min𝐶 𝑗 ≠𝐶 𝑏b𝐶 𝑗 (𝑒). Suppose that 𝑏 (𝑒) = 𝑏𝐶 ′ (𝑒) and 𝑏b(𝑒) = 𝑏b𝐶 ′′ (𝑒), for some, possibly different, clusters 𝐶 ′ and 𝐶 ′′. We have: (1 − 𝜀)𝑏 (𝑒) = (1 − 𝜀)𝑏𝐶 ′ (𝑒) ≤ (1 − 𝜀)𝑏𝐶 ′′ (𝑒) ≤ 𝑏b𝐶 ′′ (𝑒) = 𝑏b(𝑒) 13

and

𝑏b(𝑒) = 𝑏b𝐶 ′′ (𝑒) ≤ 𝑏b𝐶 ′ (𝑒) ≤ (1 + 𝜀)𝑏𝐶 ′ (𝑒) = (1 + 𝜀)𝑏 (𝑒) ,

and the desired bound follows. Now we establish a bound on the relative error for the the denominator of the estimator b 𝑠 (𝑒) b (𝑒) max{b (see Equation (5)). Define 𝑀 (𝑒) max{𝑎(𝑒), 𝑏 (𝑒)} and 𝑀 𝑎 (𝑒), 𝑏b(𝑒)}. We show that b (𝑒) − 𝑀 (𝑒))/𝑀 (𝑒)| ≤ 𝜀. Suppose that 𝑀 (𝑒) = 𝑎(𝑒), hence 𝑎(𝑒) ≥ 𝑏 (𝑒) (the case 𝑀 (𝑒) = 𝑏 (𝑒) |( 𝑀 b (𝑒) = 𝑎b(𝑒), the bound is immediate. Instead, if 𝑀 b (𝑒) = 𝑏b(𝑒), hence 𝑏b(𝑒) ≥ 𝑎b(𝑒), is analogous). If 𝑀 the bound follows since b (𝑒) = 𝑏b(𝑒) ≤ (1 + 𝜀)𝑏 (𝑒) ≤ (1 + 𝜀)𝑎(𝑒) = (1 + 𝜀)𝑀 (𝑒) , 𝑀 and

b (𝑒) = 𝑏b(𝑒) ≥ 𝑎b(𝑒) ≥ (1 − 𝜀)𝑎(𝑒) = (1 − 𝜀)𝑀 (𝑒) . 𝑀

We are now ready to obtain a bound on the absolute difference |b 𝑠 (𝑒) − 𝑠 (𝑒))|. Using the relative error bounds established above, we have that 𝑏b(𝑒) − 𝑎b(𝑒) (1 + 𝜀)𝑏 (𝑒) − (1 − 𝜀)𝑎(𝑒) (1 + 𝜀)𝑏 (𝑒) − (1 − 𝜀)𝑎(𝑒) b 𝑠 (𝑒) = ≤ = . (1 − 𝜀)𝑀 (𝑒) (1 − 𝜀)𝑎(𝑒) b (𝑒) 𝑀 Simple algebraic manipulations show that

(1 + 𝜀)𝑏 (𝑒) − (1 − 𝜀)𝑎(𝑒) 2𝜀 4𝜀 (𝑠 (𝑒) + 1) ≤ 𝑠 (𝑒) + = 𝑠 (𝑒) + , (1 − 𝜀)𝑎(𝑒) 1−𝜀 1−𝜀

where the last inequality follows since 𝑠 (𝑒) ≤ 1. With analogous calculations, we obtain that 𝑏b(𝑒) − 𝑎b(𝑒) (1 − 𝜀)𝑏 (𝑒) − (1 + 𝜀)𝑎(𝑒) (1 − 𝜀)𝑏 (𝑒) − (1 + 𝜀)𝑎(𝑒) ≥ = (1 + 𝜀)𝑀 (𝑒) (1 + 𝜀)𝑎(𝑒) b (𝑒) 𝑀 4𝜀 4𝜀 2𝜀 (𝑠 (𝑒) + 1) ≥ 𝑠 (𝑒) − ≥ 𝑠 (𝑒) − , = 𝑠 (𝑒) − 1+𝜀 1+𝜀 1−𝜀

b 𝑠 (𝑒) =

which concludes the proof.

□

We observe that Lemma 2 provides rigorous probabilistic guarantees (holding with probability at least 1 − 3𝛿/5) on the approximation error of each element in the set b 𝑠 (𝑉 ) returned by silhpps-all, holding simultaneously for all the elements in 𝑉 . Note that, for example, when 𝜀 ≤ 1/2 then 4𝜀/(1 − 𝜀) ≤ 8𝜀 = O (𝜀). To complete the proof of Theorem 1, we now bound the sample complexity and the number of distance computations performed by Algorithm 1, as follows. Lemma 3. With probability at least 1 − 2𝛿/5, over the execution of silh-pps-all, for every cluster 𝐶 𝑗 with |𝐶 𝑗 | > 𝑡 it holds that, 𝑘 𝑛𝑘 0 −2 |𝐹𝐶 𝑗 | = O ln and |𝐹𝐶 𝑗 | = O 𝜀 ln . 𝛿 𝛿

Proof. Consider a cluster |𝐶 𝑗 | with |𝐶 𝑗 | > 𝑡. |𝐹𝐶0 𝑗 | corresponds to a binomial random variable Bin(|𝐶 𝑗 |, 2/|𝐶 𝑗 | ln(5𝑘/𝛿)), implying that E[|𝐹𝐶0 𝑗 |] = 2 ln(5𝑘/𝛿)). Therefore, a simple application of the Chernoff bound suffices to show that |𝐹𝐶0 𝑗 | = O (ln(𝑘/𝛿)) with probability at least 1−𝛿/(5𝑘), hence the bound holds for all 𝐹𝐶0 𝑗 ’s with probability at least 1 − 𝛿/5. Let us now bound the size of |𝐹𝐶 𝑗 | (our result expands on the argument of Chechik et al. [2015]). |𝐹𝐶 𝑗 | is a Poisson binomial random variable with expectation ∑︁ ∑︁ E |𝐹𝐶 𝑗 | = 𝑝𝑒 ≤ 𝑡 𝛾𝑒 . 𝑒∈𝐶 𝑗

Observe that ∑︁

(

∑︁

¯ 1 𝑑 (𝑒, 𝑒) 𝛾𝑒 = max , max ¯ |𝐶 𝑗 | 𝑒∈𝐹 ¯ 𝐶0 𝑊𝐶 𝑗 (𝑒) 𝑒∈𝐶 𝑗 𝑒∈𝐶 𝑗 𝑗 ∑︁ ¯ 𝑑 (𝑒, 𝑒) max ≤1+ 0 ¯ ¯ 𝐶 𝑊𝐶 𝑗 (𝑒) 𝑒∈𝐹 𝑒∈𝐶 𝑗 𝑗 ∑︁ ∑︁ 𝑑 (𝑒, 𝑒) ¯ ≤1+ ¯ 𝑊𝐶 𝑗 (𝑒) 0 𝑒∈𝐶

𝑒∈𝐶 𝑗

)

(by definition of 𝛾𝑒 , 𝑒 ∈ 𝑉 ) (for 𝐴, 𝐵 ≥ 0 then max{𝐴, 𝐵} ≤ 𝐴 + 𝐵) (max{𝐴𝑖 } ≤

𝑗 𝑒∈𝐹 ¯ 𝐶 𝑗

=1+

∑︁ ∑︁ 𝑑 (𝑒, 𝑒) ¯

¯ 𝐶0 𝑒∈𝐶 𝑗 𝑒∈𝐹 𝑗

∑︁

𝐴𝑖 , 𝐴𝑖 ≥ 0)

(swapping finite and bounded sums)

¯ 𝑊𝐶 𝑗 (𝑒)

= 1 + |𝐹𝐶0 𝑗 |

¯ . (definition of 𝑊𝐶 𝑗 (𝑒))

Now consider the following cases: (I) |𝐹𝐶0 𝑗 | = 1 consists of a single well-positioned element; (II) |𝐹𝐶0 𝑗 | = |𝐶 𝑗 |, i.e., the initial sample size consists of all elements of 𝐶 𝑗 ; (III) 1 < |𝐹𝐶0 𝑗 | < |𝐶 𝑗 |, Í accounting for the remaining cases. Clearly, by the above analysis, case (I) yields 𝑒 𝛾𝑒 ≤ 2. Now let us define 𝛾¯𝑒 = max{1/|𝐶 𝑗 |, max𝑒 ′ ∈𝐶 𝑗 𝑑 (𝑒, 𝑒 ′)/𝑊𝐶 𝑗 (𝑒)}, that is 𝛾¯𝑒 is computed through the actual weights 𝑊𝐶 𝑗 (𝑒) over all elements of 𝐶 𝑗 . For case (II) when |𝐹𝐶0 𝑗 | = |𝐶 𝑗 | the algorithm silh-pps-all relies on the values 𝛾¯𝑒 . Let 𝛾𝑒 , 𝑒 ∈ 𝐶 𝑗 be the coefficients obtained when 𝐹𝐶0 𝑗 contains a single well-positioned element, then by Equation (6) ∑︁ there exists a constant 𝑐 such that ∑︁ 𝛾¯𝑒 ≤ 𝑐 𝛾𝑒 ≤ 2𝑐 . 𝑒∈𝐶 𝑗

𝑒∈𝐶 𝑗

Í Therefore even when |𝐹𝐶0 𝑗 | = |𝐶 𝑗 | we have that 𝑒 𝛾𝑒 ≤ O (1). Note that case (II) implies that even in the highly unlikely case of |𝐹𝐶0 𝑗 | = |𝐶 𝑗 |, the final sample 𝐹𝐶 𝑗 will be such that E[|𝐹𝐶 𝑗 |] = O (𝑡). Finally, consider case (III), when 1 < |𝐹𝐶0 𝑗 | < |𝐶 𝑗 |, and let 𝛾𝑒 be the coefficients Í Í obtained by silh-pps-all with 𝐹𝐶0 𝑗 . By contradiction assume that 𝑒 𝛾𝑒 > 𝑒 𝛾¯𝑒 , for 𝛾¯𝑒 defined as n o Í 𝑑 (𝑒,¯ 𝑒) 𝑑 (𝑒,¯ 𝑒) 1 1 in case (II), which implies that 𝑒 max |𝐶 𝑗 | , max𝑒∈𝐹 > 0. ¯ 𝑗 𝑊𝐶 (𝑒) ¯ 0 𝑊𝐶 (𝑒) ¯ − max |𝐶 𝑗 | , max𝑒∈𝐶 ¯ 𝐶𝑗

𝑗

It is simple to observe that none of the terms in the summation can be larger than 0 since 𝐹𝐶0 𝑗 ⊂ 𝐶 𝑗 , yielding the contradiction. Hence, we have that for each set 𝐹𝐶0 𝑗 , |𝐹𝐶0 𝑗 | ∈ (1, |𝐶 𝑗 |) it holds that Í Í 𝑒 𝛾𝑒 ≤ 𝑒 𝛾¯𝑒 ≤ O (1), and in turn E[|𝐹𝐶 𝑗 |] = O (𝑡). Combining the above three cases we have 15

′ 𝑒𝑚+1

−𝑅2

′ 𝑒𝑚

𝑒 1′

...

𝑒 1𝑒 2

− 𝑚1 − 𝑟 − 𝑚1 − 𝑚1 + 𝑟 −𝑟

...

𝑒𝑚

𝑒𝑚+1

𝑟

𝑅1 R

Figure 2: Instance used to prove Theorem 2 (see Section 4.1.3). 𝐶 1 = {𝑒 1, . . . , 𝑒𝑚+1 } and 𝐶 2 = ′ {𝑒 1′ , . . . , 𝑒𝑚+1 }. that any (random) choice of set 𝐹𝐶0 𝑗 will yield E[|𝐹𝐶 𝑗 |] = O (𝑡). Applying the Chernoff bound, we obtain that |𝐹𝐶 𝑗 | = O (𝑡) with probability 1 − 𝛿/(5𝑘) by the choice of 𝑡 as from Lemma 1, proving our statement. □ Observe that the number of distance computations performed by silh-pps-all is bounded by 𝑘 ∑︁ 𝑗=1

|𝐶 𝑗 |(|𝐹𝐶0 𝑗 | + |𝐹𝐶 𝑗 |) ,

(7)

where for ease of notation |𝐹𝐶0 𝑗 | = 0 when |𝐶 𝑗 | ≤ 𝑡. By combining Equation (7), together with Lemma 3, we obtain the bound on the sample size and the bound on the number of distance computations claimed in Theorem 1, which hold with probability at least 1 − 2𝛿/5—yielding Theorem 1 by Lemmas 1 and 2. 4.1.3

Proof of Theorem 2

We now show that simple Poisson sampling over each cluster 𝐶 𝑗 , 𝑗 ∈ [𝑘], with an expected sample size equal to silh-pps-all, is not sufficient to obtain an accurate estimate of weights 𝑊𝐶 𝑗 (𝑒), 𝑗 ∈ [𝑘], 𝑒 ∈ 𝑉 . That is, we will provide an instance of Problem 1 where a uniform sampling technique with a sample complexity of Θ(𝑡𝑘), i.e., roughly the same sample size of silh-pps-all, outputs estimates of 𝑠 (𝑒), 𝑒 ∈ 𝑉 with constant error Ω(1) for a large fraction of elements in 𝑉 , proving Theorem 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 𝑉 = {𝑒 1, . . . , 𝑒𝑚+1, 𝑒 1′ , . . . , 𝑒𝑚 𝑚+1 }, |𝑉 | = 2𝑚 + 2 for fixed and ′ very large 𝑚 and C = {𝐶 1, 𝐶 2 } such that 𝐶 1 = {𝑒 1, . . . , 𝑒𝑚+1 }, 𝐶 2 = {𝑒 1′ , . . . , 𝑒𝑚+1 } and fix 𝑑 to be the Euclidean distance over R ⊃ 𝑉 .5 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 𝑒𝑖 < 𝑒 𝑗 (resp. 𝑒𝑖′ > 𝑒 ′𝑗 ) if it holds 1 ≤ 𝑖 < 𝑗 ≤ 𝑚 + 1. The elements of 𝐶 1 are placed as follows: (𝑖) the 𝑚 elements 𝑒𝑖 , 1 ≤ 𝑖 ≤ 𝑚, span an interval of length 2𝑟 , 𝑟 = 𝑜 (1/𝑚), 𝑟 > 0 centered at the origin (i.e., 𝑑 (𝑒𝑖 , 0) ≤ 𝑟 ); (𝑖𝑖) 𝑒𝑚1 +1 ≥ 0 is such that 𝑑 (𝑒𝑚1 +1, 0) ≥ 𝑅1 , for some large value 𝑅1 (that we will specify later). The structure of 𝐶 2 is similar to the structure of 𝐶 1 : the 𝑚 elements 𝑒𝑖′, 1 ≤ 𝑖 ≤ 𝑚 ′ are within an interval of length 𝑟 centered at −1/𝑚, and 𝑒𝑚+1 < 0 at distance 𝑅2 > 0, 𝑅2 = Θ(𝑚) ′ from 0. Furthermore, let the points 𝑒𝑖 (resp. 𝑒𝑖 ) for 𝑖 ∈ [𝑚] be equally spaced over the intervals ′ )) equals 𝛽𝑟 with 𝛽 2/(𝑚 − 1), for 𝑖 ∈ [𝑚 − 1]. containing them, that is 𝑑 (𝑒𝑖 , 𝑒𝑖+1 ) (resp. 𝑑 (𝑒𝑖′, 𝑒𝑖+1 We display the geometry of the resulting dataset in Figure 2. Finally let 𝑅1 ≥ 𝐿′ (𝑅2 + 1) for arbitrary large 𝐿′. 5 Recall for 𝑎, 𝑏 ∈ R the Euclidean distance corresponds to 𝑑 (𝑎, 𝑏) = |𝑎 − 𝑏 |, denoted also as length.

Proposition 1. Consider dataset 𝑉 and its clustering C as constructed above. Then, for each element 𝑒𝑖 , 𝑖 ∈ [𝑚] there exists a large constant 𝐿 such that 𝑠 (𝑒𝑖 ) = (1−𝐿) ≈ −1. 𝐿

Proof. Fix an element from 𝐶 1 in the closed interval of radius 𝑟 centered at the origin, i.e., 𝑒𝑖 , 𝑖 ∈ [𝑚]. Then consider the value of 𝑎(𝑒𝑖 ), 𝑖 ∈ [𝑚] in the definition of 𝑠 (𝑒𝑖 ), we have that 𝑚 𝑊𝐶 (𝑒𝑖 ) 𝑑 (𝑒𝑚+1, 𝑒𝑖 ) 𝑚 − 1 1 ∑︁ (𝑅1 + 𝑟 ) 𝑅1 𝑎(𝑒𝑖 ) = = 𝑑 (𝑒 𝑗 , 𝑒𝑖 ) + ≤ 2𝑟 + ≤ + 3𝑟 , 𝑚 𝑚 𝑗=1,𝑗≠𝑖 𝑚 𝑚 𝑚 𝑚 by using the triangle inequality. We also have the following lower bound on 𝑎(𝑒𝑖 ),

𝑚 𝑊𝐶 (𝑒𝑖 ) 1 ∑︁ 𝑑 (𝑒𝑚+1, 𝑒𝑖 ) 𝑚 − 1 (𝑅1 − 𝑟 ) 𝑅1 + 𝑟 𝑎(𝑒𝑖 ) = = 𝑑 (𝑒 𝑗 , 𝑒𝑖 ) + ≥ 𝛽𝑟 + ≥ , 𝑚 𝑚 𝑗=1,𝑗≠𝑖 𝑚 𝑚 𝑚 𝑚

where in the first inequality we used the fact that the distance between two elements 𝑒𝑖 , 𝑒 𝑗 with 1 ≤ 𝑖 ≠ 𝑗 ≤ 𝑚 of 𝐶 1 is 𝛽𝑟, 𝛽 = 2/(𝑚−1) and that 𝑑 (𝑒𝑚+1, 𝑒𝑖 ) ≥ 𝑅1 −𝑟 . Hence 𝑎(𝑒𝑖 ) = 𝑅1 /𝑚+𝑜 (1/𝑚) by the fact that 𝑟 = 𝑜 (1/𝑚). Next for the term 𝑏 (𝑒𝑖 ), 𝑖 ∈ [𝑚] with a similar argument we have that 𝑚 ′ ,𝑒 ) 𝑑 (𝑒𝑚+1 𝑊𝐶2 (𝑒𝑖 ) 1 ∑︁ 𝑖 ′ 𝑏 (𝑒𝑖 ) = = 𝑑 (𝑒 𝑗 , 𝑒𝑖 ) + ≤ 𝑚+1 𝑚 + 1 𝑗=1 𝑚+1 (𝑅2 + 𝑟 ) 𝑟 +1 𝑅2 𝑚 1 + 2𝑟 + ≤ + 2𝑟 + . ≤ 𝑚+1 𝑚 𝑚+1 𝑚+1 𝑚+1

(e.g., by using the triangle inequality). Finally, for the lower bound,

𝑚 ′ ,𝑒 ) 𝑑 (𝑒𝑚+1 𝑊𝐶2 (𝑒𝑖 ) 1 ∑︁ 𝑖 𝑏 (𝑒𝑖 ) = = ≥ 𝑑 (𝑒 ′𝑗 , 𝑒𝑖 ) + 𝑚+1 𝑚 + 1 𝑗=1 𝑚+1 𝑚 (𝑅2 − 𝑟 ) 1 − 2𝑚𝑟 − 𝑟 𝑅2 1 ≥ − 2𝑟 + ≥ + . 𝑚+1 𝑚 𝑚+1 𝑚+1 𝑚+1

Hence 𝑏 (𝑒𝑖 ) = (𝑅2 + 1)/(𝑚 + 1) + 𝑜 (1/𝑚). Focusing on the dominating terms of 𝑎(𝑒𝑖 ) and 𝑏 (𝑒𝑖 ) and setting 𝑅1 = 𝐿(𝑅2 + 1)𝑚/(𝑚 + 1), we have that 𝑠 (𝑒𝑖 ) =

(𝑅2 + 1)/(𝑚 + 1) − 𝑅1 /𝑚 (1 − 𝐿)(𝑅2 + 1)/(𝑚 + 1) (1 − 𝐿) 1 = = = −1 + ≈ −1 . max{(𝑅2 + 1)/(𝑚 + 1), 𝑅1 /𝑚} 𝐿(𝑅2 + 1)/(𝑚 + 1) 𝐿 𝐿

□ Comparing uniform sampling with our algorithm silh-pps-all. Consider now a uniform Poisson sampling strategy (denoted as uniform sampling for short), that estimates 𝑊𝐶 𝑗 (𝑒), 𝑗 = 1, 2 as follows: 1. Collect a sample 𝑆𝐶 𝑗 , 𝑗 = 1, 2 of points from each cluster, including each point with probability 𝑝 = 𝑡/|𝐶 𝑗 | = 𝑡/(𝑚 + 1), 𝑝 > 0 for 𝑡 = 𝑜 (𝑚) (i.e., the same 𝑡 used by silh-ppsall); 17

2. Estimate 𝑊𝐶 𝑗 (𝑒), 𝑗 = 1, 2 as: b U (𝑒) = 𝑊 𝐶𝑗

∑︁ 𝑑 (𝑒, 𝑒 ′)

𝑒 ′ ∈𝑆𝐶 𝑗

𝑝

∑︁ 𝑑 (𝑒, 𝑒 ′)(𝑚 + 1)

𝑒 ′ ∈𝑆𝐶 𝑗

𝑡

b U (𝑒) are unbiased. We can also use the estimators It is easy to show that the estimators 𝑊 𝐶𝑗

b U (𝑒) to approximate the silhouette 𝑠 (𝑒) in an analogous way as done in Equation (5). That is 𝑊 𝐶𝑗 for each 𝑒 ∈ 𝐶 1 , we consider the estimate b U (𝑒)/(𝑚 + 1) − 𝑊 b U (𝑒)/𝑚 𝑊 𝐶2 𝐶1 n o . b 𝑠 U (𝑒) = b U (𝑒)/(𝑚 + 1), 𝑊 b U (𝑒)/𝑚 max 𝑊 𝐶2 𝐶1

In the next lemma, we show that for a fixed sample size 𝑡 = 𝑜 (𝑚) and with arbitrary large probability close to 1, it holds that |b 𝑠 U (𝑒) − 𝑠 (𝑒)| ≈ 2. That is, the estimation error of the uniform sampling strategy corresponds to the maximum error attainable, for all 𝑒 ∈ 𝐶 1 \ {𝑒𝑚+1 }. In contrast, the estimates provided by silh-pps-all, using the same value of 𝑡 achieve, for 𝑒 ∈ 𝐶 1 , |b 𝑠 (𝑒) −𝑠 (𝑒)| ≈ 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 b 𝑠 U (𝑒) over all elements in 𝑒 ∈ 𝐶 2 . Proposition 2. Let 𝑉 and C = {𝐶 1, 𝐶 2 } be the instance to Problem 1 constructed above. Then, for each element in 𝑒𝑖 ∈ 𝐶 1 with 𝑖 ∈ [𝑚] the value of the uniform Poisson sampling strategy yields b 𝑠 U (𝑒𝑖 ) ≈ 1 while silh-pps-all yields b 𝑠 (𝑒𝑖 ) ≈ −1 with arbitrarily large probability. Therefore, |𝑠 (𝑒𝑖 ) − b 𝑠 U (𝑒𝑖 )| ≈ 2, while |𝑠 (𝑒𝑖 ) − b 𝑠 (𝑒𝑖 )| ≈ 0.

Proof. Let 𝑎 be a small constant. For ease of presentation and given that 𝑚 is very large, our analysis considers the following case: “the uniform Poisson sampling strategy over 𝐶 𝑗 yields 𝑎𝑡 total sampled elements”. Note that since the number of samples corresponds to a Binomial distribution Bin(𝑚 + 1, 𝑡/(𝑚 + 1)) for large 𝑚, this is often the case, up to a very small constant. A similar result holds for the size of 𝐹𝐶 𝑗 , 𝑗 ∈ [𝑘] of our algorithm silh-pps-all (see the proof of Lemma 3). Therefore, for simplicity, in our analysis we assume that both |𝑆𝐶 𝑗 | = |𝐹𝐶 𝑗 | = 𝑎𝑡, exactly. Consider the uniform sampling strategy and define the following two events 𝐸 1 ′ “𝑒𝑚+1 is not sampled in 𝑆𝐶1 ” and 𝐸 2 “𝑒𝑚+1 is not sampled in 𝑆𝐶2 ”. We have that P[𝐸 1 ] = P[𝐸 2 ] = (1 − 𝑡/(𝑚 + 1)) = 1 − 𝑜 (1) for arbitrary large 𝑚, by the choice of 𝑡. Hence, conditioning on 𝐸 1 and 𝐸 2 , we obtain that for 𝑒 ∈ 𝐶 1 \ {𝑒𝑚+1 },

and

b U (𝑒) ≥ 𝑚 + 1 𝛽𝑟𝑎𝑡 ≥ (𝑚 + 1)2𝑟𝑎 , 𝑊 𝐶1 𝑡 (𝑚 − 1)

1 𝑎(𝑚 + 1) + 2𝑟 𝑎𝑡 ≤ + (𝑚 + 1)2𝑟𝑎 , 𝑡 𝑚 𝑚 𝑚+1 1 𝑎(𝑚 + 1) U b 𝑊𝐶2 (𝑒) ≥ − 2𝑟 𝑎𝑡 ≥ − 𝑎(𝑚 + 1)2𝑟 . 𝑡 𝑚 𝑚 b U (𝑒) ≤ 𝑚 + 1 𝑊 𝐶2

b U (𝑒) ≤ 𝑚 + 1 2𝑟𝑎𝑡 = (𝑚 + 1)2𝑟𝑎 𝑊 𝐶1 𝑡

b U (𝑒) approaches 0 while 𝑊 b U (𝑒) = 𝑎(𝑚 + 1)/𝑚. Hence, for 𝑚 large and 𝑟 = 𝑜 (1/𝑚) we have that 𝑊 𝐶1 𝐶2 Therefore, focusing on the dominant terms we have that b 𝑠 U (𝑒) =

1/𝑚 + 𝑜 (1) ≈ 1 . 1/𝑚

b𝐶 𝑗 (𝑒), 𝑗 = 1, 2 computed by silh-pps-all. Where we condition Next consider the estimates 𝑊 b𝐶1 (𝑒) ¯ 𝑒)/𝑊 ¯ for on the events 𝐸 1 and 𝐸 2 defined over 𝐹𝐶0 𝑗 . Then we have that 𝛾𝑒 ≥ max𝑒∈𝐹 ¯ 𝐶0 𝑑 (𝑒, 1 𝑒 = 𝑒𝑚+1 that is 𝛾𝑒 ≥ (𝑅1 − 𝑟 )/(𝑅1 + 𝑟 + 2𝑚𝑟 ) = 1 − 𝑜 (1). Therefore over 𝑡 points selected in the final sample 𝐹𝐶1 , point 𝑒𝑚+1 will be selected with probability arbitrary close to 1 by silh-pps-all, ′ by the choice of 𝑟 . A similar argument applies to 𝑒𝑚+1 . We therefore have that for 𝑒 ∈ 𝐶 1 \ {𝑒𝑚+1 } it holds ∑︁ 𝑑 (𝑒, 𝑒 ′) b𝐶1 (𝑒) = 𝑑 (𝑒, 𝑒𝑚+1 ) + 𝑊 = 𝑅1 + 𝑜 (𝑅1 ) ′ 𝑝 𝑒 ′ 𝑒 ∈𝐹 ∑︁ 𝑑 (𝑒, 𝑒 ′) 𝐶1

′ b𝐶2 (𝑒) = 𝑑 (𝑒, 𝑒𝑚+1 𝑊 )+

𝑒 ′ ∈𝐹𝐶2

𝑝𝑒′

= 𝑅2 + 𝑜 (𝑅2 )

by the facts that 𝑟 = 𝑜 (1/𝑚), 𝑡 = 𝑜 (𝑚), 𝑅2 = Θ(𝑚) and that 𝐶 2 is centered at −1/𝑚. Hence, by the choice of 𝑅1 = 𝐿(𝑅2 + 1), for a suitably large 𝐿, focusing on the dominant terms b 𝑠 (𝑒) =

𝑅2 /(𝑚 + 1) − 𝑅1 /𝑚 ≈ −1 . max{𝑅2 /(𝑚 + 1), 𝑅1 /𝑚}

The final claim follows by computing the distance of b 𝑠 U (𝑒) and b 𝑠 (𝑒) using the exact silhouette values of 𝑒𝑖 , 𝑖 ∈ [𝑚] from Proposition 1, i.e., 𝑠 (𝑒) ≈ −1. □

Proving our result. The final claim of Theorem 2 follows as a corollary of Proposition 2, and noting that |𝑉 | = 𝑛 = 2𝑚 + 2 and |𝐶 1 | = 𝑚 + 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 𝑠 (𝑒), 𝑒 ∈ 𝑉 that cannot be matched by a uniform sampling strategy, used to approximate the values of 𝑊𝐶 𝑗 (𝑒), 𝑗 ∈ [𝑘], when both approaches retain the same sample complexity. We further complement our theoretical insights with our extensive experimental evaluation in Section 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 C, as captured in our next statement.

Problem 2. Given a set of elements 𝑉 and a 𝑘-clustering C of 𝑉 and two parameters 𝜀, 𝛿 ∈ (0, 1), obtain an estimator 𝑠 ′ (C) such that |𝑠 ′ (C) −𝑠 (C)| ≤ 𝑓 (𝜀) with probability at least 1−𝛿, minimizing the number of distance computations required to compute 𝑠 ′ (C). Recall that 𝑠 (C) is largely used in practice for evaluating the quality of a clustering C—where C is often obtained optimizing popular clustering metrics such as the 𝑘-means or 𝑘-medoids 19

objective functions [Dudek, 2020]. Hence, obtaining an accurate approximation of 𝑠 (C) is of high practical utility, e.g., for identifying a good value for the parameter 𝑘 on large datasets where 𝑠 (C) cannot be computed exactly. In the next sections we will obtain several estimators for the value 𝑠 (C), solving Problem 2, where, in analogy with Problem 1, we want 𝑓 (𝜀) = Θ(𝜀). In addition, we want to output our highly accurate estimate 𝑠 ′ (C) performing significantly less distance computations required by the exact algorithm, i.e., Θ(𝑛 2 ). 4.2.1

Estimator based on the silhouette definition

Our first estimator (global-sampler, gl-s) is based on the definition of 𝑠 (C) itself (see Equation 3). That is, 𝑠 (C) corresponds to the average of the silhouette of all the elements 𝑒 ∈ 𝑉 , where each term 𝑠 (𝑒) 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 𝑉 , and using the average of their exact silhouette value (i.e., 𝑠 (𝑒)) as proxy for the overall silhouette 𝑠 (C). The estimator is computed as follows: we sample {𝑒 1, . . . , 𝑒𝑚 } random elements from 𝑉 where each element 𝑒 ℓ , ℓ ∈ [𝑚] is sampled uniformly at random from the set 𝑉 . The resulting estimator can be expressed as 𝑚 1 ∑︁ 𝑠 (𝑒 ℓ ) . (8) b 𝑠 1 (C) = 𝑚 ℓ=1 The next result quantifies the guarantees offered by estimator b 𝑠 1 (C).

Theorem 3. Let 𝑉 be a dataset of 𝑛 elements, and let C be a 𝑘-clustering of 𝑉 . If 𝑚 ≥ 𝜀22 ln 𝛿2 , then |b 𝑠 1 (C) − 𝑠 (C)| ≤ 𝜀 with probability at least 1 − 𝛿.

Proof. Note that 𝑠 (𝑒 ℓ ), ℓ = 1, . . . , 𝑚 are random variables taking values in the interval [−1, 1]. Í Moreover, for each ℓ = 1, . . . , 𝑚, we have that E [𝑠 (𝑒 ℓ )] = 𝑒∈𝑉 𝑠 (𝑒)/𝑛 = 𝑠 (C). The result follows directly from the application of Hoeffding’s bound (Theorem 6). □ The time complexity (i.e., number of distance computations) required to compute b 𝑠 1 (C) is O (𝑚𝑛), for fixed 𝑚.

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 𝑉 , a clustering C of 𝑉 and two parameters 𝜀, 𝛿 ∈ (0, 1) computes a set b 𝑠 (𝑉 ) such that for each 𝑒 ∈ 𝑉 it holds |b 𝑠 (𝑒) − 𝑠 (𝑒)| ≤ (4𝜀)/(1 − 𝜀) with probability at least 1 − 𝛿. Our second estimator (global-pps-full-sampler, gl-pps-f) approximates the silhouette 𝑠 (C) using the estimates b 𝑠 (𝑒) for each element 𝑒 ∈ 𝑉 as obtained by silh-pps-all, that is b 𝑠 2 (C) =

1 ∑︁ b 𝑠 (𝑒) , 𝑛 𝑒∈𝑉

(9)

where the values b 𝑠 (𝑒), 𝑒 ∈ 𝑉 are computed by silh-pps-all. The next theorem establishes an upper bound to the absolute error of b 𝑠 2 (C) with respect to 𝑠 (C). 20

Theorem 4. Let 𝑉 be a dataset of 𝑛 elements, and let C be a 𝑘-clustering of 𝑉 . Let b 𝑠 2 (C) be the estimate of the silhouette of the clustering 𝑠 (C) computed as in Equation (9) by running silh-pps-all for given parameters 𝜀 and 𝛿, with 0 < 𝜀, 𝛿 < 1, and for a suitable choice of the sample size 𝑡. Then, with probability at least 1 − 𝛿, it holds that |b 𝑠 2 (C) − 𝑠 (C)| ≤

4𝜀 . 1−𝜀

The proof is an immediate consequence of the definition of 𝑠 (C) (Equation (3)) and b 𝑠 2 (C) (Equation (9)), and of the guarantees provided by Lemmas 1 and 2. For a given (expected) sample size 𝑡,6 the time complexity to compute b 𝑠 2 (C) is the same required by silh-pps-all, that is O (𝑛𝑘𝑡) with probability at least 1 − 𝛿 (see Lemma 3, and the related discussion). 4.2.3

Estimator using two-phase sampling

Let us consider, for ease of notation, the set b 𝑠 (𝑉 ) = {b 𝑠 (𝑒 1 ), . . . ,b 𝑠 (𝑒𝑛 )} of estimates computed with silh-pps-all. Note that the values in the set b 𝑠 (𝑉 ) are fixed once the pps sampling phase of Algorithm 1 (i.e, Phase 1) has been performed, even if the values in b 𝑠 (𝑉 ) have not explicitly 7 been computed yet. We therefore approximate the estimator in Equation (9) as follows: 1. Execute Phase 1 of Algorithm 1 and obtain all samples 𝐹𝐶 𝑗 , 𝑗 ∈ [𝑘], so that the values of b 𝑠 (𝑒𝑖 ), 𝑖 ∈ [𝑛] are determined (i.e., the set b 𝑠 (𝑉 )), even if not computed yet; 2. Select uniformly at random 𝑚 ∈ N, 𝑚 ≥ 2 indices 𝑖 1, . . . , 𝑖𝑚 in the set {1, . . . , 𝑛}; 3. For each 𝑖 𝑗 , 𝑗 ∈ {1, . . . , 𝑚} evaluate b 𝑠 (𝑒𝑖 𝑗 );

4. Return the average 𝑠¯(C) over the 𝑚 estimates b 𝑠 (𝑒𝑖 𝑗 ) obtained above.

The resulting estimator (global-pps-sub-sampler, gl-pps-s) can be expressed as follows, 1 ∑︁ 𝑠¯(C) = b 𝑠 (𝑒𝑖 𝑗 ) . 𝑚 𝑗=1 𝑚

(10)

The following result shows that 𝑠¯(C) is an unbiased estimator of b 𝑠 2 (C).

Lemma 4. E [¯𝑠 (C)] = b 𝑠 2 (C).

Proof. For each 𝑗 ∈ {1, . . . , 𝑚}, let 𝑋 𝑗 be a random variable with value 𝑠ˆ(𝑒𝑖 𝑗 ) (observe that 𝑖 𝑗 is Í itself a random variable). Note that 𝑠¯(C) = 𝑚1 𝑚𝑗=1 𝑋 𝑗 . By taking expectations, using linearity, and the fact that r.v.’s 𝑖 𝑗 ’s are taken uniformly at random from the set {1, . . . , 𝑛} we have 𝑚 𝑚 𝑛 𝑛 1 ∑︁ 1 ∑︁ ∑︁ 1 1 ∑︁ E[¯𝑠 (C)] = E 𝑋𝑗 = b 𝑠 (𝑒𝑖 ) = b 𝑠 (𝑒𝑖 ) = b 𝑠 2 (C). 𝑚 𝑗=1 𝑚 𝑗=1 𝑖=1 𝑛 𝑛 𝑖=1

6 Recall that the number of samples in each cluster in silh-pps-all is O (𝑡) with probability at least 1 − 2𝛿/5.

□

7 After Phase 1, the values b 𝑠 (𝑒 1 ), . . . ,b 𝑠 (𝑒𝑛 ) are determined, i.e., there is no randomness once the samples 𝐹𝐶 𝑗 for

𝑗 ∈ [𝑘] are fixed.

Table 1: Summary of the various estimators to approximate 𝑠 (C), the silhouette of a clustering C, as captured by Problem 2. We label with algorithm the method used to evaluate each estimator. The sample size 𝑚(𝜀 2, 𝛿 2 ) is detailed in Lemma 5. Algorithm Estimator gl-s gl-pps-f gl-pps-s

b 𝑠 1 (C)

b 𝑠 2 (C) 𝑠¯(C)

Guarantees |b 𝑠 1 (C) − 𝑠 (C)| ≤ 𝜀

4𝜀 |b 𝑠 2 (C) − 𝑠 (C)| ≤ 1−𝜀

4𝜀 1 |¯𝑠 (C) − 𝑠 (C)| ≤ 1−𝜀 + 𝜀2 1

Distance-Computations O 𝜀𝑛2 log 𝛿1 𝑛𝑘 O 𝑛𝑘 log 2 𝜀 𝛿 O 𝑛 + 𝑚(𝜀𝜀2,𝛿2 2 )𝑘 log 𝑛𝑘 𝛿1 1

Probability > 1−𝛿

> 1−𝛿

> (1 − 𝛿 1 ) (1 − 𝛿 2 )

The following result provides the guarantees on the approximation of b 𝑠 2 (C) provided by its estimator 𝑠¯(C). 2 2 Lemma 5. Fix (𝜀 2, 𝛿 2 ) ∈ (0, 1) then if 𝑚 ≥ 𝜀 2 ln 𝛿22 then P[|b 𝑠 2 (C) − 𝑠¯(S)| ≥ 𝜀 2 ] ≤ 𝛿 2 . 2

Proof. First recall variables 𝑋 𝑗 , 𝑗 ∈ [1, 𝑚] are bounded in the range [−1, 1]. Then combining the proof of Lemma 4 showing that for each 𝑗 ∈ [𝑚], E[𝑋 𝑗 ] = b 𝑠 (C) and Hoeffding’s bound (Theorem 6) we have that 2𝑚𝜀 22 P [|¯𝑠 (C) − b 𝑠 2 (C)| ≥ 𝜀 2 ] ≤ 2 exp − ≤ 𝛿2 , 4

where the last step is by the choice of 𝑚 as in statement.

□

We now prove the resulting error obtained from approximating the silhouette 𝑠 (C) with 𝑠¯(C).

and such event holds with probability at least (1 − 𝛿 1 )(1 − 𝛿 2 ).

□

It now remains to control |𝑠 (C) − 𝑠¯(C)| ≤ 𝜀 1 + 𝜀 2 to hold with probability at least 1 − 𝛿 within user defined accuracy and confidence 𝜀, 𝛿 ∈ (0, 1). Since the final error is bounded by 𝜀 1 + 𝜀 2 we can set 𝜀 2 = 𝜀 − 𝜀 1 . Noting that 𝜀 1 ≤ 4𝜀 ′/(1 − 𝜀 ′) it is sufficient to set 𝜀 ′ = 𝜀/𝑞, 𝑞 > 4. Finally to control the confidence we can set 𝛿 1 = 𝛿/2 and 𝛿 2 = 𝛿/(2 − 𝛿), yielding to the desired confidence to be at least 1 − 𝛿—solving Problem 2 as desired. When performing Phase 1 of silh-pps-all for a fixed 𝑡, the time complexity to compute 𝑠¯(C) is at most O (𝑛 log(𝑛𝑘/𝛿) + 𝑚𝑘𝑡). That is the sum of the time complexity of Phase 1 (i.e., O (𝑛 log(𝑛𝑘/𝛿)) and the time complexity of running Phase 2 of silh-pps-all only for the 𝑚 elements 𝑒𝑖 𝑗 , 𝑗 ∈ {1, . . . , 𝑚} (i.e., O (𝑚𝑘𝑡)). We provide a summary of our sampling techniques designed to solve Problem 2 in Table 1. 22

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 [Dean and Ghemawat, 2008, Leskovec et al., 2014, Pietracaprina et al., 2012]. At the end of the section, we will briefly discuss how to simply adapt our design to the popular MPC model [Beame et al., 2017, Im et al., 2023, Karloff et al., 2010]. 4.3.1

Computational model

MapReduce adopts a functional programming model where an algorithm is specified as a sequence of rounds. In each round: 1. a multiset 𝑋 of key-value pairs is first transformed into a new multiset 𝑋 ′ of key-value pairs by applying a map function independently to each pair (possibly in parallel); 2. the multiset 𝑋 ′ is then transformed into a multiset 𝑌 of pairs by applying a reduce function independently to each subset of key-value pairs of 𝑋 ′ sharing the same key (possibly in parallel); 3. the multiset of key-value pairs 𝑌 is then used as an input to the subsequent round, if any. Throughout the description of our algorithms, a key-value pair will be denoted as (key | value1, . . . , value𝑟 ), where ⟨value1, . . . , value𝑟 ⟩ corresponds to an 𝑟 -dimensional value, for (a small) constant 𝑟 . As customary in the MapReduce literature Leskovec et al. [2014], 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 (2) 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: 𝑖) the number of rounds; 𝑖𝑖) the maximum amount 𝑀𝐿 of local memory required by any instance of the map or reduce functions; 𝑖𝑖𝑖) the maximum amount 𝑀𝐴 of aggregate memory required at any round to store the input, output and intermediate data. In our setting, consider a 𝑘-clustering C = {𝐶 1, . . . , 𝐶𝑘 } of a dataset 𝑉 = {𝑒 0, . . . 𝑒𝑛−1 } of 𝑛 elements, and assume that each element 𝑒 ∈ 𝑉 has additional information 𝐶 (𝑒) denoting the cluster to which 𝑒 ∈ 𝑉 is assigned to over C. Initially, we represent the input clustering by the following set of key-value pairs: {(𝑖 |𝑒𝑖 ) : 0 ≤ 𝑖 < 𝑛} .

Reasonably, we assume that the values 𝑛, 𝑘, and 𝑡 used in the MapReduce algorithm, and the size |𝐶 𝑗 | for each 𝑗 ∈ [𝑘], are known globally. A key ingredient of our design is that in each round, the elements of 𝑉 , and the newly created intermediate data, will be partitioned into 𝑤 subsets of size 𝑛/𝑤 each, where 𝑤 ∈ [0, 𝑛 − 1] is a suitable design parameter. To simplify the description and the analysis of our methods, we assume that 𝑤 ∈ O (𝑛𝛼 ), for some (small) 𝛼 < 1/2. While our methods can be generalized to larger values of 𝑤 (i.e., those corresponding to 𝛼 ≥ 1/2) with minimal modifications, we highlight that 𝑤 constitutes an upper√bound to the maximum parallelism exploitable by the algorithm. Thus, assuming that 𝑤 ∈ O 𝑛 does not constitute a limitation in practice, especially for modern datasets of massive size 𝑛. 23

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 𝐹𝐶 𝑗 , for every cluster 𝐶 𝑗 and 𝑗 ∈ [𝑘], while the last round computes the estimates b 𝑠 (𝑒) of 𝑠 (𝑒) of each element 𝑒 ∈ 𝑉 . Round 1. • Map. Map each pair (𝑖 |𝑒𝑖 ) into the pair (𝑖 mod 𝑤 |𝑒𝑖 , 0). In addition, if |𝐶 (𝑒𝑖 )| > 𝑡 then, with probability (2/|𝐶 (𝑒𝑖 )|) ln(2𝑘/𝛿), select 𝑒𝑖 to be part of the initial Poisson sample 𝐹𝐶(0)(𝑒 ) , 𝑖 producing 𝑤 additional pairs (ℓ |𝑒𝑖 , 1), with 0 ≤ ℓ < 𝑤. • Reduce. For 0 ≤ ℓ < 𝑤, the reducer associated with key ℓ has a copy of all initial samples 𝐹𝐶(0) (represented by the pairs (ℓ |𝑒𝑖 , 1)) and a subset 𝑉ℓ of 𝑉 (represented by the pairs 𝑗 (ℓ |𝑒𝑖 , 0)). Observe that the sets 𝑉ℓ , ℓ ∈ [0, 𝑤) form a balanced partition of 𝑉 into 𝑤 subsets. Then, for each element 𝑒𝑖 ∈ 𝐹𝐶(0)(𝑒 ) , the reducer computes the term 𝑊𝑖,ℓ , defined as the 𝑖 sum of distances from 𝑒𝑖 to all elements of 𝑉ℓ ∩ 𝐶 (𝑒𝑖 )) producing the pair (ℓ |𝑒𝑖 ,𝑊𝑖,ℓ , 1). In addition, for each pair (ℓ |𝑒𝑖 , 0) produce the pair (ℓ |𝑒𝑖 , 0, 0). Round 2. • Map. Map each pair (ℓ |𝑒𝑖 ,𝑊𝑖,ℓ , 1) into the 𝑤 pairs (0|𝑒𝑖 ,𝑊𝑖,ℓ , 1), (1|𝑒𝑖 ,𝑊𝑖,ℓ ), . . . (𝑤 − 1|𝑒𝑖 ,𝑊𝑖,ℓ , 1), and each pair (ℓ |𝑒𝑖 , 0, 0) into itself. • Reduce. For 0 ≤ ℓ < 𝑤, the reducer associated with key ℓ has the following pairs: for every element 𝑒𝑖 belonging to the initial sample 𝐹𝐶(0)(𝑒 ) of its cluster, 𝑤 pairs (ℓ |𝑒𝑖 ,𝑊𝑖,ℓ ′ , 1), 𝑖 with 0 ≤ ℓ ′ < 𝑤; and for every 𝑒𝑖 ∈ 𝑉ℓ , a pair (ℓ |𝑒𝑖 , 0, 0). The reducer first computes 𝑊𝐶 (𝑒𝑖 ) (𝑒𝑖 ), for each element 𝑒𝑖 ∈ 𝐹𝐶(0)(𝑒 ) , by summing all 𝑊𝑖,ℓ ′ ’s. Next, for each 𝑒𝑖 ∈ 𝑉ℓ it 𝑖 produces the pair (ℓ |𝑒𝑖 , 𝑝 (𝑒𝑖 )), where: 1) if |𝐶 (𝑒𝑖 )| ≤ 𝑡, then 𝑝 (𝑒𝑖 ) = 1; otherwise 2) 𝑝 (𝑒𝑖 ) is computed in terms of the 𝑊𝐶 (𝑒𝑖 ) (𝑒) with 𝑒 ∈ 𝐹𝐶0 (𝑒𝑖 ) , as specified in Line 1 of Algorithm 1. Round 3. • Map. Map each pair (ℓ |𝑒𝑖 , 𝑝 (𝑒𝑖 )) into the pair (ℓ |𝑒𝑖 , 𝑝 (𝑒𝑖 ), 0). In addition, with probability 𝑝 (𝑒𝑖 ) select 𝑒𝑖 to be part of the Poisson sample 𝐹𝐶 𝑗 producing pairs (ℓ ′ |𝑒𝑖 , 𝑝 (𝑒𝑖 ), 1), with 0 ≤ ℓ ′ < 𝑤. • Reduce. For 0 ≤ ℓ < 𝑤, the reducer associated with key ℓ has a copy of all samples 𝐹𝐶 𝑗 (represented by the pairs (ℓ |𝑒𝑖 , 𝑝 (𝑒𝑖 ), 1)) and the subset 𝑉ℓ (represented by the pairs (ℓ |𝑒𝑖 , 𝑝 (𝑒𝑖 ), 0)). Then, for each 𝑒𝑖 ∈ 𝑉ℓ , the reducer computes the estimate b 𝑠 (𝑒𝑖 ), and generates the pair (0|𝑒𝑖 ,b 𝑠 (𝑒𝑖 )) to be reported in output. Recall that Lemma 3 states that with probability at least 𝐶 𝑗 with 1 − 2𝛿/5, for every cluster 0 −2 −2 |𝐶 𝑗 | > 𝑡, |𝐹𝐶 𝑗 | = O (ln(𝑘/𝛿)) and |𝐹𝐶 𝑗 | = O 𝜀 ln(𝑛𝑘/𝛿) . Then, since 𝑡 = O 𝜀 ln(𝑛𝑘/𝛿) , it immediately follows that, with probability at least 1 − 2𝛿/5 the maximum local memory required by any map or reduce phase in any of the above rounds is bounded by 𝑀𝐿 = O 𝑛/𝑤 + 𝑤 + 𝑘𝜀 −2 ln(𝑛𝑘/𝛿) = O 𝑛/𝑤 + 𝑘𝜀 −2 ln(𝑛𝑘/𝛿) and the aggregate memory is bounded by 𝑀𝐴 = O (𝑤𝑀𝐿 ). 4.3.3

Approximating the silhouette of a clustering

In Section 4.2, we introduced three methods for computing an approximation b 𝑠 (C) to the silhouette of the entire clustering C, see Table 1. The two estimators b 𝑠 2 (C) (see Equation (9)) and 𝑠¯(C) (see Equation (10)), can be computed through simple aggregations of all, or an 𝑚-sample, 24

of the individual silhouettes estimates b 𝑠 (𝑒) for 𝑒 ∈ 𝑉 . 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 𝑤 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 Section 4.3.2. Finally, concerning the sampling method that computes b 𝑠 1 (C) based on Equation (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 O (1 + (log 𝑤/log(𝑛/𝑤))) = O (1) rounds, local memory −2 𝑀𝐿 = O 𝑛/𝑤 + 𝑘𝜀 ln(1/𝛿) and aggregate memory 𝑀𝐴 = O (𝑤𝑀𝐿 ), 4.3.4

Porting MapReduce algorithms to the MPC framework

The MPC framework models a parallel system with 𝑝 machines, each equipped with a local memory of size 𝑠. Similar to MapReduce, the computation is structured as a sequence of synchronous rounds. In a round, each machine (1) performs an arbitrary amount of computation on its local data in memory; (2) exchanges point-to-point messages with other machines. The MPC framework allows arbitrary communication patterns, subject to the only constraint that each machine may send or receive at most 𝑠 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: 𝑖) a MapReduce computation is purely functional and exhibits, at each round, a degree of (virtual) parallelism equal to the number of instances of the map and reduce functions active in the round; in contrast 𝑖𝑖) an MPC computation has fixed parallelism 𝑝 throughout all rounds. 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 𝑤 disjoint subsets, and the computation essentially involves: 𝑖) local operations within each partition; 𝑖𝑖) broadcast of individual elements to all partitions; 𝑖𝑖𝑖) summing 𝑤 contributions, one per partition. As a consequence, our algorithms can be immediately adapted to the MPC framework by associating 1/2 one machine per partition, setting 𝑝 = 𝑤 and 𝑠 = 𝑀𝐿 . Finally, since in our setting 𝑀𝐿 = Ω 𝑛 , all the required broadcast and aggregation operations can be performed in constant rounds using standard MPC techniques [Czumaj et al., 2025]. 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. Theorem 5. There exist MapReduce and MPC and algorithms providing distributed implementations of silh-pps-all, for local silhouette estimation, and of all our algorithms for global silhouette estimation from Table 1. For fixed 𝜀, 𝛿 ∈ (0, 1), 𝑤 = Θ(𝑛𝛼 ), for 𝛼 ≤ 1/2, and 𝑘 = O 𝑛 1−𝛼 /log2 (𝑛) , our distributed algorithms require a constant number of rounds, sublinear local memory 𝑀𝐿 = 1−𝛼 Θ𝑛 , and linear aggregate memory, with probability at least 1 − 2𝛿/5. 25

Experimental evaluation

In this section, we present the results of our extensive experimental evaluation. We start by describing the setup of our experiments (Section 5.1), and then we investigate the following issues. • I1. Compare our new estimators (from Section 4.2) for the global silhouette 𝑠 (C) with stateof-the-art baselines, studying the trade-offs between accuracy and efficiency exhibited by our estimators (Section 5.2). • I2. Evaluate the accuracy of our algorithm silh-pps-all for the estimation of all local silhouette values 𝑠 (𝑒), for each element 𝑒 ∈ 𝑉 (Section 5.3). • I3. Study the scalability of a distributed implementation of our algorithm silh-pps-all (Section 5.4). • I4. Illustrate applications of our new algorithms: 1. our silh-pps-all algorithm provides highly accurate local estimates that can be used to visualize accurate silhouette plots; 2. our global estimates can be used for the accurate selection of a good value of the parameter 𝑘 of a 𝑘-clustering method (Section 5.5).

5.1

Setup

5.1.1

Methods and baselines

We now introduce all methods considered in our experiments. We distinguish between methods for global (Problem 2) and local estimation (Problem 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 Section 4.1.3). The list of methods is as follows. • Algorithms gl-s, gl-pps-f, and gl-pps-s (from Table 1), which compute the estimators in Equations (8) to (10), respectively, using pps sampling to estimate individual silhouette values. • Two algorithms, dubbed gl-uni-f and gl-uni-s, which compute the estimators in Equations (9) and (10), respectively, using uniform sampling (rather than pps sampling) to estimate all individual silhouettes (see Section 4.1.3). • “FS” by Frahling and Sohler [2008], which, for each element 𝑒 ∈ 𝑉 computes the term 𝑎(𝑒) according to its definition (Equation (2)), while estimates the term 𝑏 (𝑒) as Í 1/|𝐶 𝜇𝑒 | 𝑒 ′ ∈𝐶 𝜇𝑒 𝑑 (𝑒 ′, 𝑒), where 𝐶 𝜇𝑒 corresponds to the cluster with closest center (or centroid) to 𝑒 ∈ 𝑉 , different from the cluster 𝐶 of 𝑒. Observe that the estimate of 𝑏 (𝑒) 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 𝑒 and all centers (or centroids) of the clusters in C \ {𝐶 𝜇𝑒 }—to avoid performing an exact computation of 𝑏 (𝑒) when the check fails. • “SIMPL” by Hruschka et al., the simplified silhouette approach, which evaluates the silhouette 𝑠 (𝑒) of each element 𝑒 ∈ 𝑉 using the distances of 𝑒 to the closest centers (or 26

Table 2: Datasets used in our experiments. For each dataset, the first three columns report its name and label, its size 𝑛, and the dimensionality 𝑧 of its elements (i.e., 𝑒 ∈ R𝑧 ). Column “exact” indicates whether the Θ(𝑛 2 ) exact algorithm could be run within the allotted time limit; Column “𝑘-med” indicates whether the dataset could be clustered with 𝑘-medoids, under the distances listed in Section 5.1.2); and Column “size” labels the dataset as small (S), medium (M) or large (L), used in our discussion. Name (label)

𝑛

𝑧

exact 𝑘-med size

Breast (BR) 568 Wine (WI) 6.5 𝐾 Credit (CR) 30 𝐾 Shuttle (SHU) 58 𝐾 IoT (IOT) 123 𝐾

30 11 23 9 78

✓ ✓ ✓ ✓ ✓

S S M M M

Name (label) BioKDD (BioKDD) RNA-seq (RNA) Metro (MT) PowerHouse (PH) Gowalla (GOW)

𝑛

𝑧

exact 𝑘-med size

146 𝐾 489 𝐾 1.5 𝑀 2.1 𝑀 6.4 𝑀

74 8 7 7 2

✓ ✗ ✗ ✗ ✗

✗ ✗ ✗ ✗ ✗

L L L L L

centroids) as proxies of the average distances. Formally, let 𝜇 1, . . . , 𝜇𝑘 the centers of clusters 𝐶 1, . . . , 𝐶𝑘 ∈ C then for an element 𝑒 in some cluster 𝐶 the method SIMPL evaluates 𝑎(𝑒) = 𝑑 (𝑒, 𝜇) where 𝜇 is the center of 𝐶, and 𝑏 (𝑒) = min𝐶 𝑗 ≠𝐶 𝑑 (𝑒, 𝜇 𝑗 ).

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 Table 1, where the estimation error can be controlled as a function of the input parameters. Local estimates. When solving Problem 1 we will compare our method silh-pps-all with the uniform-sampling based approach detailed in Section 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 Problem 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 Table 2. We provide more details of the processing of each dataset in Section C.1. We use small-size datasets (those marked with S in Table 2) and CR to showcase applications of Problems 1 and 2 (see Section 5.5). As for medium-size datasets, they are mainly used in Section 5.2 and Section 5.3, and are processed as follows. We cluster each dataset using both the 𝑘-means and the 𝑘medoids objectives. For the 𝑘-means objective we consider the Euclidean distance 𝑑 eucl (𝑥, 𝑦) = √︁ Í 2 𝑖 (𝑥𝑖 − 𝑦𝑖 ) . For the 𝑘-medoids objective, in addition to 𝑑 eucl , we also consider several widely used distance functions: √︃Í √︃Í Í 2 2 • 𝑑 cos , the cosine distance 𝑑 cos (𝑥, 𝑦) = 1 − 𝑖 (𝑥𝑖 − 𝑦𝑖 )/( 𝑖 𝑥𝑖 𝑖 𝑦𝑖 ); Í • 𝑑 manh , the Manhattan distance 𝑑 manh (𝑥, 𝑦) = 𝑖 |𝑥𝑖 − 𝑦𝑖 |; Í • 𝑑 can , the Canberra distance 𝑑 can (𝑥, 𝑦) = 𝑖 |𝑥𝑖 − 𝑦𝑖 |/(|𝑥𝑖 | + |𝑦𝑖 |) .

For each clustering objective (either 𝑘-means or 𝑘-medoids) and each distance, we use different values of the parameter 𝑘 ∈ {2, 5, 10, 15, 20} to obtain different clusterings C. Hence, 27

each distinct configuration (dataset, value of 𝑘, distance, objective) leads to a different clustering C, which we use as an input to either Problem 2 or Problem 1. Finally, for the large datasets from Table 2 we proceeded similarly to the medium-sized datasets, but we only clustered such data using the 𝑘-means objective to obtain different clusterings C varying 𝑘, with 𝑘 ∈ {2, 5, 10, 15, 20}. We do not use 𝑘-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 Section C.2. The code to reproduce our results is available online.8 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 𝑡 controlling the expected sample size, rather than obtaining 𝑡 as a function of the other parameters (see Algorithm 1). When comparing the algorithms from Table 1 for Problem 2, and their uni variant, we will set their parameters to have a comparable runtime: we fix 𝑡 for gl-pps-f and gl-uni-f and we set the sample size 𝑚 of gl-s to 𝑘𝑡. Unless otherwise stated we set 𝛿 = 0.01. 5.1.4

Exact silhouette scores

For the small and medium sized datasets from Table 2 we were able to compute the exact silhouette scores on all the various configurations described in Section 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 𝑠 (𝑒) for each element 𝑒 ∈ 𝑉 , we then proceed as follows. We perform five independent executions of silh-pps-all and its uniform variant (uni) from Section 4.1.3, obtaining respectively b 𝑠 (𝑒)𝑖 and b 𝑠 (𝑒) with 𝑖 ∈ {1, . . . , 5} with a high sample size of 𝑡 = 800. We Í U 𝑖 then estimate 𝑠 (𝑒) = (1/10) 𝑖 (b 𝑠 (𝑒)𝑖 + b 𝑠 U (𝑒)𝑖 ) 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 Section 5.1.1. We compare all methods considering the average error |b 𝑠 (C) − 𝑠 (C)| ∈ [0, 2] over ten independent runs (see Section 5.1.3 for details), where b 𝑠 (C) corresponds to the approximation provided by the method. In our discussion we also mention the relative approximation error that corresponds to |b 𝑠 (C) − 𝑠 (C)|/𝑠 (C). For medium-sized datasets we set the parameters of the various methods (see Table 3) as: 𝑡 = 16 for both gl-pps-f and gl-uni-f, 𝑚 = ⌈𝑛/20⌉, 𝑡 = 16 for gl-pps-s and gl-uni-s, and the sample size of gl-s to 𝑚′ = 𝑡𝑘 depending on C (see Section 5.1.3). Our choice of a small value of 𝑡 is motivated by the size of our datasets, while the choice of 𝑚 enforces the computation of the estimator in Equation (10) on a strictly sublinear number of elements in 𝑛. 8 https://github.com/iliesarpe/ScalableSilhouetteComputation.

0.1

0.0

0.00

Speedup over exact

SHU, 𝑑 eucl, 𝑘 = 15, 𝑠 (C) = 0.24, 𝑟𝑡 EX = 25.4

Speedup over exact

SHU, 𝑑 cos, 𝑘 = 2, 𝑠 (C) = 0.76, 𝑟𝑡 EX = 156.1

101

104 103 102 101

0.05

0.0

0.1

101 0.0

0.1

Average error

Speedup over exact

IOT, 𝑑 manh, 𝑘 = 10, 𝑠 (C) = 0.85, 𝑟𝑡 EX = 1522.4

Speedup over exact

IOT, 𝑑 manh, 𝑘 = 5, 𝑠 (C) = 0.74, 𝑟𝑡 EX = 1516.8

105

103

101 0.000

0.025

0.050

102 101 100

0.10

0.00

0.05

0.10

IOT, 𝑘 = 5, 𝑠 (C) = 1.0, 𝑟𝑡 EX = 1743.3

IOT, 𝑑 manh, 𝑘 = 2, 𝑠 (C) = 0.57, 𝑟𝑡 EX = 1516.0

105 103 101 0.0005

105 103 101

0.0

0.0010

0.2

Average error

IOT, 𝑑 manh, 𝑘 = 15, 𝑠 (C) = 0.84, 𝑟𝑡 EX = 1512.6

IOT, 𝑑 cos, 𝑘 = 10, 𝑠 (C) = 0.96, 𝑟𝑡 EX = 2498.4

105

103

101 0.00

Average error

103

Average error

0.0000

Average error

103

0.05

104

Average error

0.2

Average error 105

100 0.00

Average error

103

101

0.1

Average error 105

102

Speedup over exact

101

103

SHU, 𝑘 = 20, 𝑠 (C) = 0.39, 𝑟𝑡 EX = 30.3

Speedup over exact

100

102

104

gl-uni-s

0.05

Average error

Speedup over exact

101

103

gl-pps-s

SHU, 𝑘 = 2, 𝑠 (C) = 0.6, 𝑟𝑡 EX = 30.1

Speedup over exact

102

104

gl-uni-f

Speedup over exact

103

gl-pps-f

CR, 𝑑 can, 𝑘 = 10, 𝑠 (C) = 0.14, 𝑟𝑡 EX = 35.8

Speedup over exact

104

0.0

gl-s

SIMPL Speedup over exact

CR, 𝑑 can, 𝑘 = 2, 𝑠 (C) = 0.33, 𝑟𝑡 EX = 35.8

105

103

101 0.00

0.01

Average error

Figure 3: Trade-off between accuracy and efficiency over medium sized datasets (see Table 2). The 𝑥 axis is associated with the average error over 10 independent runs, while the 𝑦 axis is associated with the average speed-up over the exact computation. For each plot that corresponding to a different configuration, we report: the dataset, the distance function (when 𝑘-medoids was used to obtain the input clustering C), the value 𝑘, the exact silhouette 𝑠 (C), and the runtime of the exact algorithm 𝑟𝑡 EX . 5.2.1

Medium datasets

We show some representative results in Figure 3, 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×, with respect to the exact approach, is the one relative to the CR dataset with 𝑑 can and 𝑘 = 10. Such behavior is not surprising since, for small 𝑘, it is likely that some clusters have roughly linear size, hence FS will need a quadratic number of distance computations to compute the exact 𝑎(𝑒)’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 𝑏 (𝑒), 𝑒 ∈ 𝑉 often corresponds to the cluster with the closest center in C to 𝑒. 29

gl-s gl-pps-f

gl-uni-f gl-pps-s, 𝑚 = 𝑛/100 GOW, 𝑘-means, 𝑘 = 2, 𝑠 (C) = 0.7777 𝑡 = 64

0.015

0.0075 0.0050

0.005

0.0025

100

101

100

Average runtime (sec) GOW, 𝑘-means, 𝑘 = 10, 𝑠 (C) = 0.5713 𝑡 = 64

0.03 0.02

100

0.00

101

0.005

100

101

Average runtime (sec) PH, 𝑘-means, 𝑘 = 2, 𝑠 (C) = 0.9397

𝑡 = 32

𝑡 = 64

0.005

0.00

100

10−1

Average runtime (sec)

10−1

PH, 𝑘-means, 𝑘 = 10, 𝑠 (C) = 0.3806 𝑡 = 64

0.020

0.06

100

𝑡 = 128

0.02

0.01

0.010 0.02

0.005

0.00

0.000

100

101

100

101

Average runtime (sec)

101

𝑡 = 32

GOW, 𝑘-means, 𝑘 = 20, 𝑠 (C) = 0.5404

100

101

100

𝑡 = 64

101

𝑡 = 128

0.008 0.006

0.02

0.004

0.01

101

0.002 0.000

0.00

101

102

𝑡 = 64

0.005

0.000

100

𝑡 = 128

0.010

0.005

100

101

100

Average runtime (sec) PH, 𝑘-means, 𝑘 = 20, 𝑠 (C) = 0.3343

𝑡 = 32

102

0.015

10−1

101

Average runtime (sec) PH, 𝑘-means, 𝑘 = 5, 𝑠 (C) = 0.72

0.010

0.06

0.015

0.04

100

Average runtime (sec)

0.02

0.000

101

Average error

𝑡 = 32

101

0.015

0.005

10−1

100

𝑡 = 32

0.010

0.01

0.010 0.005

𝑡 = 128

0.010

0.02

0.015

0.01

100

0.015

0.03

0.02

102

𝑡 = 128

0.020 0.02

0.03

0.010

0.01

0.00

𝑡 = 64

0.03

𝑡 = 128

0.015

0.02

0.01

GOW, 𝑘-means, 𝑘 = 5, 𝑠 (C) = 0.6778

0.020

0.03

0.04

101

gl-uni-s, 𝑚 = 𝑛/10

0.04

0.00

Average error

𝑡 = 32

Average error

0.0100

0.010

10−1

Average error

𝑡 = 32

0.0125

0.010

0.005

𝑡 = 128

Average error

0.020

gl-uni-s, 𝑚 = 𝑛/100 gl-uni-s, 𝑚 = 𝑛/20

Average error

𝑡 = 32

gl-pps-s, 𝑚 = 𝑛/20 gl-pps-s, 𝑚 = 𝑛/10

𝑡 = 64

101

𝑡 = 128

0.03 0.015

0.04

0.02

0.01

0.00

0.010

100

101

0.005

100

101

Average runtime (sec)

100

101

102

Figure 4: Methods comparison. The 𝑥 axis (in log scale) is associated with the average runtime, and the 𝑦 axis with the average error (and its standard deviation) over 10 independent runs. For each clustered dataset we considered different values of the expected sample size 𝑡 ∈ {32, 64, 128}.

Next, we note that the SIMPL method is often the fastest approach among all methods we compared (except on the SHU data clustered with 𝑘-means and 𝑘 = 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 𝑘-medoids and distance 𝑑 eucl the exact silhouette is 𝑠 (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 𝑘-medoids, 𝑑 can and 𝑘 = 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 𝑑 can and 𝑘 = 2), SHU (with 𝑑 cos and 𝑘 = 2), or IOT (with 𝑑 cos , 𝑘 = 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 𝑘-means and 𝑘 = 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 𝑘-means and distance 𝑑 eucl . We did not consider baseline FS for its poor scalability and the SIMPL method for its poor quality estimates (see Section 5.2.1). In addition, we used highly accurate estimates of 𝑠 (C) as ground truth, given the poor scalability of the exact computation (see Section 5.1.4). The results are reported in Figure 4. First, as expected, all algorithms exhibit increasing accuracy as the sample size 𝑡 grows. For example, on dataset GOW with 𝑘 = 5, the error with 𝑡 = 128 for gl-pps-f is less than half the error obtained with 𝑡 = 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 gluni-f respectively, using only a sample of size 𝑚 of the total elements 𝑛, we observe that a small sample (𝑚 = 𝑛/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× 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 Algorithm 1. 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 31

Average maximum error (𝜂¯𝐵 )

GOW, 𝑘 = 2, 𝑡 = 32 0.25 0.20 0.15 0.10 0.05 0.00 0.03

0.13

0.23

0.33

0.44

0.54

0.64

pps

0.74

0.85

0.95

GOW, 𝑘 = 10, 𝑡 = 64 0.4 0.3 0.2 0.1 0.0 -0.18 -0.08 0.03

0.13

0.23

0.33

0.44

0.54

0.64

0.74

Silhouette value (bucket)

IOT, 𝑑 can, 𝑘 = 5, 𝑡 = 64

IOT, 𝑑 manh, 𝑘 = 10, 𝑡 = 64

Average maximum error (𝜂¯𝐵 )

uni

0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00 -0.28 -0.18 -0.08 0.03 0.13 0.23 0.33 0.44 0.54 0.64 0.74 0.85 0.95

Silhouette value (bucket)

0.85

0.95

0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00 -0.38 -0.28 -0.18 -0.08 0.03 0.13 0.23 0.33 0.44 0.54 0.64 0.74 0.85 0.95

Silhouette value (bucket)

Figure 5: Plots of the average maximum error and its standard deviation obtained on all elements 𝑒 ∈ 𝑉 , grouped by bucket, for pps and uni, over four configurations. When we report the distance, it indicates that the dataset has been clustered with 𝑘-medoids. For ease of visualization we only display the average maximum error over non-empty buckets. 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 𝑚 compared to 𝑛—yielding a significant speedup on large datasets. Therefore, for large datasets, where the exact computation of 𝑠 (C) is impractical even for a single and fixed value of 𝑘, 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 𝑚 = 𝑛/10 and 𝑡 = 128 output very accurate silhouette estimates for all values of 𝑘 ∈ {2, 5, 10, 20} in less than five minutes and one minute, respectively. 32

5.3

Local silhouette estimation

We now validate our algorithm silh-pps-all that estimates all values of 𝑠 (𝑒), 𝑒 ∈ 𝑉 , 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 b 𝑠 U (𝑒), 𝑒 ∈ 𝑉 using uniform Poisson sampling over each cluster 𝐶 ∈ C (see Section 4.1.3) . For each dataset and value of 𝑘 considered in Section 4.2, we evaluate the accuracy of the local estimates in two settings. First, we consider a bucketing of the values of silhouette 𝑠 (𝑒), 𝑒 ∈ 𝑉 . That is, we assign each element 𝑒 ∈ 𝑉 to a bucket based on its silhouette value 𝑠 (𝑒) where each bucket 𝐵𝑖 , with 𝑖 = 1, . . . , 40, corresponds to elements 𝑒 with 𝑠 (𝑒) ∈ [𝐿𝑖 , 𝑈𝑖 ), where 𝐿𝑖 = −1 + (𝑖 − 1)/20 and 𝑈𝑖 = −1 + 𝑖/20.9 We use the buckets 𝐵𝑖 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 𝐶 1, . . . , 𝐶𝑘 in C. 5.3.1

Accuracy over buckets

For each bucket 𝐵𝑖 let 𝐸𝑖 ⊆ 𝑉 be the elements of 𝑉 with 𝑠 (𝑒) ∈ 𝐵𝑖 . We computed, for each independent run, the maximum error max𝑒∈𝐸𝑖 {|b 𝑠 (𝑒) − 𝑠 (𝑒)|} (where b 𝑠 (𝑒) is obtained either with pps or uni), and averaged the maximum error over ten independent runs, to obtain the average maximum error (𝜂¯𝐵𝑖 ) and its standard deviation, shown in Figure 5. Note that 𝜂¯𝐵𝑖 is a very precise metric to capture the local estimation performances, given by the fact that it considers the maximum error over each bucket 𝐵𝑖 , which can be very large, i.e., at most 2. From Figure 5 we first observe that both pps and uni approaches achieve smaller maximum error (𝜂¯𝐵𝑖 ) for buckets corresponding to elements with higher 𝑠 (𝑒). Such behavior is likely caused by the fact that for elements with high value of 𝑠 (𝑒), the difference between 𝑏 (𝑒) − 𝑎(𝑒) (see Equation (2)) is high, hence even rough estimates of terms 𝑊𝐶 𝑗 (𝑒) yield a good approximation of 𝑠 (𝑒). The previous observation also explains the larger error over elements with silhouette 𝑠 (𝑒) close to 0. Comparing the average error 𝜂¯𝐵𝑖 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 𝑑 can and 𝑘 = 5 the estimates reported by pps have half of the maximum error achieved by uni. Similarly, on the GOW data for 𝑘 = 2, 𝑡 = 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 𝑘 = 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 Theorem 1. We provide some additional results in Section D. 5.3.2

Accuracy over clusters

Next we investigate the local approximation error of pps and uni across each cluster of C. That is, for each cluster we computed both the maximum error max𝑒∈𝐶 {|b 𝑠 (𝑒) − 𝑠 (𝑒)|} and the average Í error (1/|𝐶 |) 𝑒∈𝐶 |b 𝑠 (𝑒) − 𝑠 (𝑒)| (where b 𝑠 (𝑒) is computed either with uni or pps). In Figure 6, 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 9The bucket corresponding to index 𝑖 = 40 also contains the right endpoint 𝐿

40 = 1.

uni

pps

0.15 0.10 0.05 0.00

0.4

0.4 0.3 0.2 0.1 0.0

Cluster index

Maximum error

0.20

GOW, 𝑘 = 10, 𝑡 = 64

0.5

Average error

Maximum error

Average error

GOW, 𝑘 = 5, 𝑡 = 32 0.15

0.10

0.05

0.00

Cluster index

0.050 0.025 0.000

0.20 0.15 0.10 0.05

Cluster index

0.15 0.10 0.05 0.00

0.00

Maximum error

0.075

0.20

Average error

Maximum error

Average error

0.100

0.1

Cluster index

PH, 𝑘 = 10, 𝑡 = 128

0.25

0.125

0.2

0.0

MT, 𝑘 = 2, 𝑡 = 64 0.150

0.3

Cluster index

0.25 0.20 0.15 0.10 0.05 0.00

0.30

Cluster index

Figure 6: Plots of the average and maximum errors within each cluster achieved by pps and uni, over ten independent runs, for various input configurations. Each dot represents a distinct run. significantly smaller average or maximum error compared to uni. As an example, on the MT dataset with 𝑘 = 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 𝑘 = 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 Section D. Summary for issue I2. We investigated the accuracy for the extremely challenging task of computing the local silhouette value of all elements 𝑒 ∈ 𝑉 (i.e., Problem 1). Our experiments show that our silh-pps-all algorithm provides highly-accurate estimates for all elements 𝑒 ∈ 𝑉 , with significantly reduced variance compared to the uni baseline. That is, the estimates of pps are more accurate independently of the actual silhouette values (Section 5.3.1), and independently of the cluster to which the elements belong (Section 5.3.2), which are desirable properties for practical applications.

5.4

Distributed implementation

As discussed in Section 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 silhpps-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 Section 4.3 (see Section C.2 for details of the platform used, and Section C.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 𝑘 ∈ {2, 5, 10, 15, 20}, and varying 𝑡 ∈ {512, 1024}. We considered our largest datasets from Table 2. For each configuration we perform five runs, and recorded the speedup over the sequential implementation. That is, let 𝑇𝑖 be the average time (over the five runs) to execute silh-pps-all with 𝑖 ∈ {1, 2, 4, 8, 16} threads, i.e., 𝑇1 is associated with the sequential implementation. We will study the average speedup 𝑇1 /𝑇𝑖 . 34

𝑘 =2

𝑘 =5

𝑘 =10

PH, 𝑡 = 512

PH, 𝑡 = 1024

2 1

4 2 1

1 2

Threads

MT, 𝑡 = 1024

Avg. speedup

4 2 1

1 2

𝑥 =𝑦

MT, 𝑡 = 512

𝑘 =20

Avg. speedup

𝑘 =15

4 2 1

1 2

Threads

1 2

Threads

Figure 7: Average parallel speedup over five independent runs by running silh-pps-all with a parallel implementation. Each configuration shows different values for the sample size 𝑡. The results are reported in Figure 7. We observe a remarkable speed-up over the sequential implementation of the parallel algorithm up to 8 threads on all 𝑘 values (except for 𝑘 = 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 𝑘 and 𝑡, since our implementation parallelizes the work over different clusters in C (see Section C.3). That is, for smaller values of 𝑘 and 𝑡 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 𝑡 leads to a higher speedup (e.g., on MT with 𝑡 = 1024 the speedup with 16 threads is about 14 for 𝑘 = 20). Remarkably, when using 𝑡 = 1024 on PH with 𝑘 ≥ 5 and 16 threads, we obtain a speedup of at least 10× 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 Table 4 in Section D, 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 C that closely matches the exact one (e.g., see an example of silhouette plot in Figure 1b). In this setting, we focus on BR, our smallest dataset, for ease of visualization,10 and for the fact that we can compute the exact values of 𝑠 (𝑒), for all 𝑒 ∈ 𝑉 .

10To avoid cluttered plots, it is good practice to visualize only a small subset of the elements of 𝑉 over silhouette

plots; even on large datasets.

−0.4

Cluster index

BR, 𝑘-medoids, 𝑘 = 2, 𝑑 can, 𝑡 = 16

0.0

−0.2

0.2

0.4

0.6

−0.4

Silhouette value 𝑠 (𝑒)

0.0

−0.2

0.2

0.4

0.6

Silhouette value b 𝑠 PPS (𝑒)

−0.4

BR, 𝑘-means, 𝑘 = 10, 𝑡 = 16

1 −0.50

−0.25

Cluster index

10 9 8 76 5 4 3 2

Cluster index

10 9 8 76 5 4 3 2

1 0.00

0.25

0.50

Silhouette value 𝑠 (𝑒)

0.75

1.00

−0.50

0.0

−0.2

0.2

0.4

0.6

Silhouette value b 𝑠 UNI (𝑒)

10 9 8 76 5 4 3 2 1

−0.25

0.00

0.25

0.50

Silhouette value b 𝑠 PPS (𝑒)

0.75

1.00

−0.50

−0.25

0.00

0.25

0.50

Silhouette value b 𝑠 UNI (𝑒)

0.75

1.00

Figure 8: Silhouette plot construction. (Left): Silhouette plot obtained from the exact algorithm. (Center): Silhouette plot obtained by silh-pps-all, where elements are ordered according to the Left plot. (Right): Silhouette plot obtained by uni, where elements are ordered according to the Left plot. For all plots, the dashed line represents the average value, computed over all elements. Figure 8 reports an example of two clusterings C, obtained respectively with 𝑘-medoids (top) and with 𝑘-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 𝑠 (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 𝑘

Let 𝑉 be a dataset and let 𝐾 = {𝑘 1, . . . , 𝑘 |𝐾 | } be a set containing different values of 𝑘 used to cluster 𝑉 . The global silhouette is often used to select the value of 𝑘 ∈ 𝐾 that yields the best clustering for 𝑉 . Specifically, for each 𝑘 ∈ 𝐾 we obtain a clustering C𝑘 (e.g., using 𝑘-means) and select C𝑘 ∗ as the best clustering for 𝑉 , where 𝑘 ∗ = arg max𝑘∈𝐾 𝑠 (C𝑘 ) [Shahapure and Nicholas, 2020]. In this subsection we compare the effectiveness of gl-pps-s versus gl-uni-s in detecting the best value of 𝑘 to be used for clustering. We focus on these two estimators, since, as argued before, 36

gl-pps-s, 𝑚 = 𝑛/10

CR, 𝑡 = 10

1.0

Cumulative accuracy

gl-uni-s, 𝑚 = 𝑛/10

0.9

0.8

9 10 11 12 13 14

WI, 𝑡 = 10

1.0 0.9 0.8 0.7 0.6

ℓ

9 10 11 12 13 14

ℓ

𝒌 is either 𝒌 pps or 𝒌 uni . For Figure 9: Plots of Cumulative accuracy(ℓ, b 𝒌), where ℓ ∈ [1, 14] and b each value of ℓ, the plots report the cumulative accuracy of each estimation method over 200 runs: values closer to 1, denote a higher accuracy for the identification of the best value of 𝑘. they are accurate and their computation is scalable (see summary for Issue I1 in Section 5.2). Thus, they can be used to compute high-quality approximations of 𝑠 (C𝑘 ) for many different candidate values of 𝑘. We considered two datasets for which we can compute 𝑠 (C𝑘 ) exactly and efficiently, namely CR and WI. We clustered each dataset with 𝑘-means, with different values of 𝑘 ∈ 𝐾, with 𝐾 = {2, . . . , 15}, and for each resulting clustering C𝑘 we computed the exact global silhouette 𝑠 (C𝑘 ). Then, for each of the two algorithms gl-pps-s and gl-uni-s, both run with 𝑚 = 𝑛/10 and 𝑡 = 10, we performed 𝑁 = 200 independent runs for each distinct value of 𝑘. Let b 𝑠 (C𝑘 ) 𝑗 be the estimate of either gl-pps-s or gl-uni-s on dataset 𝑉 clustered with value 𝑘 ∈ [2, 15] at iteration 𝑗 ∈ [1, 𝑁 ]. We let 𝒌 exact be the vector of the 14 values of 𝑘 ∈ 𝐾, sorted by non-increasing values of the corresponding global silhouette, that is, 𝑠 (C𝒌1exact ) ≥ 𝑠 (C𝒌2exact ) ≥ · · · ≥ 𝑠 (C𝒌14exact ). Similarly, we let 𝒌 pps ( 𝑗) (resp. 𝒌 uni ( 𝑗)) denote the vector of the 14 values of 𝑘 ∈ 𝐾, sorted by non-increasing values of the estimates of the corresponding global silhouette provided by gl-pps-s (resp. gl-uni-s) in the 𝑗-th run, for 𝑗 ∈ [1, 𝑁 ]. Note that, given a run 𝑗 ∈ [1, 𝑁 ] the best value of 𝑘 ∈ 𝐾 returned by an approximate method corresponds to 𝒌 1pps ( 𝑗), or 𝒌 1uni ( 𝑗). To assess the accuracy of gl-pps-s and gl-uni-s for selecting a reasonable value of 𝑘 ∈ 𝐾 with respect to the exact approach, we now define a cumulative accuracy function of each approximate method. Given ℓ ∈ [1, |𝐾 |], the cumulative accuracy evaluated at ℓ corresponds to the fraction of runs in which the approximate method selects as best value 𝑘 ∈ 𝐾 a value in the top-ℓ positions of 𝒌 exact (i.e., 𝑘 ∈ {𝒌 1exact, . . . , 𝒌 ℓexact }). More formally given the set of 𝑁 sorted vectors b 𝒌 ( 𝑗), where b 𝒌 ( 𝑗) is pps uni either 𝒌 ( 𝑗) or 𝒌 ( 𝑗), 𝑗 ∈ [1, 𝑁 ], the cumulative accuracy corresponds to, Cumulative accuracy(ℓ, b 𝒌)

𝑁 ℓ i 1 ∑︁ ∑︁ hb 1 𝒌 1 ( 𝑗) = 𝒌𝑖exact . 𝑁 𝑗=1 𝑖=1

Where 1[·] corresponds to the indicator function. Observe that, as a function of ℓ, the cumulative accuracy eventually hits the value of 1 for ℓ = |𝐾 |. 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 𝑘, according to 𝒌 exact . Results are reported in Figure 9, 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 37

the selection as best value of 𝑘 one of the top-3 values in 𝒌 exact on dataset CR, and one of the top-2 values in 𝒌 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 𝑘 in 𝒌 exact on the CR and less than 92% of the runs detect one of the top-2 values in 𝒌 exact on the WI dataset. We observe that the global accuracy is affected by how the exact silhouette values vary across the top-ℓ elements of 𝒌 exact (we show the ordering of 𝒌 exact in Table 5). For example, in the above experiments the silhouette values associated with the top-2 values of 𝒌 exact do not differ much, making the detection of the best 𝑘 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 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 𝑘 yielding higher 𝑠 (C𝑘 )) for a given dataset 𝑉 , improving over existing methods. Our experiments showcase the large applicability of the methods designed in this work for both local and global estimation scenarios.

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: 𝑖) the (approximate) construction of silhouette plots; and 𝑖𝑖) the selection of the best clustering granularity 𝑘 for a dataset. These are some open issues deserving further investigation: 𝑖) devising a more flexible, dynamic version of the algorithm where the desired accuracy can be incrementally refined by reusing (part of) the previously sampled elements; 𝑖𝑖) designing data-dependent sampling approaches, that adapt their sample complexity to the input datasets; 𝑖𝑖𝑖) extending our methods to other popular clustering validation measures, which we discuss in more detail in Section B.

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.

References C. C. Aggarwal and C. K. Reddy, editors. Data Clustering. Chapman and Hall/CRC, sep 2014. doi: 10.1201/9781315373515. F. Altieri, A. Pietracaprina, G. Pucci, and F. Vandin. Scalable Distributed Approximation of Internal Measures for Clustering Evaluation, pages 648–656. Society for Industrial and Applied Mathematics, Jan. 2021. ISBN 9781611976700. doi: 10.1137/1.9781611976700.73. O. Arbelaitz, I. Gurrutxaga, J. Muguerza, J. M. Pérez, and I. Perona. An extensive comparative study of cluster validity indices. Pattern Recognition, 46(1):243–256, Jan. 2013. ISSN 0031-3203. doi: 10.1016/j.patcog.2012.07.021. B. Bahmani, B. Moseley, A. Vattani, R. Kumar, and S. Vassilvitskii. Scalable k-means++. Proceedings of the VLDB Endowment, 5(7):622–633, mar 2012. doi: 10.14778/2180912.2180915. M. F. Balcan, S. Ehrlich, and Y. Liang. Distributed k-means and k-median clustering on general topologies. NIPS, June 2013. A. Barger and D. Feldman. k-means for streaming and distributed big sparse data. In Proceedings of the 2016 SIAM International Conference on Data Mining. Society for Industrial and Applied Mathematics, jun 2016. doi: 10.1137/1.9781611974348.39. P. Beame, P. Koutris, and D. Suciu. Communication steps for parallel query processing. J. ACM, 64(6):40:1–40:58, 2017. M. Ceccarello, C. Fantozzi, A. Pietracaprina, G. Pucci, and F. Vandin. Clustering uncertain graphs. Proceedings of the VLDB Endowment, 11(4):472–484, dec 2017. doi: 10.1145/3186728.3164143. M. Ceccarello, A. Pietracaprina, and G. Pucci. Solving k-center clustering (with outliers) in MapReduce and streaming, almost as accurately as sequentially. Proceedings of the VLDB Endowment, 12(7):766–778, mar 2019. doi: 10.14778/3317315.3317319. S. Chechik, E. Cohen, and H. Kaplan. Average distance queries through weighted samples in graphs and metric spaces: High scalability with tight statistical guarantees. In N. Garg, K. Jansen, A. Rao, and J. D. P. Rolim, editors, Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques, APPROX/RANDOM 2015, August 24-26, 2015, 39

Princeton, NJ, USA, volume 40 of LIPIcs, pages 659–679. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2015. doi: 10.4230/LIPIcs.APPROX-RANDOM.2015.659. E. Cohen, S. Chechik, and H. Kaplan. Clustering small samples with quality guarantees: Adaptivity with one2all PPS. In S. A. McIlraith and K. Q. Weinberger, editors, Proceedings of the Thirty-Second AAAI Conference on Artificial Intelligence, (AAAI-18), the 30th innovative Applications of Artificial Intelligence (IAAI-18), and the 8th AAAI Symposium on Educational Advances in Artificial Intelligence (EAAI-18), New Orleans, Louisiana, USA, February 2-7, 2018, pages 2884–2891. AAAI Press, 2018. URL https://www.aaai.org/ocs/index.php/AAAI/AAAI18/ paper/view/17084. A. Czumaj, G. Gao, M. Ghaffari, and S. H. Jiang. Fully scalable MPC algorithms for euclidean k-center. In Proceedings of the 52nd International Colloquium on Automata, Languages, and Programming, ICALP 2025, volume 334 of LIPIcs, pages 64:1–64:20, 2025. J. Dean and S. Ghemawat. MapReduce. Communications of the ACM, 51(1):107–113, jan 2008. doi: 10.1145/1327452.1327492. D.-T. Dinh, T. Fujinami, and V.-N. Huynh. Estimating the Optimal Number of Clusters in Categorical Data Clustering by Silhouette Coefficient, pages 1–17. Springer Singapore, 2019. ISBN 9789811512094. doi: 10.1007/978-981-15-1209-4_1. D. P. Dubhashi and A. Panconesi. Concentration of Measure for the Analysis of Randomized Algorithms. Cambridge University Press, Cambridge, UK, 2009. ISBN 978-0-521-88427-3. A. Dudek. Silhouette Index as Clustering Evaluation Tool, pages 19–33. Springer International Publishing, 2020. ISBN 9783030523480. doi: 10.1007/978-3-030-52348-0_2. A. Ene, S. Im, and B. Moseley. Fast clustering using MapReduce. In Proceedings of the 17th ACM SIGKDD international conference on Knowledge discovery and data mining - KDD '11. ACM Press, 2011. doi: 10.1145/2020408.2020515. D. Feldman and M. Langberg. A unified framework for approximating and clustering data. In Proceedings of the 43rd annual ACM symposium on Theory of computing - STOC '11. ACM Press, 2011. doi: 10.1145/1993636.1993712. G. Frahling and C. Sohler. A Fast k-Means Implementatino Using Coresets. International Journal of Computational Geometry & Applications, 18(06):605–625, dec 2008. doi: 10.1142/ s0218195908002787. J. Han, M. Kamber, and J. Pei. Data mining: concepts and techniques. Elsevier, 2011. doi: 10.1016/c2009-0-61819-5. B. A. Hassan, N. B. Tayfor, A. A. Hassan, A. M. Ahmed, T. A. Rashid, and N. N. Abdalla. From a-to-z review of clustering validation indices. Neurocomputing, 601:128198, Oct. 2024. ISSN 0925-2312. doi: 10.1016/j.neucom.2024.128198. C. Hennig, M. Meila, F. Murtagh, and R. Rocci. Handbook of Cluster Analysis. CRC PRESS, S.l, 2015. ISBN 9780367570408. 40

M. S. Hossain and R. A. Angryk. GDClust: A graph-based document clustering technique. In Seventh IEEE International Conference on Data Mining Workshops (ICDMW 2007). IEEE, oct 2007. doi: 10.1109/icdmw.2007.104. E. Hruschka, L. de Castro, and R. Campello. Evolutionary algorithms for clustering geneexpression data. In Fourth IEEE International Conference on Data Mining (ICDM'04). IEEE. doi: 10.1109/icdm.2004.10073. A. M. Ikotun, A. E. Ezugwu, L. Abualigah, B. Abuhaija, and J. Heming. K-means clustering algorithms: A comprehensive review, variants analysis, and advances in the era of big data. Information Sciences, 622:178–210, Apr. 2023. ISSN 0020-0255. doi: 10.1016/j.ins.2022.11.139. S. Im, R. Kumar, S. Lattanzi, B. Moseley, and S. Vassilvitskii. Massively parallel computation: Algorithms and applications. Foundations and Trends® in Optimization, 5(4):340–417, 2023. H. J. Karloff, S. Suri, and S. Vassilvitskii. A model of computation for mapreduce. In M. Charikar, editor, Proceedings of the Twenty-First Annual ACM-SIAM Symposium on Discrete Algorithms, SODA 2010, Austin, Texas, USA, January 17-19, 2010, pages 938–948. SIAM, 2010. doi: 10.1137/ 1.9781611973075.76. G. Karypis and V. Kumar. A fast and high quality multilevel scheme for partitioning irregular graphs. SIAM Journal on Scientific Computing, 20(1):359–392, Jan. 1998. ISSN 1095-7197. doi: 10.1137/s1064827595287997. L. Kaufmann and P. Rousseeuw. Clustering by means of medoids. Data Analysis based on the L1-Norm and Related Methods, pages 405–416, 01 1987. J. Kleinberg. An impossibility theorem for clustering. Advances in neural information processing systems, 15, 2002. H. Lai, T. Huang, B. Lu, S. Zhang, and R. Xiaog. Silhouette coefficient-based weighting k-means algorithm. Neural Computing and Applications, 37(5):3061–3075, Dec. 2024. ISSN 1433-3058. doi: 10.1007/s00521-024-10706-0. L. Lenssen and E. Schubert. Medoid silhouette clustering with automatic cluster number selection. Information Systems, 120:102290, Feb. 2024. ISSN 0306-4379. doi: 10.1016/j.is.2023.102290. J. Leskovec, A. Rajaraman, and J. D. Ullman. Mining of Massive Datasets. Cambridge University Press, 2014. doi: 10.1017/cbo9781139924801. L. Liu, R. Jin, C. C. Aggarwal, and Y. Shen. Reliable clustering on uncertain graphs. In M. J. Zaki, A. Siebes, J. X. Yu, B. Goethals, G. I. Webb, and X. Wu, editors, 12th IEEE International Conference on Data Mining, ICDM 2012, Brussels, Belgium, December 10-13, 2012, pages 459–468. IEEE Computer Society, 2012. doi: 10.1109/ICDM.2012.11. Y. Liu, Z. Li, H. Xiong, X. Gao, and J. Wu. Understanding of internal clustering validation measures. In 2010 IEEE International Conference on Data Mining. IEEE, Dec. 2010. doi: 10.1109/ icdm.2010.35. 41

R. Lletı, M. Ortiz, L. Sarabia, and M. Sánchez. Selecting variables for k-means cluster analysis by using a genetic algorithm that optimises the silhouettes. Analytica Chimica Acta, 515(1): 87–100, July 2004. ISSN 0003-2670. doi: 10.1016/j.aca.2003.12.020. S. Lloyd. Least squares quantization in PCM. IEEE Transactions on Information Theory, 28(2): 129–137, mar 1982. doi: 10.1109/tit.1982.1056489. L. Lovmar, A. Ahlford, M. Jonsson, and A.-C. Syvänen. Silhouette scores for assessment of snp genotype clusters. BMC Genomics, 6(1), Mar. 2005. ISSN 1471-2164. doi: 10.1186/ 1471-2164-6-35. G. Malkomes, M. J. Kusner, W. Chen, K. Q. Weinberger, and B. Moseley. Fast distributed kcenter clustering with outliers on massive data. In C. Cortes, N. D. Lawrence, D. D. Lee, M. Sugiyama, and R. Garnett, editors, Advances in Neural Information Processing Systems 28: Annual Conference on Neural Information Processing Systems 2015, December 7-12, 2015, Montreal, Quebec, Canada, pages 1063–1071, 2015. URL https://proceedings.neurips.cc/ paper/2015/hash/8fecb20817b3847419bb3de39a609afe-Abstract.html. A. Mazzetto, A. Pietracaprina, and G. Pucci. Accurate mapreduce algorithms for k-median and k-means in general metric spaces. In P. Lu and G. Zhang, editors, 30th International Symposium on Algorithms and Computation, ISAAC 2019, December 8-11, 2019, Shanghai University of Finance and Economics, Shanghai, China, volume 149 of LIPIcs, pages 34:1–34:16. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2019. doi: 10.4230/LIPIcs.ISAAC.2019.34. M. Mitzenmacher and E. Upfal. Probability and Computing. Cambridge University Press, Cambridge, UK, 2nd edition, 2017. ISBN 9781316651124. D. Moulavi, P. A. Jaskowiak, R. J. G. B. Campello, A. Zimek, and J. Sander. Density-based clustering validation. In M. J. Zaki, Z. Obradovic, P. Tan, A. Banerjee, C. Kamath, and S. Parthasarathy, editors, Proceedings of the 2014 SIAM International Conference on Data Mining, Philadelphia, Pennsylvania, USA, April 24-26, 2014, pages 839–847. SIAM, 2014. doi: 10.1137/1.9781611973440. 96. R. T. Ng and J. Han. Efficient and effective clustering methods for spatial data mining. In J. B. Bocca, M. Jarke, and C. Zaniolo, editors, VLDB’94, Proceedings of 20th International Conference on Very Large Data Bases, September 12-15, 1994, Santiago de Chile, Chile, pages 144–155. Morgan Kaufmann, 1994. URL http://www.vldb.org/conf/1994/P144.PDF. A. Pietracaprina, G. Pucci, M. Riondato, F. Silvestri, and E. Upfal. Space-round tradeoffs for mapreduce computations. In U. Banerjee, K. A. Gallivan, G. Bilardi, and M. Katevenis, editors, International Conference on Supercomputing, ICS’12, Venice, Italy, June 25-29, 2012, pages 235– 244. ACM, 2012. doi: 10.1145/2304576.2304607. A. Pranjal and M. F. Balcan. Center based clustering: A foundational perspective, 2013. P. J. Rousseeuw. Silhouettes: A graphical aid to the interpretation and validation of cluster analysis. Journal of Computational and Applied Mathematics, 20:53–65, nov 1987. doi: 10.1016/ 0377-0427(87)90125-7. 42

E. Schubert. Stop using the elbow criterion for k-means and how to choose the number of clusters instead. ACM SIGKDD Explorations Newsletter, 25(1):36–42, June 2023. ISSN 1931-0153. doi: 10.1145/3606274.3606278. E. Schubert and L. Lenssen. Fast k-medoids clustering in rust and python. Journal of Open Source Software, 7(75):4183, 2022. E. Schubert and P. J. Rousseeuw. Fast and eager 𝑘-medoids clustering: O (𝑘) runtime improvement of the pam, clara, and clarans algorithms. Information Systems, 101:101804, Nov. 2021. ISSN 0306-4379. doi: 10.1016/j.is.2021.101804. T. Sellam, R. Cijvat, R. Koopmanschap, and M. L. Kersten. Blaeu: Mapping and navigating large tables with cluster analysis. Proc. VLDB Endow., 9(13):1477–1480, 2016. doi: 10.14778/3007263. 3007288. URL http://www.vldb.org/pvldb/vol9/p1477-sellam.pdf. K. R. Shahapure and C. Nicholas. Cluster quality analysis using silhouette score. In 2020 IEEE 7th International Conference on Data Science and Advanced Analytics (DSAA), pages 747–748. IEEE, Oct. 2020. doi: 10.1109/dsaa49011.2020.00096. P.-N. Tan, M. Steinbach, and V. Kumar. Introduction to Data Mining. Pearson Addison Wesley, Boston, 2005. ISBN 9780321321367. C. Tomasini, L. R. Emmendorfer, E. N. Borges, and K. S. Machado. A methodology for selecting the most suitable cluster validation internal indices. In S. Ossowski, editor, Proceedings of the 31st Annual ACM Symposium on Applied Computing, Pisa, Italy, April 4-8, 2016, pages 901–903. ACM, 2016. doi: 10.1145/2851613.2851885. M. Van der Laan, K. Pollard, and J. Bryan. A new partitioning around medoids algorithm. Journal of Statistical Computation and Simulation, 73(8):575–584, Aug. 2003. ISSN 1563-5163. doi: 10.1080/0094965031000136012. F. Wang, H. Franco-Penya, J. D. Kelleher, J. Pugh, and R. J. Ross. An analysis of the application of simplified silhouette to the evaluation of k-means clustering validity. In P. Perner, editor, Machine Learning and Data Mining in Pattern Recognition - 13th International Conference, MLDM 2017, New York, NY, USA, July 15-20, 2017, Proceedings, volume 10358 of Lecture Notes in Computer Science, pages 291–305. Springer, 2017. doi: 10.1007/978-3-319-62416-7\_21. URL https://doi.org/10.1007/978-3-319-62416-7_21. C. Wiwie, J. Baumbach, and R. Röttger. Comparing the performance of biomedical clustering methods. Nature Methods, 12(11):1033–1038, sep 2015. doi: 10.1038/nmeth.3583. H. Xiong and Z. Li. Clustering validation measures. In C. C. Aggarwal and C. K. Reddy, editors, Data Clustering: Algorithms and Applications, pages 571–606. CRC Press, 2013. D. Xu and Y. Tian. A comprehensive survey of clustering algorithms. Annals of Data Science, 2 (2):165–193, June 2015. ISSN 2198-5812. doi: 10.1007/s40745-015-0040-1.

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 [Mitzenmacher and Upfal, 2017, Theorem 4.12], and is crucial for the analysis of our estimators. Theorem 6 (Hoeffding’s bound). Let 𝑋 1, . . . , 𝑋𝑚 be independent random variables taking values in [𝑎, 𝑏] ⊆ R such that for all 𝑗 ∈ [1, 𝑚], E[𝑋 𝑗 ] = 𝜇. Then " # 1 ∑︁ 2𝑚𝜀 2 . P 𝑋 𝑗 − 𝜇 ≥ 𝜀 ≤ 2 exp − 𝑚 𝑗 (𝑏 − 𝑎) 2 The next results (belonging to the large farmily of bounds known as Chernoff-Hoeffding bounds) are used in the proof of Theorem 1. The form of the bounds that we use is from the textbook by Dubhashi and Panconesi [2009, Theorem 1.1]. Theorem 7 (Chernoff-Hoeffding bounds). Let 𝑋 1, . . . , 𝑋𝑚 be independent random variables taking Í values in [0, 1] and let 𝑋 = 𝑚 𝑖=1 𝑋𝑖 . Then for any 𝜀 ∈ (0, 1) it holds that 2 2 𝜀 𝜀 and P [𝑋 ≤ (1 − 𝜀)E[𝑋 ]] ≤ exp − E[𝑋 ] . P [𝑋 ≥ (1 + 𝜀)E[𝑋 ]] ≤ exp − E[𝑋 ] 3 2

Estimating other clustering validation measures

Our strategy that estimates the terms 𝑊𝐶 𝑗 (𝑒), 𝑗 ∈ [𝑘], 𝑒 ∈ 𝑉 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 𝑉 . 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 C Tan et al. [2005]. More precisely, consider a 𝑘-clustering C = {𝐶 1, . . . , 𝐶𝑘 } and define the cohesion and separation of C as Í𝑘 Í ′ ′′ 1 𝑗=1 𝑒 ′,𝑒 ′′ ∈𝐶 𝑗 𝑑 (𝑒 , 𝑒 ) , Coh(C) = Í𝑘 |𝐶 𝑗 | 2 𝑗=1 2 Í Í Í ′ ′′ 1≤𝑗1 < 𝑗2 ≤𝑘 𝑒 ′ ∈𝐶 𝑗1 𝑒 ′′ ∈𝐶 𝑗2 𝑑 (𝑒 , 𝑒 ) Í Sep(C) = , 1≤𝑗1 < 𝑗2 ≤𝑘 (|𝐶 𝑗−1 ||𝐶 𝑗2 |) respectively. These measures have also been used to evaluate the average cluster reliability on networks, where distances correspond to connection probabilities [Ceccarello et al., 2017, Liu et al., 2012]. We can rewrite the above measures in terms of the sums 𝑊𝐶 (𝑒) defined in Equation (4), as follows Í Í Í𝑘 Í 1≤𝑗1 < 𝑗2 ≤𝑘 𝑒∈𝐶 𝑗1 𝑊𝐶 𝑗2 (𝑒) 1 𝑗=1 𝑒∈𝐶 𝑗 𝑊𝐶 𝑗 (𝑒) Coh(C) = and Sep(C) = Í . Í𝑘 |𝐶 𝑗 | 2 1≤𝑗1 < 𝑗2 ≤𝑘 (|𝐶 𝑗1 ||𝐶 𝑗2 |) 𝑗=1

Clearly, highly accurate approximations of terms 𝑊𝐶 𝑗 (𝑒), 𝑗 ∈ [𝑘], 𝑒 ∈ 𝑉 with low relative errors yield approximations with low relative error for Coh(C) and Sep(C). Specifically, let d d Coh(C) and Sep(C) as the approximation of Coh(C) and Sep(C) respectively. Where each b𝐶 (𝑒) computed approximation is obtained by substituting the term 𝑊𝐶 (𝑒) with the value 𝑊 within silh-pps-all. The following result is an immediate consequence of Lemma 1.

d Theorem 8. Let 𝑉 be a dataset of 𝑛 elements, and let C be a 𝑘-clustering of 𝑉 . Let Coh(C) b𝐶 (𝑒) d and Sep(C) be the approximations to Coh(C) and Sep(C) respectively based on the values 𝑊 computed within silh-pps-all for parameters 0 < 𝜀, 𝛿 < 1, and for a suitable choice of constant 𝑐 > 0 in the definition of the sample size 𝑡. Then with probability at least 1 − 𝛿 d Coh(C) − Coh(C) ≤𝜀 Coh(C)

Reproducibility

C.1

Datasets

and

d Sep(C) − Sep(C) ≤𝜀 . Sep(C)

We consider various datasets from the UCI and SNAP repositories, that are publicly available online.11 We provide a script to download and process each dataset with our code.12 For the 𝑘-means objective we used the implementation from scikit-learn,13 which uses Lloyd’s algorithm coupled with an initialization obtained through 𝑘-means++. For the 𝑘-medoids objective dataset we used the state-of-the-art approach by Schubert and Lenssen [2022].14

C.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 𝑘-means objective, we used the default implementation from the Scikit-learn library,15 while for the 𝑘-medoids objective we used the 𝑘-medoids library [Schubert and Rousseeuw, 2021].16 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.

C.3

Parallel implementation

To parallelize our algorithm silh-pps-all we rely on the OpenMP library. Our implementation adopts the following parallelization strategy: 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/

Table 3: Algorithms compared in our experimental evaluation in Section 5.2, and their required parameters: 𝑡 is used in Phase 1 of our Algorithm 1; 𝑚 is used to evaluate the estimators in Equations (9) and (10); while we denote with 𝑚′ the value of the sample size used to evaluate Equation (8). Both parameters 𝑡 and 𝑚 apply to pps and its uniform Poisson-sampling variant uni (see Section 4.1.3). Algorithm name

𝑡

𝑚 𝑚′

gl-s ✗ ✗ gl-pps-f, gl-uni-f ✓ ✗ gl-pps-s, gl-uni-s ✓ ✓

✓ ✗ ✗

Table 4: Average runtime (in seconds) required by the sequential execution of silh-pps-all (i.e., 𝑇1 ) under the setting of Section 5.4. Configuration 𝑘 = 2 𝑘 = 5 𝑘 = 10 𝑘 = 15 𝑘 = 20 MT, 𝑡 = 512 MT, 𝑡 = 1024 PH, 𝑡 = 512 PH, 𝑡 = 1024

35.8 72.8 44.4 87.7

72.5 223.9 109.6 354.6

254.3 525.7 356.8 756.3

392.6 702.9 561.0 1075.2

506.8 984.8 736.9 1427.7

• We parallelize the for loops in Line 1 of Algorithm 1, inspired by Round 1 and Round 2 of our Mapreduce algorithm; • We parallelize the for loop in Line 1 of Algorithm 1, 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 Section 5.4.

Additional results

In this section we complement our extensive experimental evaluation, by presenting additional results to those of Section 5. • In Figure 10 we show, for all the missing large dataset of Table 2. The results of the experiments in Section 5.2.2 that compare the global estimates of the various randomized algorithms. • In Figure 11 we show additional results for the accuracy of pps and uni approaches over various buckets—extending results from Section 5.3.1. • In Figure 12 we show additional results on the average and maximum error on each cluster C of pps and uni for the setting of Section 5.3.2. 46

Table 5: Values of 𝑘 ordered by 𝑠 (C𝑘 ) (ordering over 𝒌 exact ). The setting is from Section 5.5.2. CR 𝑘 𝒌 1exact 𝒌 2exact 𝒌 3exact 𝒌 4exact 𝒌 5exact 𝒌 6exact 𝒌 7exact 𝒌 8exact 𝒌 9exact exact 𝒌 10 exact 𝒌 11 exact 𝒌 12 exact 𝒌 13 exact 𝒌 14

𝑘 =2 𝑘 =4 𝑘 =3 𝑘 =5 𝑘 =7 𝑘 =6 𝑘 =8 𝑘 = 10 𝑘 =9 𝑘 = 11 𝑘 = 15 𝑘 = 14 𝑘 = 12 𝑘 = 13

WI 𝑠 (C𝑘 ) 0.552 0.463 0.451 0.419 0.391 0.387 0.38 0.364 0.343 0.336 0.333 0.314 0.312 0.304

𝑘 𝑘 =2 𝑘 =3 𝑘 =4 𝑘 =5 𝑘 =6 𝑘 =7 𝑘 =8 𝑘 =9 𝑘 = 10 𝑘 = 11 𝑘 = 12 𝑘 = 14 𝑘 = 15 𝑘 = 13

𝑠 (C𝑘 ) 0.51 0.504 0.45 0.418 0.392 0.367 0.357 0.342 0.334 0.329 0.326 0.325 0.325 0.323

• In Table 4 we show the average sequential runtime required by silh-pps-all considering the setting of Section 5.4. • In Table 5 we show the values of 𝑠 (C𝑘 ) for each 𝑘 ∈ 𝒌 under the setting of Section 5.5.

gl-s gl-pps-f

BioKDD, 𝑘-means, 𝑘 = 2, 𝑠 (C) = 0.5579 0.04

0.050 0.025

0.00

10−1

100

𝑡 = 32

BioKDD, 𝑘-means, 𝑘 = 10, 𝑠 (C) = 0.1902

10−1

𝑡 = 64

0.04

0.00

100

10−1

100

101

MT, 𝑘-means, 𝑘 = 2, 𝑠 (C) = 0.4212 𝑡 = 64

0.02

0.01

100

10−1

0.010 0.005

𝑡 = 32

100

10−1

Average runtime (sec)

MT, 𝑘-means, 𝑘 = 10, 𝑠 (C) = 0.4712 𝑡 = 64

0.020

0.04

0.010

0.02

0.00

101

𝑡 = 32

RNA, 𝑘-means, 𝑘 = 2, 𝑠 (C) = 0.5441

0.03

100

101

𝑡 = 64

0.02

0.01

𝑡 = 32

0.000

100

RNA, 𝑘-means, 𝑘 = 10, 𝑠 (C) = 0.4531 𝑡 = 64

0.005

0.000

10−1

100

0.01 0.00

0.00

100

101

𝑡 = 64

0.02

0.010

0.01

0.005

100

10−1

101

100

Average runtime (sec)

101

𝑡 = 128

10−1

100

101

0.000

MT, 𝑘-means, 𝑘 = 20, 𝑠 (C) = 0.4309 𝑡 = 64

0.015

0.0075

0.010

0.0050

0.02 0.01

0.005

100

101

Average runtime (sec) 𝑡 = 64

100

0.03

0.015

0.010

0.005

0.01

0.005

100

𝑡 = 32

RNA, 𝑘-means, 𝑘 = 20, 𝑠 (C) = 0.4801

10−1

100

0.000

10−1

Average runtime (sec)

0.03

𝑡 = 64

0.02

𝑡 = 128

100

10−1

𝑡 = 128

0.005

0.01

100

102

0.010

0.02

0.01

101

0.015

0.020

0.03

𝑡 = 128

RNA, 𝑘-means, 𝑘 = 5, 𝑠 (C) = 0.4937

0.02

101

0.0025

100

0.04

10−1

100

Average runtime (sec)

0.000

10−1

100

Average runtime (sec)

0.0100

0.04

0.005

0.01

0.02

0.020

𝑡 = 128

0.010

0.02

0.03

100

0.015

0.03

0.04

101

𝑡 = 128

MT, 𝑘-means, 𝑘 = 5, 𝑠 (C) = 0.5186

0.00

10−1

Average runtime (sec)

0.020

0.04

101

100

𝑡 = 32

0.005

10−1

100

10−1

0.03

0.04

𝑡 = 128

0.010

0.00

10−1

𝑡 = 64

0.06

0.02

0.00

101

0.015

0.03

0.02

10−2

100

Average runtime (sec) 0.04

100

BioKDD, 𝑘-means, 𝑘 = 20, 𝑠 (C) = 0.1297

0.04

10−1

0.00

0.000

100

10−1

𝑡 = 32

0.01

100

10−1

0.005

10−1

0.05

𝑡 = 128

0.005

0.01

0.005

Average runtime (sec)

0.10

101

0.02

0.015

0.03

0.010

0.00

0.000

100

10−1

0.010

0.06

0.015

0.00

100

𝑡 = 32

Average error

0.02

0.01

𝑡 = 128

0.020

0.04

0.02

0.00

101

Average error

0.06

Average error

100

Average runtime (sec)

𝑡 = 32

𝑡 = 128

0.015

𝑡 = 32

0.01

0.02

0.03

𝑡 = 128

0.02

𝑡 = 64

0.020

10−1

0.03

0.06

BioKDD, 𝑘-means, 𝑘 = 5, 𝑠 (C) = 0.311

𝑡 = 32

100

10−1

Average runtime (sec)

Average error

10−2

gl-uni-s, 𝑚 = 𝑛/10

0.00

Average error

0.02

0.00

gl-uni-s, 𝑚 = 𝑛/100 gl-uni-s, 𝑚 = 𝑛/20

0.04

0.02

0.08

Average error

𝑡 = 128

Average error

𝑡 = 64

0.075

0.000

Average error

gl-pps-s, 𝑚 = 𝑛/20 gl-pps-s, 𝑚 = 𝑛/10

Average error

𝑡 = 32

gl-uni-f gl-pps-s, 𝑚 = 𝑛/100

100

101

0.000

Average runtime (sec)

100

101

Figure 10: Methods comparison. The 𝑥 axis is associated with the average runtime, and the 𝑦 axis 48 with the average error (and its standard deviation) over 10 independent runs. For each clustered dataset we considered different values of the expected sample size 𝑡 ∈ {32, 64, 128}, as illustrated.

Average maximum error (𝜂¯𝐵 )

IOT, 𝑑 eucl, 𝑘 = 5, 𝑡 = 128 0.6 0.5 0.4 0.3 0.2 0.1 0.0 -0.59-0.49-0.38-0.28-0.18-0.08 0.03 0.13 0.23 0.33 0.44 0.54 0.64 0.74 0.85 0.95

0.8

0.6

0.4

0.2

0.0 -0.38 -0.28 -0.18 -0.08 0.03 0.13 0.23 0.33 0.44 0.54 0.64 0.74 0.85 0.95

Silhouette value (bucket)

PH, 𝑘 = 15, 𝑡 = 128

PH, 𝑘 = 2, 𝑡 = 128

0.3

0.2

0.1

0.0 -0.18

IOT, 𝑑 manh, 𝑘 = 5, 𝑡 = 128

Silhouette value (bucket)

0.4

-0.28

pps

Average maximum error (𝜂¯𝐵 )

uni

-0.08

0.03

0.13

0.23

0.33

0.44

0.54

0.64

0.74

0.08 0.06 0.04 0.02 0.00 0.64

MT, 𝑘 = 2, 𝑡 = 64 0.25 0.20 0.15 0.10 0.05 0.00 -0.18

-0.08

0.03

0.13

0.23

0.33

0.44

0.54

0.74

0.85

0.95

Silhouette value (bucket)

Average maximum error (𝜂¯𝐵 )

Silhouette value (bucket)

0.64

0.74

Silhouette value (bucket)

MT, 𝑘 = 20, 𝑡 = 64 0.4 0.3 0.2 0.1 0.0 -0.18

-0.08

0.03

0.13

0.23

0.33

0.44

0.54

0.64

0.74

0.85

Silhouette value (bucket)

Figure 11: Plots of the average maximum error and its standard deviation obtained on all elements 𝑒 ∈ 𝑉 , grouped by bucket, for pps and uni, over four configurations. When we report the distance, it indicates that the dataset has been clustered with 𝑘-medoids. For ease of visualization we only display the average maximum error over non-empty buckets.

uni

pps PH, 𝑘 = 5, 𝑡 = 128

PH, 𝑘 = 2, 𝑡 = 64

0.0100 0.0075 0.0050 0.0025 0.0000

0.08 0.06 0.04 0.02

Cluster index

0.08 0.06 0.04 0.02 0.00

0.00

0.30

0.10

Maximum error

0.0125

Average error

Maximum error

Average error

0.0150

Cluster index

0.10 0.05 0.00

Cluster index

0.2

0.1

Cluster index

0.15 0.10 0.05 0.00

Cluster index

0.30 0.25 0.20 0.15 0.10 0.05 0.00

Cluster index

RNA, 𝑘 = 5, 𝑡 = 64

Cluster index

RNA, 𝑘 = 15, 𝑡 = 128

0.20

0.10 0.05 0.00

0.3

Average error

0.15

0.2

0.1

0.0

Cluster index

Maximum error

0.12

Maximum error

Average error

0.05

0.20 0.3

0.0

0.10

Maximum error

0.15

BioKDD, 𝑘 = 10, 𝑡 = 128 Average error

Maximum error

Average error

0.20

0.00

BioKDD, 𝑘 = 2, 𝑡 = 64 0.25

0.25

0.10 0.08 0.06 0.04 0.02 0.00

Cluster index

0.25 0.20 0.15 0.10 0.05 0.00

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

Cluster index

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

Cluster index

Figure 12: For each cluster in C we show the average and maximum errors achieved by silh-ppsall (pps) and the uniform sampling based approach from Section 4.1.3 (uni). Each dot represent an independent run.

Record · ID 332491 · SHA-256 879939d37d587726
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.