$\textrm{GPUSCAN}^{++}$: Efficient Structural Graph Clustering on GPUs

Let $G=(V,E)$ be an undirected and unweighted graph, where $V(G)$ is the set of vertices and $E(G)$ is the set of edges, respectively. We denote the number of vertices as $n$ and the number of edges as $m$, i.e., $n=|V(G)|$, $m=|E(G)|$. For a vertex $v\in V(G)$, we use ${\mathsf{nbr}}(v,G)$ to denote the set of neighbors of $v$. The degree of a vertex $v\in V(G)$, denoted by ${\mathsf{deg}}(v,G)$ is the number of neighbors of $v$, i.e., ${\mathsf{deg}}(v,G)=|{\mathsf{nbr}}(v,G)|$. For simplicity, we omit $G$ in the notations when the context is self-evident.

Note that the $\epsilon$-neighborhood of a given vertex $u$ contains $u$ itself as $\sigma(u,u)=1$. When the number of $\epsilon$-neighbors of a vertex is large enough, it becomes a core vertex:

Given a core vertex $u$, the structurally reachable vertex of $u$ is defined as:

Intuitively, a cluster is a set of vertices that can be structurally reachable from any core vertex. Formally,

Following the above definitions, it is clear that two vertices that are not adjacent in the graph have no effects on the clustering results even if they are similar. Therefore, we only need to focus on the vertex pair incident to an edge. In the following, for an edge $(u,v)$, we use the edge similarity of $(u,v)$ and the similarity of $u$ and $v$ interchangeably for brevity. We summarize the notations used in the paper in Table I.

For the thread safety of GPU programming, CUDA provides programmers with atomic operations, including $\mathsf{atomicAdd}$ and $\mathsf{atomicCAS}$. $\mathsf{atomicAdd(\emph{address},\emph{val})}$ reads the 16-bit, 32-bit, or 64-bit word old located at the address address in global or shared memory, computes $\mathsf{(\emph{old}+\emph{val})}$, and stores the result back to memory at the same address. These three operations are performed in one atomic transaction and returns old value. $\mathsf{atomicCAS(\emph{address},\emph{old},\emph{new})}$ (Compare And Swap) reads the 16-bit, 32-bit, or 64-bit word old located at the address address in global or shared memory, computes $\mathsf{(\emph{old}==\emph{old}?\emph{new}:\emph{old})}$, and stores the result back to memory at the same address. These three operations are performed in one atomic transaction and return old.

Structural graph clustering ($\mathsf{SCAN}$) is first proposed in [40]. After that, a considerable number of follow-up works are conducted in the literature. For the online algorithms regarding a fixed parameter $\epsilon$ and $\mu$, $\mathsf{SCAN}$++ [33], pSCAN [7] and $\mathsf{anySCAN}$ [43] speed up $\mathsf{SCAN}$ by pruning unnecessary structural similarity computation. $\mathsf{SparkSCAN}$ [44] conducts the structural clustering on the Spark system. $\mathsf{pmSCAN}$ [30] considers the case that the input graph is stored on the disk. $\mathsf{ppSCAN}$ [8] parallelizes pruning-based algorithms and uses vectorized instructions to further improve the performance.

For the index-based algorithms aiming at returning the clustering result fast for frequently issued queries, GS*-Index [39] designs a space-bounded index and proposes an efficient algorithm to answer a clustering query. [38] parallelizes GS*-Index and uses locality-sensitive hashing to design a novel approximate index construction algorithm. [26] studies the maintenance of computed structural similarity between vertices when the graph is updated and proposes an approximate algorithm that achieves $O(\log^{2}|V|+\log|V|\cdot\log\frac{M}{\delta^{*}})$ amortized cost for each update for a specified failure probability $\delta^{*}$ and every sequence of $M$ updates. [22] studies structural clustering on directed graphs and proposes an index-based approach to computing the corresponding clustering results. However, all these algorithms are CPU-based and cannot be translated into efficient solutions on GPUs due to the inherently unique characteristics of a GPU different from a CPU.

Besides $\mathsf{SCAN}$, other graph clustering methods are also proposed in the literature. However, all these methods are not based on structural clustering, and thus do not share the effectiveness of $\mathsf{SCAN}$ in detecting clusters, hubs, and outliers. [27, 1] provides a comprehensive survey on this topic.

