{'fraction_non_alphanumeric': 0.038203283842197885, 'fraction_numerical': 0.013395069111038972, 'mean_word_length': 4.381574108127988, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
ParChain: A Framework for Parallel Hierarchical Agglomerative Clustering using Nearest-Neighbor Chain

This paper studies the hierarchical clustering problem, where the goal is to produce a dendrogram that represents clusters at varying scales of a data set. We propose the ParChain framework for designing parallel hierarchical agglomerative clustering (HAC) algorithms, and using the framework we obtain novel parallel algorithms for the complete linkage, average linkage, and Ward's linkage criteria. Compared to most previous parallel HAC algorithms, which require quadratic memory, our new algorithms require only linear memory, and are scalable to large data sets. ParChain is based on our parallelization of the nearest-neighbor chain algorithm, and enables multiple clusters to be merged on every round. We introduce two key optimizations that are critical for efficiency: a range query optimization that reduces the number of distance computations required when finding nearest neighbors of clusters, and a caching optimization that stores a subset of previously computed distances, which are likely to be reused. Experimentally, we show that our highly-optimized implementations using 48 cores with two-way hyper-threading achieve 5.8--110.1x speedup over state-of-the-art parallel HAC algorithms and achieve 13.75--54.23x self-relative speedup. Compared to state-of-the-art algorithms, our algorithms require up to 237.3x less space. Our algorithms are able to scale to data set sizes with tens of millions of points, which existing algorithms are not able to handle.


INTRODUCTION
Clustering is an unsupervised machine learning method that has been widely used in many fields including computational biology, computer vision, and finance to discover structures in a data set [2,6,24,32,39,46,65,69]. To group similar objects at all resolutions, a hierarchical clustering can be used to produce a tree that represents clustering results at different scales. The resulting hierarchical cluster structure is called a dendrogram, which is a tree representing the agglomeration of clusters, as shown in Figure 1 There is a rich literature on designing hierarchical agglomerative clustering (HAC) algorithms [46]. Unfortunately, exact HAC algorithms usually require Ω( 2 ) work, since the distances between all pairs of points have to be computed. To accelerate exact HAC algorithms due to their significant computational cost, there have been several parallel exact HAC algorithms proposed in the literature [23,31,33,41,42,58,67,70], but most of them maintain a distance matrix, which requires quadratic memory, making them unscalable to large data sets. The only parallel exact algorithm that works for the metrics that we consider and uses subquadratic space is by Zhang et al. [70], but it has not been shown to scale to large data sets. In this paper, we propose the ParChain framework for designing parallel exact HAC algorithms that use linear memory, based on the classic nearest-neighbor chain algorithm.
The nearest-neighbor chain (NNC) algorithm [45] is a popular algorithm that can be used for a wide range of HAC metrics [5,19,28,34,53]. A nearest-neighbor chain (NNC) is a linked list of nodes, where each node represents a cluster and all except at most one node have a pointer to its nearest neighbor (its successor). The chain can start from an arbitrary cluster. If a node does not have a pointer, its nearest neighbor is not yet computed, and this node is called a terminal node. If we follow the pointers on the nodes, we obtain a "chain" of clusters, which either terminates at a terminal node, or at a reciprocal nearest neighbor (R-NN) pair, which is a pair of clusters that are each other's nearest neighbor. The sequential NNC algorithm [5,19] works by iteratively adding a node to a single chain through finding the nearest neighbor of the terminal node until an R-NN pair is found. Each point is initially a singleton cluster and a terminal node of a single-node chain. The sequential algorithm picks an arbitrary node to start growing from. After an R-NN pair is found, the R-NN pair is then merged, and the chain is grown again to find another R-NN pair to merge. After −1 merges, the algorithm finishes, producing a hierarchy of clusters. Example. We now give an example of the definitions above by briefly describing running our ParChain framework for parallel HAC on the small data set in Figure 1. This example uses the complete linkage metric, where the distance between two clusters is the distance of the farthest pair of points, one from each cluster.
Our framework is based on the key insight that all R-NN pairs can be merged simultaneously, which provides parallelism. On each round, it merges all R-NN pairs in parallel (breaking ties lexicographically 1 to prevent cycles). Before the first round, each point in { , . . . , } is represented by a chain with only one node, and all points are singleton terminal nodes. The R-NN pairs are found by finding the nearest neighbors of all terminal nodes, which by definition are the clusters whose nearest neighbors are unknown at the beginning of a round. On the first round, we find the nearest neighbors for all points in parallel. Now  , and so we merge it and create a dendrogram node for cluster { , , }. We do not need to find the nearest neighbor of in this round, because it is not a terminal node and we know that its nearest neighbor will not change (due to the reducibility property which will be defined more formally in Section 2). On the third round, { , , }, { , }, and are terminal nodes. We find the nearest neighbors for them and merge the R-NN  Figure 1: (a) The four rounds of the nearest-neighbor chain algorithm on a six point data set using the complete linkage metric. The black circles are clusters containing more than one point. An arrow from point to point means that 's cluster's nearest neighbor is 's cluster. The orange arrows are new neighbors found on this round; the green arrow means the nearest neighbor is not updated on this round due to the reducibility property. On Round 2, we avoided a nearest neighbor search for cluster { }. The numbers on the arrows between pairs of points are the distances between the clusters that the points belong to according to the complete linkage metric (furthest distance between a pair of points, one from each cluster). Ties are broken lexicographically. (b) The dendrogram for complete linkage clustering. The label on each internal node corresponds to the furthest point pair in the two clusters that are merged in the algorithm, and its distance is equal to the node's height in the dendrogram. ParChain achieves high space efficiency and parallelism, which enables it to scale HAC to large data sets that are orders of magnitude larger than those used in previous work. There are two challenges in achieving both space efficiency and high parallelism. The first challenge is to maintain all chains and merge reciprocal nearest neighbor clusters correctly and efficiently in parallel. Unlike Jeon and Yoon's algorithm [33], which is based on locks (and has limited parallelism for large core counts), we use lock-free approaches based on filtering and atomic operations (Section 3). The second challenge is to efficiently find the nearest neighbors of clusters when growing the chain, without storing the distance matrix. We introduce a range query optimization that significantly reduces the number of distance computations used to find the nearest neighbor of a cluster for low-dimensional data sets in Euclidean space (Section 4), as well as a new caching technique that stores a subset of previous distance computations that are likely to be reused to further accelerate nearest neighbor searches (Section 5). In the example in Figure 1, the range query optimization avoids computing the distance between clusters { } and { , } in round 2 when { } searches for its nearest neighbor, because only clusters { , } and { } will be within the range. The caching technique avoids storing all pairs of distances among the six points. In contrast, many previous methods [23,31,33,51,52,58] require a quadratic-space distance matrix and compute distances to all other clusters when searching for the nearest neighbor of a cluster.
We apply ParChain to develop new linear-space parallel HAC algorithms for the complete, Ward's, and average linkage criteria. Our framework can be applied for any linkage criteria that satisfies the reducibility property, which ensures that the nearest neighbor distance of clusters can never be smaller as clusters merge (defined more formally in Section 2).
Though the worst case time complexity of our algorithms is ( 3 ), we observe that the running time is close to quadratic in practice on low-dimensional data sets because the range query is able to filter out many clusters. Many spatial, sensor, and computer vision data sets, where HAC is applicable, are low dimensional. In Section 6, we show experimentally on a variety of real-world and synthetic data sets (up to 16 dimensions) that our algorithms achieve 13.75-54.23x self-relative speedup on a 48-core machine with two-way hyper-threading. We also achieve 5.8-110.1x speedup over the state-of-the-art parallel implementations. Our algorithms use up to 237.3x less space than existing implementations, and are able to scale to larger data sets with tens of millions of points, which existing algorithms are not able to handle.
We summarize our contributions below:

BACKGROUND
The input to the hierarchical agglomerative clustering (HAC) problem is a data set to be clustered and a linkage criteria that specifies how distances between clusters are computed. The output of HAC is a tree called a dendrogram, where the height of each dendrogram node represents the dissimilarity between the merged two clusters according to the desired linkage criteria. A flat clustering, which assigns the same ID to every object in the same cluster and different IDs to objects in different clusters, can be obtained by cutting the dendrogram at some height. Thus, cutting the dendrogram at different heights gives clusterings at different scales. An example of a dendrogram is shown in Figure 1(b). In the rest of the section, we present our notations, the three linkage criteria considered in this paper, and some relevant techniques used by our algorithm. Notation. Let be a length-vector in -dimensional space, and let ∥ ∥ denote the 2 norm of , i.e., ∥ ∥ =

Linkage Criteria
We now formally define the linkage criteria considered in this paper. We use the Euclidean distance metric for all linkage criteria. For Table 1: Definitions, work, radius value, and optimizations used in our range query for different linkage criteria.

Relevant Techniques
d-tree. A d-tree is a binary spatial tree where each internal node contains a splitting hyperplane that partitions the points contained in the node between its two children. The root node contains all of the points, and the d-tree is constructed by recursing on each of its two children after splitting, until a leaf node is reached. A leaf node contains at most points for a predetermined constant . The d-tree can be constructed in parallel by performing the split and constructing each child in parallel. The bounding box of a node is the smallest rectangular box that encloses all of its points. Nearest-Neighbor Query. A nearest-neighbor query takes a set of points P and a query point , and returns for its nearest neighbor in P (besides itself if ∈ P). An all-nearest-neighbor query takes a set of points P, and returns for all points in P its nearest neighbor in P besides itself. The all-nearest-neighbor query can be performed efficiently using a dual-tree traversal [14,15,47], which we have parallelized. Range Query. A range query takes a set of points P, constructs a data structure to store the points, and reports or counts all points in some range . In this paper, we use balls to represent the ranges, and we use d-trees to store the points. Other Parallel Primitives. A parallel filter takes an array and a predicate function , and returns a new array containing ∈ for which ( ) is true, in the same order that they appear in . A parallel reduce takes as input a sequence [ 1 , 2 , . . . , ] and an associative binary operator ⊕, and returns the overall sum (using ⊕) of the elements ( 1 ⊕ 2 ⊕ . . . ⊕ ). A parallel hash table stores key-value pairs, and supports concurrent insertions, updates, and finds. WriteMin is a priority concurrent write that takes as input two arguments, where the first argument is the location to write to and the second argument is the value to write; on concurrent writes, the smallest value is written [61]. WriteMax is similar but writes the largest value.

THE PARCHAIN FRAMEWORK
In this section, we present our framework ParChain for parallelizing the nearest-neighbor chain (NNC) algorithm, which works for all linkage criteria that satisfy the reducibility property explained in Section 2.1. The NNC algorithm exposes more parallelism than the naive generic algorithm, where only the R-NN pair with minimum distance is merged, by allowing multiple R-NN pairs to be merged simultaneously. Hence, our framework grows multiple chains and merges all R-NN pairs simultaneously in parallel.
Jeon and Yoon's algorithm [33] uses a similar approach for to grow multiple chains in parallel, but it does merges R-NN pairs asynchronously. It designates some threads for updating chains, and other threads for updating the cluster distances. Their algorithm also uses locks and requires quadratic memory for maintaining the distance matrix. In contrast, our algorithm proceeds in rounds where on each round, all chains are grown and all R-NN pairs are merged. Our algorithm is lock-free, and uses linear space as we avoid using the distance matrix. In addition, Jeon and Yoon's algorithm searches for the nearest cluster naively by computing the distances to all other clusters, whereas we have optimizations for finding the nearest clusters when growing the chain, which will be discussed in Sections 4 and 5.

ParChain Framework
We now formally describe the ParChain framework (Algorithm 1). ParChain gives rise to fast and space-efficient HAC algorithms. The main idea of ParChain is to avoid storing most cluster distances, and compute them on the fly using an optimized range search that considers only a small number of neighboring clusters. We also cache some of the cluster distances to reduce computational cost.
The input to the algorithm is a set of points P, a structure that is used to compute the distances between clusters based on the linkage criteria, and an integer ≥ 0 for the cache size. We store (cache) only ( ) cluster distances for an integer ≥ 0 chosen at run time. The highlighted parts of Algorithm 1 (Lines 4 and [16][17][18] are only required for > 0, and we will discuss them in Section 5. is able to compute the distance between two clusters, and may maintain some extra data to accelerate distance computations, such as the means and variances of clusters. It also specifies the Lance-Williams formula if > 0, which will be used for updating the entries of cached distances between clusters. Initialization. On Lines 1-6, the algorithm initializes the required data structures. It first creates dendrogram nodes (Line 1) and L, a set of chain nodes (Line 2). These nodes are used for the singleton clusters at the beginning. We create a set of active clusters A, initialized to contain all of the singleton clusters (Line 3). We also create a parallel hash table for each cluster to cache cluster distances if > 0 (Line 4). In each chain node, we store its successor (succ), its predecessor (pred), and the distance to its predecessor (pred. ) if there is one. All chain nodes initially do not have any successor or predecessor. Z represents the set of terminal nodes at the beginning of the round, and is initialized to contain the singleton chains (Line 5). The algorithm also initializes a d-tree on the points . The d-tree (Line 6) is used to accelerate nearest cluster searches. Chain Growing and Merging. On Lines 8-15, in parallel we grow all of the chains using the information in Z. We merge a node in Z with its nearest neighboring cluster if they form an R-NN pair. Specifically, on Line 8, to grow the chains we find the nearest neighbors of all current terminal nodes in Z using a d-tree range search optimization, which will be described in Section 4. ParChain can quickly compute the distances of a cluster to other clusters by considering only a small number of candidates, without needing to maintain a distance matrix. The nearest neighbors and the distances are stored in E. On Lines 10-12, we update the successor and predecessor of each terminal node in parallel to maintain the chains. If a terminal node is the nearest neighbor of multiple clusters, the WriteMin ensures that its predecessor is the cluster closest to it. Then on Line 13, we find all R-NN pairs using a parallel filter by checking for each terminal node if its successor has a successor that is itself. All R-NN pairs are stored in an array M with their distances. On Line 15, we create a new dendrogram node , ,new to represent the merged cluster for each R-NN pair in M, which will have and as children, and store the distance between the merged clusters.
After the merges, we need to update the other data structures to prepare for next round. On Lines 16-18, we update the cache tables with new distances Section 5. On Line 19, we update the extra data structures, such as the d-tree and clusters' mean and variance. On Line 20, we update the set of active clusters by including active clusters not in M (not merged this round) and the newly merged clusters from this round. Finally, on Line 21, we obtain the new set of terminal nodes Z using a parallel filter on the active clusters A. Work and Space Complexity. We summarize the time and space complexity of the state-of-the-art algorithms in Table 3. ParChain is the only algorithm that requires linear memory and works for Table 3: Worst-case work and space bounds of state-of-the-art HAC algorithms. * The authors of [43,57,64] do not report the work complexity.

Algorithm
Work Space Restrictions ParChain ( 3 ) ( ) Reducibility NN-Chain [33,52,66] ( 2 ) ( 2 ) Reducibility Generic [31,52,58,59] ( 2 log ) ( 2 ) Lance-Williams Althaus et al. [3] ( Disk-based a broad set of linkage criteria (any that satisfy the reducibility property). The main computational cost in our framework is in finding the nearest neighbors of the terminal nodes on each round, and updating the cache tables and other data structures maintained by the distance structure . Sections 4 and 5 present our novel approaches for efficiently computing nearest neighbors efficiently with low space. We now analyze the work of our framework. Let Z and A be the sets Z and A at the beginning of round , respectively. The initialization (Lines 1-6) take ( log ) work, dominated by the dtree construction work on Line 6. Lines 10-18 and Lines 20-21 take ( |Z |) work across all rounds, plus the cost of all cluster distance computations, denoted by . Line 19 takes ( |A | log |A |) work because we need to re-construct the d-tree of cluster centroids in this step. Finally, Line 8 takes ( |A ||Z |) work because for each terminal node, we need a range query on the d-tree of cluster centroids. Thus, we have that the work of ParChain is ( + ), where = |A |(|Z | + log |A |). In the worst case, the work is ( 3 ), but we show in Section 6 that in practice both and are close to quadratic and ParChain is orders of magnitude faster than the ( 2 ) work algorithms [52], even using a single thread. We expect our algorithm to give improvements on most low-dimensional data sets.
The space usage of our framework is ( (1 + )) because all data structures except the caches require linear memory, and the caches require ( ) memory. The d-tree requires memory linear in the number of points in the tree.

NEAREST-NEIGHBOR FINDING
We will now describe how to efficiently perform nearest-neighbor finding (Line 8 of Algorithm 1) for the three linkage criteria: complete, Ward's, and average linkage. We assume that we compute distances between clusters on the fly. We describe an optimization in Section 5 that uses cache tables to store some of the distances.
While a standard nearest-neighbor search is done on points, we are searching for nearest neighbors of clusters with distances based on the linkage criteria. Our d-trees store centroids of clusters of points, which we use to find nearby clusters to our query cluster. We then perform exact distance computations from our query cluster to these clusters. Unlike in standard nearest-neighbor searches, it is harder to prune in our case as the distances between clusters centroids do not necessarily correspond to distances between clusters. Instead, we compute a different search area for each cluster based on an upper bound on the distance between the query cluster and its nearest neighbor. This upper bound can be a distance between the query cluster and any other cluster. We provide a novel heuristic for finding a good upper bound on the distance to the nearest cluster, and only search within this distance in Sections 4.1-4.3. In ParChain-no-range ParChain Figure 2: Runtimes on 48 cores with two-way hyper-threading of using our optimized range query compared to not using the range query and computing the distances to all other clusters on the fly to find the nearest neighbor.

Algorithm 2: Find Nearest Neighbor
Input: d-tree T, distance structure , chains L, and set of terminal nodes Z Output: nearest neighbors of nodes in Z 1 Initialize E with a (Null, ∞) entry for each terminal node.  Figure 2, we present the performance of using our optimized range query compared to the naive method of computing the distances to all other clusters to find the nearest neighbor. We see that our optimized range query gives a 7.8-1892.4x speedup on the two example data sets. Algorithm. Given the d-tree T built on the centroid of clusters, distance structure , chain nodes L, and set of terminal nodes Z, Algorithm 2 finds the nearest cluster and the distance to it for each terminal node's cluster and stores them in E.
When finding the nearest neighbor of cluster , we search all points within some ball (¯, ) obtained from ( , ), which is a ball centered at centroid¯with radius . The radius depends on the linkage method of and an distance between and another cluster (Lines 2-9 of Algorithm 2). Since we use centroid distances, rebuilds T to be a d-tree of only the centroids of current clusters at the end of each round of Algorithm 1 (Line 19 of Algorithm 1). If has a predecessor, we set the distance to be the distance between and its predecessor (Lines 4-5 of Algorithm 2). Otherwise, we find the distance to another cluster for computing the radius for the search. Specifically, we use the distance to the cluster whose centroid is the closest to the current cluster, which can be computed using a parallelized nearest-neighbor query on the d-tree of centroids [9] (Lines 6-7 of Algorithm 2).
For the range query on Line 9 of Algorithm 2, we use the parallel range query in Algorithm 3. Given a query cluster , a d-tree node , a ball representing the range, a linkage function , and a set E of pairs of nearest neighbor candidates and distances of terminal nodes, the algorithm processes all of 's points that are in the ball to update the nearest neighbor candidates in E. Since we only process the points in the ball, on Line 1, the range query terminates if the bounding box of the tree node does not overlap with the query range. Otherwise, the range query will either process all of the points both in the node and in the ball using the update_nearest_neighbor subroutine (Algorithm 4) if it is a leaf node (Lines 2-3 of Algorithm 3) or recurse on its two children in parallel (Lines 4-6 of Algorithm 3).
In each update_nearest_neighbor( , ) call, we check if some cluster is closer to than its current nearest neighbor candidate, and if so we update 's nearest neighbor in E with Algorithm 4. We also update 's nearest neighbor to be if is closer to than its current nearest neighbor candidate. In Algorithm 4, if = 0, we will compute the distance between and on the fly (Lines 7, 9, and 10). If > 0, we will first check the cache and use a cached distance if possible (we describe more details in Section 5).
As an optimization for the first round, we know that the distances between clusters is exactly the same as the distances between their centroids in the first round, and thus we can efficiently prune searches in the d-trees. Therefore, we use an all-nearest-neighbor query for the first round, which we implemented by parallelizing the dual-tree traversal algorithm by March et al. [47]. At a high level, our algorithm processes recursive calls of the dual-tree traversal in parallel and uses WriteMin to update the nearest neighbors of points. A dual-tree traversal allows more pruning than when running individual nearest neighbor queries for each point.
In the rest of the section, we will describe the radius of the search ball for each linkage method. We will show that a cluster's nearest neighbor must have its centroid inside the ball.

Ward's Linkage
For the range query, we can use a ball with radius = √︃ | |+ min 2 min | | , where is the distance between and some cluster and min is the size of the smallest current cluster. We can obtain min using a parallel reduce on the sizes of all clusters. Figure 3(a) illustrates the range search for Ward's linkage.
Since = Δ Ward ( , ), any cluster that is closer to than must have ∥¯−¯∥ 2 ≤ 2 | |+ | | 2 | | | | . The right-hand side of the inequality becomes smaller for larger | |, thus we can upper bound the distance between 's centroid and 's centroid (i.e., The green arrows specify the predecessor of , and the distance between them is . The red boxes are the bounding boxes of a d-tree node and its four descendants. In (a), Δ( , ) is computed because¯is in the ball; however, Δ( , ) is not computed because¯is outside the ball. min is the size of the smallest current cluster. In (b), Δ( , ) is computed because is completely in the ball. Δ( , ) will not be computed, because the three shaded bounding boxes of d-tree nodes do not intersect with the ball, and so some of 's points will not be included in the count. We only compute the distance to a cluster if all of its points are included in the count.

∥¯−¯∥) by = √︃
| |+ min 2 min | | . Therefore, no cluster outside of the ball centered at¯with radius can be closer to than , and thus we only need to search for 's nearest neighbor inside this ball.

Average Linkage
In average linkage, the distance between two clusters is the average distance between all pairs of points, one from each cluster. For the range query, we use a ball centered at¯with radius = and = √︁ for Euclidean distance metric and squared Euclidean distance metric, respectively. As before, is the distance between and some cluster . Now, we show that the nearest neighbor of must have its centroid inside the ball. First, we show this for the Euclidean distance metric. Observe that if we have ∥¯− ∥ ≤ Δ avg-1 ( , ), then we have ∥¯−¯∥ ≤ Δ avg-1 ( , ) ≤ Δ avg-1 ( , ) = = , and thus the centroid¯is within the ball.

Complete Linkage
In complete linkage, the distance between two clusters is the maximum distance between a pair of points, one from each cluster. For the range query, we use a ball with radius = centered at centroid¯, where is the distance between and some cluster. By definition of the complete linkage function, the cluster distance must be no smaller than distance between their centroids, and so the nearest neighbor of has its centroid within the search ball. Range Query Optimization. For complete linkage, we can reduce the number of cluster distance computations by only computing the distance to a cluster if it is completely within the search ball. With this observation, we can optimize the algorithm by keeping the d-tree to be T , the d-tree of all points, and avoiding updating it to be the d-tree of centroids on every round. Figure 3(b) illustrates the optimized range search for complete linkage. We will prove the correctness of this optimization at the end of the subsection.
Since now T is always T , we need to slightly modify Algorithm 2 and Algorithm 3. On Line 7 of Algorithm 2, we search for the point ∉ closest to¯in T , and let be the distance between and the cluster of this point. We can use a parallel union-find structure [22] to ignore points in . For Algorithm 3, the range query might be able to terminate before Line 2 if the tree node satisfies some conditions. For each range search, we keep a count that eventually upper bounds the number of points within the ball for each cluster. In each for-loop on Line 3 of Algorithm 3, we now loop over points instead of centroids, and we atomically increment the count for 's cluster by 1 because this means we have found one more point in this cluster that is within the ball. The cluster IDs can also be maintained and queried using the parallel union-find data structure. Right before Line 2 of Algorithm 3, if all points in the d-tree node are from the same cluster , we atomically increment the count of cluster by the size of the node and prune the search; otherwise, we continue the search and recurse on the children. This gives an upper bound on the number of points in the cluster within the ball, because the ball lies inside the d-tree bounding boxes traversed.
We preprocess the tree such that in the range search we can determine in constant time if all points in the node are from the same cluster, and if so which cluster it is. Specifically, we mark the d-tree nodes with a cluster ID if all points in the node are from the same cluster, or with NULL if the points in the node belong to multiple clusters. This can be computed by recursively checking the ID of the two children of a node starting from the root, and storing the cluster ID of the children if all of their points are from the same cluster. We update this information on every round.
After processing a point or a node, if we incremented the count of a cluster , we check if the count of is equal to the size of . If so, this means that all of 's points may be within distance = . In this case, we compute the distance between the and this cluster, and use a WriteMin to update the nearest neighbor of in E (Lines 7, 9, and 10 of Algorithm 4). Finally, we show below that 's nearest neighbor must be a cluster completely within search area by claiming that clusters with Euclidean distance metric using 48-cores with two-way hyper-threading for varying cache sizes (values of ). The data sets are labeled on the -axis and are described in Section 6 (caption of Table 4).
points outside the ball must have a distance larger than to . Since is the distance between and some cluster, must have a distance no larger than to . Suppose the distance of the furthest point pair between and is Δ comp ( , ) = ( , ). Since the average Euclidean distance between points in two clusters is not smaller than the distance between their centroids (shown in Section 4.2), applying this property to and { ′ } for any point ′ ∈ , we see there must exists some ′ ∈ such that (¯, ′ ) ≤ ( ′ , ′ ). Since ( , ) is the furthest point pair, we have that (¯, ′ ) ≤ ( ′ , ′ ) ≤ ( , ). Thus, if Δ comp ( , ) = ( , ) ≤ , then all points in must be within comp (¯, ). As a result, we only need to consider a cluster as the nearest neighbor candidate of and actually compute the distance to it if all of its points are inside the ball. Dual-Tree Traversal. When computing cluster distances (Line 7 of Algorithm 4) for complete linkage, we use our parallel dual-tree traversal algorithm described earlier in the section. We need to find the distance of the farthest pair of points, and so we use WriteMax instead of WriteMin for storing the farthest distance seen. In order to perform the dual-tree traversals, creates a d-tree for each cluster at the end of each round (Line 19 in Algorithm 1).

CACHING INTER-CLUSTER DISTANCES
For some linkage function and metric combinations, such as average linkage with the Euclidean distance metric, computing inter-cluster distances can be expensive. We can avoid some recomputations of cluster distances by caching some previously computed distances for each cluster using a cache table H , represented using a parallel hash table. Users can specify a constant size of each cache based on the available memory. The total memory usage is ( (1 + )), which is less than the quadratic memory required by the distance matrix approaches. Sometimes, a larger table will lead to faster computations because we can cache more distances and avoid more recomputations. Due to the optimizations in Section 4, the distances that we compute will tend to be close to , and hence stored in H . These distances are more likely to be reused in future nearest neighbor queries.
We present a comparison of running times of average linkage with the Euclidean distance metric on several data sets using different cache sizes in Figure 4. We see that using caching improves the running times by up to a factor of 8.98x compared to not using caching. We found similar trends on other data sets. We will discuss more about our implementation's memory usage in Section 6.
In the rest of the section, we assume > 0 and describe how to query cluster distances from the cache tables, insert new entries after computing cluster distances during nearest neighbor queries, and update the tables after merging clusters. Querying and Inserting Distances between Clusters. The cache tables can be used to reduce cluster distance computations because we can insert the computed distances to the tables and query for them if we want to use them again. Now we describe how the cache is used to update the nearest neighbor candidate in the nearest neighbor search (Algorithm 4). With > 0, we might have already cached the distance Δ( , ) in one or both of the tables H and H when we find in 's range. Therefore, we first query for the distance in the cache tables (Line 2), and only compute the distance if the return value is NOT_FOUND; otherwise we can directly use the queried distance to update the nearest neighbor candidate in E (Lines 3-6). If we compute the distance (Line 7), we will attempt to insert it into both of the tables (Line 8). The insertion may fail for a cache table if it is full, i.e., it already contains entries. Since we insert distances between and the clusters that are within its search range in all rounds so far, the distances stored in H are likely to be between and nearby clusters. Thus in later rounds, these cached distances are more likely to be queried. On Lines 9-10, WriteMin updates the nearest neighbors of and in E. When querying Δ( , ) with get_cached_dist( , ) (Line 2), we search for the entry with key in H , and the entry with key in H . If in a cache table, the key does not exist, then the query fails. If the queries in both tables fail, we return NOT_FOUND. If the search is successful in one of the tables, we return the distance stored in the table. We search in both cache tables since the caches are of limited size, and so the distance could potentially be stored in just one of the two tables. Updating Cache Tables after Merging Clusters. We now describe how to update the entries in the cache tables after clusters are merged (Lines 16-18 of Algorithm 1). If during a round and are merged into a new cluster , we will try to compute the distance between and all clusters whose subclusters' distance(s) with or are stored in H ∪ H . These distances can be used to accelerate the computation of Δ( , ) using the Lance-Williams formula described in Section 2.1.
The update_cached_dists function called on Line 18 of Algorithm 1 is presented in Algorithm 5. On Line 2, we loop over the distances ′ in the cache tables of H and H . Without loss of generality, assume ′ = Δ( , ℓ ) is a distance in H between and some cluster ℓ (the case for an entry in H is similar). Line 3 skips over the entries that represent distances between and , since they are now merged. Otherwise, there are two cases. In case (1), ℓ is also a cluster merged in this round (Lines 4-6), and we let be the cluster that ℓ merged into. We compute the new distance  (2), ℓ is not a new cluster merged in this round (Lines 7-9), and we have = ℓ . We can also use ′ to accelerate the computation of Δ( , ). Figure 5 illustrates one loop of the algorithm where the entry being considered is ′ = Δ( , ℓ ) ∈ H (shown in the bottom gray box). In both cases, we store entry = Δ( , ) computed from ′ into both H and H on Line 10. Now we describe the update rule for computing . For case (2), we can just directly apply the Lance-Williams formula [38] introduced in Section 2 and compute Δ( ∪ , ) from Δ( , ), Δ( , ), and Δ( , ). For case (1), we can apply the Lance-Williams formula and compute Δ( where ℓ ′ is the cluster that ℓ merges with to form . To compute Δ( , ℓ ) and Δ( , ℓ ′ ), which are not cached since is a newly merged cluster, we can apply the Lance-Williams formula again since = ∪ . Below we give the update rule for average linkage with Euclidean distance metric as an example. For case (2), where = ℓ , we have For case (1), let ℓ and ℓ ′ be 's children. We have If the distances between ℓ ′ ( ) and one or both of and are also cached in case (1) (case (2)), we can also query them and accelerate the computation of by avoiding some distance recomputation. For example, in update rule (1) is also cached, we can query for it and compute in constant time; otherwise, we need | || | point distance computations to find , which is less than the | || | point distance computation required by a brute force method.
During the nearest neighbor range search, two clusters and might find each other as nearest neighbor candidates, and both want to compute Δ( , ) in parallel. A similar situation can happen when updating the cache entries for and . Our implementation avoids these duplicate distance computations by having each cluster first insert a special entry into the hash table H min( , ) (or H max( , ) if H min( , ) is full), and then only compute the distance if the insertion was successful. The special entry can only be successfully inserted once for each pair of clusters and , and so Δ( , ) will only be computed once.

EXPERIMENTS
Testing Environment. We perform experiments on a c5.24xlarge machine on Amazon EC2, with 2 Intel Xeon Platinum 8275CL (3.00GHz) CPUs for a total of 48 hyper-threaded cores, and 192 GB of RAM. By default, we use all cores with hyper-threading. We use the g++ compiler (version 7.5) with the -O3 flag, and use Cilk Plus, which is supported in g++, for parallelism in our code [40].
For parallel experiments, we use numactl -i all to balance the memory allocation across nodes. We also perform three runs of each parallel experiment and report the smallest running time. We allocate a maximum of 15 hours for each run of a running time test, and do not report the times for tests that exceed this limit.
We test the following implementations for HAC. We refer to complete linkage as comp, Ward's linkage as Ward, average linkage with Euclidean distance metric as avg-1, and average linkage with squared Euclidean distance metric as avg-2.
• PC Our parallel ParChain framework is implemented in C++, using the range query and caching optimizations. • PC-mr A parallel C++ NNC implementation that uses a distance matrix and merges all R-NNs in each round. All cluster distances are obtained from the distance matrix. It uses the range query optimization to find the nearest neighbors. • PC-m A naive parallel C++ NNC implementation that uses a distance matrix and merges all R-NNs in each round. All cluster distances are obtained from the distance matrix. A parallel reduce is used over the distance matrix rows to find the nearest neighbor instead of using a range query. • scipy (sc) [51] Scipy's serial implementation in Cython, which uses the NNC algorithm with a distance matrix for all of the linkage criteria tested. • sklearn (sk) [59] Scikit-learn's serial implementation in Cython, which uses a distance matrix. It has a heap for all distances and merges the global closest neighbor pair on each round. • fastcluster [52] A serial C++ implementation of HAC with a Python interface. It contains two approaches of implementations of HAC-one generic implementation (fc-gen) uses the naive algorithm where the global R-NN pair is merged in each round, and the other (fc-nnc) is based on the NNC algorithm using a distance matrix. For Ward's linkage, fastcluster has a linear space implementation for the naive method that work by computing the cluster distances on the fly using the cluster centroids. We also wrote linear-space implementations for the NNC algorithm that compute distances on the fly for Ward's linkage and average linkage with squared Euclidean distance metric. We report the running time of the linear space implementation when available. • fastprotein (fp) [31] A parallel C++ implementation that parallelizes the naive NNC algorithm by computing the global R-NN and updating the distance matrix in parallel on each of the − 1 rounds. It only supports complete and average linkage.
• Jeon (Je) [33] A parallel C++ implementation of the parallel NNC algorithm by Jeon and Yoon. [33], which only supports average linkage with the Euclidean distance metric. • Althaus (Al) [3] Our parallel C++ implementation of Althaus et al. 's complete linkage algorithm that uses linear memory. Data Sets. We use both synthetic and real-world datasets, all of which fit in the RAM of our machine. Let be the number of points. The GaussianDisc data set contains points inside a bounding hypergrid with side length 5

√
. 90% of the points are equally divided among five clusters, each with a Gaussian distribution. Each cluster has its mean randomly sampled from the hypergrid, a standard deviation of 1/6, and a diameter of √ . The remaining points are randomly distributed. The UniformFill data set contains points distributed uniformly at random inside a bounding hypergrid with side length √ . We generate the synthetic data sets with 10 million points for dimensions = 2 and = 5.
We also use two existing simulation datasets. UCI1 [7] is a 10dimensional data set with 19, 020 data points. This data set is generated to simulate registration of high energy gamma particles [29]. UCI4 [35] is a 10-dimensional data set with 100, 000 data points. This data set has a pseudo-periodic time series for each of its dimension, and hence is likely to form long chains.
We use the following real-world data sets. GeoLife [71] is a 3dimensional data set with 24, 876, 978 data points. This data set contains user location data, and is extremely skewed. HT [30] is a 10-dimensional data set with 928, 991 data points containing home sensor data. CHEM [25] is a 16-dimensional data set with 4, 208, 261 data points containing chemical sensor data.
When referring to the data sets in this section, we use a prefix to indicate its dimensionality and suffix to indicate its size. To obtain smaller data sets, we randomly sample from the corresponding larger data sets. The letter "U" indicates UniformFill and "G"indicates "GaussianDisc". Cache Sizes. In our experiments, we use a cache size of = 64 for avg-1 and = 0 for complete, Ward, and avg-2 except otherwise noted. The choice of = 64 will be discussed in more detail in Section 6.3 along with the benefit of our range query and caching optimizations. We use = 0 for complete, Ward, and avg-2 to show our framework achieves good performance on cheaper linkage criteria even without caching. Figure 6 shows the running times vs. number of threads for all of the serial and parallel implementations on three small data sets (2D-GaussianDisc-10K, 10D-UCI1-19K, and 10D-UCI4-100K). Implementations with a single data point are serial. We only compare them on the small data sets because the algorithms that require quadratic memory run out of memory for larger data sets.

Comparison with Other Implementations
We see that our implementation PC almost always outperforms existing implementations across all thread counts. Even the version of our algorithm using the distance matrix without the range query optimization (PC-m) is faster than all other implementations at higher thread counts. Unlike the existing parallel implementations, fastprotein (fp) and Jeon (Je), our implementations are more scalable since we merge all R-NN pairs on each round, and do not use locks. On the small data sets, using 48 cores with hyper-threading,  Table 4 and Figure 7 for running times and scalability of PC on larger data sets. Table 4: Runtimes (seconds) and self-relative parallel speedups of fastcluster and PC. "PC-1" is our runtime on 1 thread and "PC-48h" is our runtime using 48 cores with two-way hyper-threading. PC on avg-1 timed out on several larger data sets. For fastcluster on several larger data sets, we only have runtimes for Ward and avg-2 because it requires quadratic memory for comp and avg-1 and runs out of memory. For Ward, we report fc-gen because it is faster than fc-nnc. "-" means that the running time exceeds 15 hours. For the data set names, the first number indicates dimension, the letter "U" indicates UniformFill, "G" indicates "GaussianDisc", "1M" indicates 1 million data points, and "10M" indicates 10 million data points.  PC is 5.8-88.0x faster than fp and 37.5-110.1x faster than Je. Table 4 shows the running times for PC, fc-gen, and fc-nnc on larger data sets (Je does not scale to these data sets due to its quadratic memory requirement), and we see that PC is 64.77-733.90x faster than fastcluster on these data sets. On a single thread, we find that PC is 2.19-47.92x faster than the next fastest implementation (except on 10D-UCI1-19K for complete linkage, where PC is 1.58x slower than fc-nnc and fc-gen).
In Figure 6, PC shows limited scalability on higher thread counts, because these data sets are small and the overhead of using more threads is high relative to the work of the algorithm. However, in the next subsection, we show that PC is able to achieve higher parallel scalability on larger data sets.

Scalability
Scalability with Thread Count. Table 4 and Figure 7 present the runtime and scalability of PC on different numbers of threads for larger data sets, which most existing implementations do not scale to. For average linkage with the Euclidean distance metric, the speedups for several data sets are not shown since the singlethreaded experiments timed out. We see that using 48 Figure 7, we can see that on most data sets, our algorithm keeps scaling up until 48 threads. Scalability with Data Size. Figure 8 shows the runtimes of our algorithm PC on three data sets of varying sizes using 48 cores with two-way hyper-threading. We observe that PC scales well with data set size. The scalability is better for comp, Ward, and avg-2 than for avg-1 because avg-1 always requires quadratic work to compute cluster distances, while Ward and avg-2 require constant time to compute them and comp usually requires less than quadratic time to compute them due to pruning. Runtime Decomposition. We now describe the breakdown of running time across different steps of ParChain, as well as the scalability of each step. Figure 9 shows the speedups and running times of steps of avg-1 using 48 cores with hyper-threading on different data sets. From Algorithm 1, "init" corresponds to Lines 1-6 where we initialize the data structures; "nn" corresponds to Lines 8-12 where we find the nearest neighbors of all terminal nodes and update the chains; "merge" corresponds to Lines 13-18 where we merge the R-NNs and update the cache tables if > 0; and "update" corresponds to Lines 19-21 where we update the data structures to prepare for next round.
From Figure 9 (top), we see that the "nn" and "merge" steps are more scalable with respect to thread count than the "update" and "init" steps. The reason is that in the "nn" step, we find the nearest neighbor of all terminal nodes in parallel, and in the "merge" step, we merge all R-NNs in parallel. The numbers of terminal nodes and R-NNs are usually much larger than the number of available threads, and thus there is a lot of opportunity for parallelism. The "init" and "update" step are less scalable because they have less work to be divided across threads.
From Figure 9 (bottom), we see that the "nn" and "merge" steps are less scalable with respect to data size than the "update" and "init" steps. This is because the "nn" and "merge" steps asymptotically dominate the work of the whole algorithm.

Analysis of Our Framework
We now discuss the effects of our range query and caching optimizations and show our running time is close to quadratic in practice. From Figure 6, PC-mr is 1.67x faster on average than PC-m on 48 cores with two-way hyper-threading, which shows that the benefit of using our optimized range query is larger than its overhead even on these small datasets. PC is 15.06x faster on average than PC-mr because although PC needs to compute some distances on the fly, PC avoids the overhead of computing the distance matrix avg-1 using = 256 for all clusters (we used = 128 for GeoLife due memory limitations). We use a larger cache size than our previous experiments so that fewer clusters hit the size limit, which gives us a more accurate analysis of cache usage. We found that the maximum average cache usage ranges from 6.7-84.5 slots. This explains why runtimes stop decreasing for cache sizes larger than 64 in Figure 4. If the cache is too large and many of the entries are empty, the runtime could be slower than using a smaller cache because when we merge the caches, we need to filter out the non-empty entries, and thus larger caches incur more overhead. We use = 64 for avg-1 in all other experiments to show that our framework can achieve good performance with a relatively small cache so that our memory overhead is minimal. Figure 4 shows that = 64 gives good performance across different data sets. Section 6.4 shows our memory usage is very small with = 64.
We now show that our running time is close to quadratic in practice. As described in Section 3, the work of our framework is bounded by = |A |(|Z | + log |A |) plus the cost of distance computations, . Figure 10 (left) shows that and are quadratic in the number of points in practice. Figure 10 (right) shows that only a very small fraction of clusters are included in the range queries and an even smaller fraction of distance computations between clusters are required in practice. Table 5 shows the memory usage (in megabytes) vs. data set size for 2D-GaussianDisc data sets for the different implementations. We measure the memory usage using the Valgrind Massif heap profiler tool [56]. For PC and PC-m, we use 48 cores with two-way hyperthreading. For fastprotein (fp) we use 36 threads, for Jeon (Je) we use 4 threads, and for Althaus (Al) we use 36 threads. The number of threads are chosen for best performance based on Figure 6.

Memory Usage
Our algorithm PC shows linear memory increase vs. data size on all methods while most other methods, except Al, fc-gen for Ward, and fc-nnc for Ward and avg-2 (which require linear space),   : Self-relative parallel speedup and runtimes (seconds) of steps of avg-1 using 48 cores with two-way hyper-threading on different thread counts and data set sizes. From Algorithm 1, "init" corresponds to Lines 1-6; "nn" corresponds to Lines 8-12; "merge" corresponds to Lines 13-18; and "update" corresponds to Lines 19-21. show quadratic memory increase vs. data size, which is consistent with the fact that they use a distance matrix. PC uses less memory than all other implementations, except that Al uses less memory for complete linkage, and Je uses less memory for the small 1K data set for avg-1. However, Al only works for complete linkage and is orders of magnitude slower than PC ( Figure 6). PC uses less than 2x of the memory used by Al. Overall, PC uses up to 237.3x less memory than existing implementations.

RELATED WORK
There is a rich literature in designing HAC algorithms. In the most naive algorithm, a distance matrix is used to maintain all pairwise distances between clusters. On each iteration, the matrix is searched to find the closest pair of clusters, which are then merged, and distances to this newly merged cluster are computed. The algorithm runs for − 1 iterations, after which a single cluster remains. A straightforward implementation of this algorithm gives ( 3 ) time, but it can be improved to ( 2 log ) time by storing matrix entries in heap-based priority queues [52,58]. The two popular Python libraries scipy [66] and scikit-learn [59] both provide sequential algorithms for HAC. The two libraries' implementations both compute and store a distance matrix. Fastcluster [52] contains three implementations of HAC-two heap-based naive algorithms for general linkage functions, where one uses the distance matrix and the other compute cluster distances on the fly, and an NNC algorithm that uses the distance matrix. Lopez-Sastre et al. [44] propose a sequential NNC algorithm that speeds up the chain construction using a dynamic slicing strategy that only searches for the nearest neighbor within some slices. Their algorithm only works for linkage functions where the distance can be expressed using centroids and variances.
There have been implementations that focus on reducing the in-memory space usage of HAC from quadratic to linear by writing the quadratic-space distance matrix to disk, and loading it into memory in smaller chunks [43,57]. These algorithms are sequential, and only merge one pair of clusters at a time. In contrast, our algorithm is parallel, and also does not require writing or loading additional information to and from disk. Moreover, the algorithms above are designed to take advantage of sparse distance metrics, where only some distances between data points are defined while other distances are considered to be "missing" and the points have "large" dissimilarity between them, making them less suitable for the Euclidean distance or squared Euclidean distance metrics.
There have also been many parallel algorithms developed for HAC, although it is difficult to parallelize in theory [27]. Olson [58] gives parallel algorithms, some of which parallelize the NNC algorithm by finding the nearest neighbor in parallel on each round, but still only merges one pair per round, and so there will always be − 1 rounds. Li [41] gives parallel HAC algorithms that store the distance matrix, based on an older theoretical model for a SIMD machine with distributed memory. Li and Fang [42] give parallel HAC algorithms on hypercube and butterfly network topologies. Du and Lin [23] give a parallel HAC algorithm on a cluster of compute nodes. Zhang et al. [70] propose a distributed algorithm for HAC that partitions the datasets using d-trees or quadtrees, and then for each leaf node, finds a region where the R-NN pairs might exist. In parallel, each compute node finds the local R-NN pairs in a region, and then global R-NN pairs are found from the local pairs. This method merges multiple R-NN pairs, but their paper does not specify how the distances between clusters are updated or computed after merges. Fastprotein [31] is a naive parallelization of fastcluster. Sun et al. [64] develop a parallel version of algorithms that write the distance matrix to disk and load chunks of it into memory [43,57]. However, they still only merge one pair of clusters at a time. Jeon and Yoon [33] present a parallel NNC algorithm using a distance matrix, which we discussed earlier. Althaus et al. [3] present a parallel complete linkage algorithm that uses linear main memory; however their algorithm requires − 1 rounds because they only merge the global R-NN pair on each round. Sumengen et al. [63] developed a distributed algorithm for HAC that recovers the exact HAC solution, and has good theoretical bounds under assumptions about the underlying cluster structure of the data. Similar to our algorithm, it merges reciprocal nearest-neighbor pairs in parallel. Their algorithm can be accelerated in practice by using a graph-building approach, and has good empirical performance on billion-scale datasets. In contrast, our work is focused on designing an optimized shared-memory framework for running HAC on the original points and the ( 2 ) distances between them, while accelerating the distance computation steps and distance storage steps. It would be interesting to compare both algorithms in a shared-memory setting in future work.
Besides the linkage criteria considered in this paper, other popular criteria for HAC include single, centroid, and median linkage. Single linkage with the Euclidean metric is closely related to the Euclidean minimum spanning tree problem, and can be solved efficiently using variants of minimum spanning tree algorithms [47,67]. Centroid and median linkage do not satisfy the reducibility property and cannot take advantage of the NNC algorithm. There has also been work on other hierarchical clustering methods, such as partitioning hierarchical clustering algorithms and algorithms that combine agglomerative and partitioning methods [10,18,41,48,60]. Finally, there has been work on analyzing the cost function of the HAC problem [13,16,17,21,50] and approximating the HAC problem on various linkage criteria and metrics [1,4,11,12,20,26,36,37,49].

CONCLUSION
In this paper, we presented ParChain, a framework that supports fast and space-efficient parallel HAC algorithms based on the nearestneighbor chain method. We introduced two key optimizations for efficiency, a range query optimization and a caching optimization. Using ParChain, we designed new parallel HAC algorithms for complete, average, and Ward linkage that outperform existing parallel implementations by 5.8-110.1x, while using up to 237.3x less space. It would be interesting future work to study how to improve the efficiency of ParChain by allowing approximation.