For the GPU-based approach, $\mathsf{GPUSCAN}$ [37] is the state-of-the-art GPU-based $\mathsf{SCAN}$ solution, which conducts the clustering in three separate phases: (1) $\epsilon$-neighborhood identification. In Phase 1, $\mathsf{SCAN}$ computes the similarity of each edge and determines the core vertices in the graph. (2) cluster detection. After identifying the core vertices, $\mathsf{GPUSCAN}$ detects the clusters based on Definition 2.1 in Section 2. (3) hub and outlier classification. With the clustering results, the hub and outlier vertices are classified in Phase 3. The following example illustrates the three phases of $\mathsf{GPUSCAN}$.

Fig. 2 (a)-(b) show the $\epsilon$-neighborhood identification phase. $\mathsf{GPUSCAN}$ uses $E.u$ and $E.v$ arrays to represent the edges of the graph, a ${\mathsf{flag}}$ array in which each element stores whether an edge is similar or not, and $A.u$ and $A.v$ arrays to store the neighbors of each vertex. $\mathsf{GPUSCAN}$ uses a warp to computing the similarity of each edge $(u,v)$. According to Definition 2.1, the key to compute the similarity of $(u,v)$ is to find the common neighbors in $N[u]$ and $N[v]$. To achieve this goal, for each $w\in{\mathsf{nbr}}(u)$, $\mathsf{GPUSCAN}$ uses 8 threads in a warp to check whether $w$ is also in 8 continuous elements in ${\mathsf{nbr}}(v)$. As a warp usually contains 32 threads, four elements in ${\mathsf{nbr}}(u)$ can be compared concurrently. After the comparisons of these four elements are complete, if the largest vertex id of these four elements in ${\mathsf{nbr}}(u)$ is bigger than the largest one in the 8 continuous elements ${\mathsf{nbr}}(v)$, then, the next continuous 8 elements in ${\mathsf{nbr}}(v)$ are used to repeat the above comparisons. If the largest vertex id of these four elements in ${\mathsf{nbr}}(u)$ is smaller than the largest one in the 8 continuous elements ${\mathsf{nbr}}(v)$, then, the next 4 elements in ${\mathsf{nbr}}(u)$ are used to repeat the above comparisons. Otherwise, the next 4 elements in ${\mathsf{nbr}}(u)$ and the next 8 elements in ${\mathsf{nbr}}(v)$ are used to repeat the above comparisons. The common neighbors are found when all the elements in ${\mathsf{nbr}}(u)$ and ${\mathsf{nbr}}(v)$ have been explored. Take edge $(v_{0},v_{1})$ as an example. Fig. 2 (b) shows the procedure of $\mathsf{GPUSCAN}$ to compute their common neighbors. As ${\mathsf{nbr}}(v_{0})=\{v_{1},v_{2},v_{3},v_{4},v_{5},v_{6},v_{7}\}$, the first 4 elements $\{v_{1},v_{2},v_{3},v_{4}\}$ in ${\mathsf{nbr}}(v_{0})$ are compared first. For the 32 threads in a warp, threads $T_{0}-T_{7}$ handle the comparison for $v_{1}$. The remaining threads are assigned to handle the comparison for vertices $v_{2}$, $v_{3}$, and $v_{4}$ similarly. As ${\mathsf{nbr}}(v_{1})=\{v_{0},v_{2},v_{4},v_{7}\}$, threads $T_{0}-T_{7}$ finds $v_{1}\notin{\mathsf{nbr}}(v_{1})$. Similarly, threads $T_{8}-T_{15}$, $T_{16}-T_{23}$, $T_{24}-T_{31}$ find $v_{2},v_{4}\in{\mathsf{nbr}}(v_{1})$. After that, they find the largest vertex id $v_{4}$ of elements in ${\mathsf{nbr}}(v_{0})$ that were compared is smaller than the largest one $v_{7}$ in ${\mathsf{nbr}}(v_{1})$. Then, the remaining 3 elements $\{v_{5},v_{6},v_{7}\}$ in ${\mathsf{nbr}}(v_{0})$ are used to repeat the comparisons similarly. In this comparison, threads $T_{16}-T_{23}$ find that $v_{7}$ is the common element of ${\mathsf{nbr}}(v_{0})$ and ${\mathsf{nbr}}(v_{1})$. Therefore, $\sigma(v_{0},v_{1})=(2+3)/\sqrt{5*8}=0.79>0.6$, which means $v_{1}$ and $v_{0}$ are similar, and the corresponding element in the $\mathsf{flag}$ array is set as 1.

Fig. 2 (c) shows the graph clustering phase of $\mathsf{GPUSCAN}$. Following Definition 2.1, only similar edges could be in the final clusters. Therefore, $\mathsf{GPUSCAN}$ builds two new arrays $E^{\prime}.u$ and $E^{\prime}.v$ by retrieving the similar edges in $E.u$ and $E.v$. Consider the initialization step in Fig. 2 (c), as $v_{0}$ and $v_{1}$ are similar while $v_{0}$ and $v_{2}$ are not similar, only $v_{0}$ and $v_{1}$ are kept in $E^{\prime}.u$ and $E^{\prime}.v$. Then Based on $E^{\prime}.u$ and $E^{\prime}.v$, $N_{\epsilon}[u]$ can be obtained easily. For example, for $v_{0}$, $N_{\epsilon}[v_{0}]=\{v_{0},v_{1},v_{4},v_{7}\}$, $|N_{\epsilon}[v_{0}]|=4>3$, thus, $v_{0}$ is a core vertex. Similarly, $v_{4}$ and $v_{7}$ are also core vertices. Then, it removes the edges in $E^{\prime}.u$ and $E^{\prime}.v$ not incident to core vertices. Moreover, $\mathsf{GPUSCAN}$ represents all clusters by a spanning forest through a $\mathsf{parent}$ array. In the $\mathsf{parent}$ array, each element represents the parent of the corresponding vertex in the forest, and the corresponding element for the vertices in $E^{\prime}.u$ and $E^{\prime}.v$ is initialized as the vertex itself, which is also shown in Initialization of Fig. 2 (c).

After that, $\mathsf{GPUSCAN}$ detects the clusters in the graph iteratively. In each odd iteration, for each vertex $u$ in $E^{\prime}.u$, it finds the smallest neighbor $v$ of $u$ in $E^{\prime}.v$ and takes the smaller one between $u$ and $v$ as the parent of $u$ in the $\mathsf{parent}$ array. Then, it replaces each vertex in $E^{\prime}.u$ and $E^{\prime}.v$ by its parent and removes the edges with the same incident vertices. Consider Iteration 1.1-1.4 in Fig. 2 (c), for $v_{1}$, its smallest neighbor in $E^{\prime}.v$ is $v_{0}$, which is smaller than itself, thus, the parent of $v_{1}$ is replaced by $v_{0}$ in the $\mathsf{parent}$ array. The other vertices are processed similarly (Iteration 1.1 and 1.2). Then, $v_{1}$ is replaced by $v_{0}$ in $E^{\prime}.u$ and $E^{\prime}.v$ in Iteration 1.3, and only edges $(v_{0},v_{1})$ and $(v_{1},v_{0})$ are left after removing the edges with the same incident vertices in Iteration 1.4. In each even iteration, $\mathsf{GPUSCAN}$ conducts the clustering similarly to the odd iteration except the largest neighbor is selected as its parent. Iteration 2.1-2.4 in Fig. 2 (c) show the even case. Note that to obtain the $E^{\prime}.u$ and $E^{\prime}.v$ at the end of each iteration, $\mathsf{GPUSCAN}$ uses $\mathsf{sort()}$ and $\mathsf{partition()}$ functions in the $\mathsf{Thrust}$ library provided by Nvidia [23]. When $E^{\prime}.u$ and $E^{\prime}.v$ become empty, $\mathsf{GPUSCAN}$ explores each vertex and sets its corresponding element in the $\mathsf{parent}$ array as the root vertex in the forest. For example, in Iteration 2.2, the parent of $v_{4}$ is $v_{0}$, the parent of $v_{0}$ is $v_{1}$, and $v_{1}$ is the parent of itself, thus, the parent of $v_{4}$ is set as $v_{1}$. After Phase 2 finishes, the vertices with the same parent are in the same cluster. In Fig. 2 (c), as $v_{0}$, $v_{1}$, $v_{2}$, $v_{4}$, and $v_{7}$ have the same parent, they are in the same cluster.

In Phase 3, $\mathsf{GPUSCAN}$ classifies the hub and outlier vertices. It first retrieves edge $(u,v)$ from $A.u$ and $A.v$ such that $u\in A.u$ is not in a cluster and $v\in A.v$ is in a cluster and keeps them in new arrays $E^{\prime}.u$ and $E^{\prime}.v$. Then, it replaces the vertices in $E^{\prime}.v$ by their parents. Then, for a vertex $u\in E^{\prime}.u$, it has two more different neighbors, then it is a hub vertex. Otherwise, it is an outlier vertex. Consider the example shown in Fig. 2 (d), $v_{3}$, $v_{5}$, $v_{6}$, and $v_{8}$ are not in a cluster, thus, the edges that they form with their neighbors within the cluster are retrieved and stored in $E^{\prime}.u$ and $E^{\prime}.v$ (Step 1). Then, $v_{0}$, $v_{2}$, and $v_{9}$ are replaced with their parents (In Step 2, the parent of $v_{9}$ is $v_{12}$, which is not shown in Fig. 2 due to limited space). Since $v_{8}$ has two neighbors $v_{1}$ and $v_{12}$ in $E^{\prime}.v$, it is a hub vertex. $v_{3}$, $v_{5}$, $v_{6}$ only have one neighbor in $E^{\prime}.v$, they are outlier vertices. $\Box$

• •

  High extra parallelization cost. As shown in [7], the time complexity of the best serial $\mathsf{SCAN}$ algorithm is $O(m\cdot{\mathsf{deg}}_{max})$. On the other hand, Theorem 2 shows that the work of $\mathsf{GPUSCAN}$ is $O(\Sigma_{(u,v)\in E(G)}({\mathsf{deg}}(u)+{\mathsf{deg}}(v))+c\cdot m\cdot\log m)$. Obviously, lots of extra costs are introduced in $\mathsf{GPUSCAN}$ for parallelization.

• •

  Low scalability. $\mathsf{GPUSCAN}$ assumes that the input graphs can be fit in the GPU memory. However, the sizes of real graphs such as social networks and web graphs are huge while the GPU memory is generally in the range of a few gigabytes. The assumption that the sizes of data graphs cannot exceed the GPU memory makes it unscalable to handle large graphs in practice.

$\mathsf{GPUSCAN}$ simply uses the edge array and adjacent array to represent the input graph. The shortcomings of this representation are two-fold: (1) the degree of a vertex can not be easily obtained, which is a commonly used operation during the clustering, (2) the relation between the edge array and adjacent array for the same edge is not directly maintained, which complicates the processing when identifying the cluster vertices in Phase 2. To avoid these problems, we adopt a CSR (Compressed Sparse Row)-enhanced graph layout, which contains four parts:

• •

  Vertex array. Each entry keeps the starting index for a specific vertex in the adjacent array.

• •

  Adjacent array. The adjacent array contains the adjacent vertices of each vertex. The adjacent vertices are stored in the increasing order of vertex ids.

• •

  Eid array. Each entry corresponds to an edge in the adjacent array and keeps the index of the edge in the degree-oriented edge array.

• •

  Degree-oriented edge array. The degree-oriented edge array contains the edges of the graph and the edges with the same source vertex are stored together. For each edge $(u,v)$ in the degree-oriented edge array, we guarantee ${\mathsf{deg}}(u)<{\mathsf{deg}}(v)$ or ${\mathsf{deg}}(u)={\mathsf{deg}}(v)$ and ${\mathsf{id}}(u)<{\mathsf{id}}(v)$. In this way, the workload for edges incident to a large degree vertex can be computed based on the vertex with smaller degree.

• •

  Similarity array. Each element corresponds to an edge in the degree-oriented edge array and maintain the similarity of the corresponding edge.

To reduce the high parallelization cost, the designed algorithm should avoid unnecessary computation during the clustering and fully exploit the unique characteristics of GPU structure. As shown in Section 3.2, $\mathsf{GPUSCAN}$ first computes the structural similarity for each edge in Phase 1 and then identifies the roles of each vertex accordingly. However, based on Theorem 2, computing the structural similarity is a costly operation. To avoid unnecessary structural similarity computation, we adopt a progressive approach in which the lower bound ($\underline{N_{\epsilon}}[v]$) and upper bound ($\overline{N_{\epsilon}}[v]$) of $N_{\epsilon}[v]$ are maintained. Clearly, we have the following lemma on $\underline{N_{\epsilon}}[v]$ and $\overline{N_{\epsilon}}[v]$:

According to Lemma 1, we can maintain the $\underline{N_{\epsilon}}[v]$ and $\overline{N_{\epsilon}}[v]$ progressively and delay the structural similarity computation until necessary. Following this idea, our new algorithm, $\mathsf{GPUSCAN^{++}}$, is shown in Algorithm 1. It first initializes $\underline{N_{\epsilon}}[v]$ as $1$, $\overline{N_{\epsilon}}[v]$ as ${\mathsf{deg}}(v)+1$ for each vertex $v$ (lines 1-2, $\mathsf{VA}$ refers to the vertex array). The role of each vertex and the structural similarity of each edge are initialized as unknown ($\large\footnotesize{\bf{?}}⃝$) (lines 1-4). Then, it identifies the core vertices (line 5), detects the clusters follow the core vertices (line 6), and classifies the hub vertices and outlier vertices (line 7).

To identify the core vertices, it assigns a warp for each edge $(u,v)$ to determine its structural similarity (line 1), which is similar to $\mathsf{GPUSCAN}$. However, the detailed procedures to compute the similarity are significantly different from that of $\mathsf{GPUSCAN}$. Specifically, if the role of $u$ or $v$ has been known, it delays the computation of structural similarity between $u$ and $v$ to the following steps as the role of these two vertices has been determined and it is possible that we can obtain the clustering results without knowing the structural similarity of $u$ and $v$ in the following steps (line 2). Then, it computes the structural similarity between $u$ and $v$ by invoking $\mathsf{checkSim}$. If $u$ and $v$ are similar, then the similarity indicator between $u$ and $v$ is set to $\large\footnotesize{\bfS}⃝$ (line 6) and ${\underline{N_{\epsilon}}[u]}$/${\underline{N_{\epsilon}}[v]}$ is increased by $1$ as $u$ and $v$ are similar (line 7). After that, if ${\underline{N_{\epsilon}}[u]}$/${\underline{N_{\epsilon}}[v]}\geq\mu$, the role of $u$/$v$ is set as $\large\footnotesize{\bfC}⃝$ (lines 8-9). If $u$ and $v$ are dissimilar, ${\overline{N_{\epsilon}}[u]}$/${\overline{N_{\epsilon}}[v]}$ and the role of $u$/$v$ are set accordingly (lines 11-14).

Procedure $\mathsf{checkSim}$ is used to determine the structural similarity between $u$ and $v$ based on the given $\epsilon$. Without loss of generality, we assume ${\mathsf{deg}}(u)<{\mathsf{deg}}(v)$. For each neighbor $w$ of $u$, it assigns a thread to check whether $w$ is also a neighbor of $v$ (line 17). Since the neighbors of $v$ are stored in increasing order based on their id in $\mathsf{AdjA}$, it looks up $w$ among the neighbors of $v$ in a binary search manner (lines 18-26). Specifically, three indices $\mathsf{low}$, $\mathsf{high}$, and $\mathsf{mid}$ are maintained and initialized as ${\mathsf{VA}}[v]$, ${\mathsf{VA[v+1]}}$, and $({\mathsf{low}}+{\mathsf{high}})/2$, respectively. If $w$ is less than ${\mathsf{AdjA}}[{\mathsf{mid}}]$, $\mathsf{high}$ is updated as $\mathsf{mid}$ (lines 21-22) ($\mathsf{AdjA}$ refers to the adjacent array). If $w$ is larger than ${\mathsf{AdjA}}[{\mathsf{mid}}]$, then $\mathsf{low}$ is updated as $\mathsf{mid}$ +1 (lines 23-24). Otherwise, the common neighbors of $u$ and $v$, which are recorded in $\mathsf{sum}$, are increased by 1 (line 26). Last, whether $u$ and $v$ are similar returns based on Definition 2.1 (line 27-29).

Specifically, Algorithm 3 uses the $\mathsf{parent}$ array to store the cluster id that each vertex belongs to (negative cluster id value indicates the corresponding vertex is not in a cluster currently). Similar to $\mathsf{GPUSCAN}$, Algorithm 3 also uses a spanning forest to represent the clusters through the $\mathsf{parent}$ array. However, Algorithm 3 uses a different algorithm to construct the forest, which makes $\mathsf{GPUSCAN^{++}}$ much more efficient than $\mathsf{GPUSCAN}$ as verified in our experiment. It initializes each core vertex $u$ as a separate cluster with cluster id $u$ and the non-core vertices with a negative cluster id (lines 1-3). Then, for each core vertex $u$, it first finds the root id of the cluster by the procedure $\mathsf{root}$ (line 5). After that, for the neighbors $v$ of $u$, if we have already known $u$ and $v$ are similar (line 7), we know $u$ and $v$ are the same cluster following Definition 2.1. Thus, we union the subtrees represented by $u$ and $v$ in $\mathsf{parent}$ through procedure $\mathsf{union}$ (lines 6-9).

After having explored the sub-clusters that can be determined by the computed information, for each core vertices $u$, Algorithm 3 further explores the neighbors $v$ of $u$ such that $v$ is also a core vertex but the similarity of $u$ and $v$ are unknown (lines 10-13). If $u$ and $v$ are not in the same cluster currently (line 15), then it computes the similarity between $u$ and $v$ (line 16). If $u$ and $v$ are similar, we union the subtrees represented by $u$ and $v$ in the $\mathsf{parent}$ array (line 17). Otherwise, $u$ and $v$ are marked as dissimilar (line 18). When the clusters of the core vertices are determined, the root cluster id is assigned to each core vertex (lines 19-20). At last, Algorithm 3 explores the non-core neighbors $v$ of each core vertex $u$, if the similarity between $u$ and $v$ is unknown, their similarity is computed (lines 24-25). If $u$ and $v$ are similar, $v$ is merged in the cluster represented by $u$ (lines 26-27). For the procedure $\mathsf{root}$ and $\mathsf{union}$, they are used to find the root of the tree in the forest and union two sub-trees in the forest. As they are self-explainable, we omit the description. Note the $u.{\mathsf{height}}$/$v.{\mathsf{height}}$ in line 34 represents the height of tree rooted at $u/v$.

$\mathsf{GPUSCAN^{++}}$ first processes the core vertices in parallel. To avoid redundant computation, when clustering core vertices, only neighbors with smaller ids than their own ids are considered. At the same moment, for $v_{1}$, $(v_{1},v_{0})$ are similar, then the parent of $v_{1}$ becomes $v_{0}$ and for $v_{4}$, $(v_{4},v_{1})$ are similar, then the parent of $v_{4}$ becomes $v_{1}$, which are shown in Fig. 5 (a). After that, $\mathsf{GPUSCAN^{++}}$ detects the clusters related to $v_{4}$ and $v_{7}$ as the similarities between some of their neighbors and themself have not been determined yet. For the neighbor $v_{0}$ of $v_{4}$, $\mathsf{GPUSCAN^{++}}$ finds that $v_{0}$ and $v_{4}$ have been in the same cluster by procedure $\mathsf{root}$, they do not need further process, which is shown in Fig. 5 (b). For $v_{7}$, due to the atomic operation, the thread which explores its neighbor $v_{0}$ executes first. It finds that they are not in the same cluster based on procedure $\mathsf{root}$. Therefore, after determining that they are similar by procedure $\mathsf{checkSim}$, it changes the parent of $v_{7}$ to $v_{0}$, which is shown in Fig. 5 (c). Then, $v_{7}$ and $v_{1}$ are processed by another thread similarly, which is shown in Fig. 5 (d). When the cluster of core vertices $v_{0}$, $v_{1}$, $v_{4}$, and $v_{7}$ is detected, $\mathsf{GPUSCAN^{++}}$ explores these vertices and sets their corresponding element in the $\mathsf{parent}$ array as the root vertices in the forest. Fig. 5 (e) illustrates this process, in which the $\mathsf{parent}$ of $v_{4}$ changes to the traced root vertex $v_{0}$. At last, $\mathsf{GPUSCAN^{++}}$ assigns the parent of $v_{1}$ to $v_{2}$ based on the similarity of core vertex $v_{1}$ and non-core vertex $v_{2}$ as shown in Fig. 5 (f). $\Box$