title: GOSH: Embedding Big Graphs on Small Hardware
authors: Taha Atahan Akyildiz Amro Alabsi Aljundi Kamer Kaya
text: 
  
  
  
  
  Both authors contributed equally to this research.
  aakyildiz@sabanciuniv.edu
  Sabanci University
  
  
  [1]
  amroa@sabanciuniv.edu
  Sabanci University
  
  
  Sabanci University
  kaya@sabanciuniv.edu
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  In graph embedding, the connectivity information of a graph is used to represent each vertex as a point in a d-dimensional space. Unlike the original, irregular structural information, such a representation can be used for a multitude of machine learning tasks. Although the process is extremely useful in practice, it is indeed expensive and unfortunately, the graphs are becoming larger and harder to embed. Attempts at scaling up the process to larger graphs have been successful but often at a steep price in hardware requirements. We present , an approach for embedding graphs of arbitrary sizes on a single GPU with minimum constraints. utilizes a novel graph coarsening approach to compress the graph and minimize the work required for embedding, delivering high-quality embeddings at a fraction of the time compared to the state-of-the-art. In addition to this, it incorporates a decomposition schema that enables any arbitrarily large graph to be embedded using a single GPU with minimum constraints on the memory size. With these techniques, is able to embed a graph with over 65 million vertices and 1.8 billion edges in less than an hour on a single GPU and obtains a 93% AUCROC for link-prediction which can be increased to 95% by running the tool for 80 minutes.=-1
  
  
  
  <ccs2012>
  <concept>
  <concept_id>10010147.10010169.10010170.10010171</concept_id> 
  <concept_desc>Computing methodologies Shared memory algorithms</concept_desc>
  <concept_significance>500</concept_significance>
  </concept>
  <concept>
  <concept_id>10002950.10003624.10003633.10010917</concept_id>
  <concept_desc>Mathematics of computing Graph algorithms</concept_desc>
  <concept_significance>500</concept_significance>
  </concept>
  </ccs2012>
  
  
  [500]Parallel computing methodologies Parallel algorithms
  [300]Parallel computing methodologies Shared memory algorithms
  [500]Mathematics of computing Discrete mathematics
  [300]Mathematics of computing Graph theory
  [100]Mathematics of computing Graph algorithms
  
  
  
  GOSH: Embedding Big Graphs on Small Hardware
      Kamer Kaya
      December 30, 2023
  ============================================
  
  
  
  
  
  § INTRODUCTION
  
  
  
  Graphs are widely adopted to model the interactions within real-life data such as social networks, citation networks, web data, etc. Recently, using machine learning (ML) tasks such as link prediction, node classification, and anomaly detection on graphs became a popular area with various applications from different domains. 
  The raw connectivity information of a graph, represented as its adjacency matrix, does not easily lend itself to be used in such ML tasks; regular d-dimensional representations are more appropriate for learning valid correlations between graph elements. Unfortunately, the connectivity information does not have such a structure. Recently, there has been a growing interest in the literature in the graph embedding problem which focuses on representing the vertices of a graph as d-dimensional vectors while embedding its structure into a d-dimensional space.
  
  
  
  Various graph embedding techniques <cit.> have been proposed in the literature. 
  However, these approaches do not usually scale to large, real-world graphs. For example,  <cit.>, deepwalk <cit.>, node2vec <cit.> and  <cit.> require hours of CPU training, even for small- and medium-scale graphs. 
  Although these approaches can be parallelized, a multi-core CPU implementation of takes more than two hours on 16 CPU cores for a graph with 2 million vertices and 20 million edges. 
  There are other attempts to increase graph embedding performance. 
   <cit.> and  <cit.> apply graph coarsening, a process in which a graph is compressed into smaller graphs, to make the process faster, but they do not have a parallel implementation. 
  Accelerators such as GPUs can be used to deal with large-scale graphs. However, to the best of our knowledge, the only GPU-based tool in the literature is . 
  Although is faster than the CPU counterparts, its use is limited by the device memory and it cannot embed large graphs with a single GPU.=-1 
  
  
  In this paper, we present [The code is publicly available at <https://github.com/SabanciParallelComputing/GOSH>]; an algorithm that performs parallel coarsening and many-core parallelism on GPU for fast embedding. The tool is designed to be fast and accurate and to handle large-scale graphs even on a single GPU. However, it can easily be extended to the multi-GPU setting. performs the embedding in a multilevel setting by using a novel, parallel coarsening algorithm which can shrink graphs while preserving their structural information and trying to avoid giant vertex sets during coarsening. With coarsening, the initial graph is iteratively shrunk into multiple levels. Then, starting from the smallest graph, unsupervised training is performed on the GPU. The embedding obtained from the current level is directly copied to the next one by using the coarsening information obtained from the above level. The process continues with the expanded embedding for the next level until the original graph is processed on the GPU and the final embedding is obtained. The contributions can be summarized as follows:
  
      
    * To the best of our knowledge, the only GPU-based graph embedding tool, , cannot handle embedding large graphs on a single GPU when the total size of the embedding is larger than the total available memory of the GPU. On the contrary, applies a smart decomposition, update scheduling, and synchronization to perform the embedding even when the memory requirement exceeds the memory available on a single GPU.  
     
      
    * Thanks to multilevel coarsening and smart work distribution across levels, 's embedding is ultra-fast. 
      For instance, on the graph com-lj, , a state-of-the-art GPU-based embedding tool, spends around 11 minutes to reach 98.33% AUCROC score. On the same graph, spends only 2.5 minutes and obtains 98.46% AUCROC score. Furthermore, based on the numbers given in <cit.>, the embedding takes more than 20 hours on 4 Tesla P100 GPUs for a graph of 65 million vertices and 1.8 billion edges (the link prediction scores are not reported in <cit.> for this graph). On a single Titan X GPU, reaches 95% AUCROC score within 1.5 hours.=-1 
  
      
    * The dimension, i.e., the number of features, used for the embedding process can vary with respect to the application. For different d values, the best possible GPU-implementation, which utilizes the device cores better, also differs. For different d values, performs different parallelization strategies to further increase the performance of embedding especially for small d values.=-1 
     
      
    * The multilevel setting for graph embedding has been previously applied by  <cit.> and  <cit.>. Since CPU-based embedding takes hours, the coarsening literature suggests that the time required for coarsening is negligible in comparison. In this work, we show that a parallel coarsening algorithm is necessary since is orders of magnitude faster than CPU-based approaches. To overcome this bottleneck, we propose an efficient and parallel coarsening algorithm which is empirically much faster than that of and . 
  
  
  The rest of the paper is organized as follows: In Section <ref>, the notation used in the paper is given. 
  Section <ref> describes in detail including the coarsening algorithm and the techniques designed and implemented to handle large graphs. 
  The experimental results are presented in Section <ref> and the related work is summarized comparatively in Section <ref>. 
  Section <ref> concludes the paper.=-1 
  
  
  
  § NOTATION AND BACKGROUND
  
  
  A graph G = (V, E) has V as the set of nodes/vertices and E ⊆ (V × V) as the set of edges among them. For undirected graphs, an edge is an unordered pair while in directed graphs, the order is significant. 
  
  
  
  
  An embedding of a graph G = (V,E) is a |V| × d matrix 𝐌, where d is the dimension of the embedding. The vector 𝐌[i] corresponds to a vertex i ∈ V and each value j in the vector 𝐌[i][j] captures a different feature of vertex i. The embedding of a graph can be used in many machine learning tasks such as link prediction <cit.>,  node classification <cit.> and anomaly detection <cit.>.
  
  
  There are many algorithms in the literature for embedding the nodes of a graph into a d-dimensional space.
  
  implements the embedding method of  <cit.>; a method which, in addition to having fast run-time and low memory overhead, is highly generalizable as it can produce embeddings that reflect any vertex-to-vertex similarity measure Q.
  To elaborate, this approach defines two distributions for each vertex v: The first, sim_Q^v, is obtained from the similarity values between v and every other vertex in G computed based on Q. The second, sim_E^v, is derived from the embedding by using the cosine similarities of v's embedding vector and those of every other vertex in G. A soft-max normalization is applied to these values as a post-processing step so they sum up to 1. The problem then becomes the minimization of the Kullback-Leibler (KL) divergence between the two distributions for every vertex. In this paper, we choose Q to be the adjacency similarity measure described in <cit.>=-1 . 
  
  
  
  
  
  
  
  
  
   
  The training procedure employs Noise Contrastive Estimation for convergence of the objective above as described in <cit.>. The process trains a logistic regression classifier to separate vertex samples drawn from the empirical distribution Q and vertex samples drawn from a noise distribution N, with the corresponding embedding vectors being the parameters of this classifier. More precisely, all the vertices are processed a total number of e times, i.e., epochs, where processing a vertex v ∈ V consists of drawing a single positive sample u from sim_Q^v and n_s negative samples s_1, s_2, …, s_n_s from N. In the meantime, logistic regression is used to minimize the negative log-likelihood of observing u and not observing s_1, s_2, …, s_n_s by updating the corresponding embedding vectors of v and all the other samples. A single update is shown in Algorithm <ref>, where 𝐌[v] is the embedding vector of v, 𝐌[sample] is the embedding vector of the sample, b is a binary number that is 1 if the sample is positive (drawn from Q) or negative (drawn from N), ⊙ is the dot-product operation, σ is the sigmoid function and lr is the learning rate of the classifier.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  The notation used in the paper is given in Table <ref>. 
  
  
  
  
  
  
  
  § EMBEDDING ON SMALL HARDWARE
  
  
  This section is organized as follows: First, a high-level explanation of GOSH is provided. Then, Section 3.1 describes the GPU implementation for the embedding process in detail. Following the embedding process, in Section 3.2, a new coarsening approach is introduced and the parallel implementation details are provided. Finally, Section 3.3 describes the partitioning schema which enables to handle graphs that do not fit on the GPU memory.
  
  
  
  Given a graph G_0, , shown in Algorithm <ref>, computes the embedding matrix 𝐌_0. This is done in two stages; 
  
      
    * creating a set 𝒢 = {G_0, G_1, …, G_D-1} of graphs coarsened in an iterative manner (as in the left of Figure <ref>) where one or more than one nodes in G_i-1 are represented by a super node in G_i (Line <ref>),
      
    * starting from G_D - 1, training the embedding matrix 𝐌_i for the  graph G_i and projecting it to G_i-1 to later train 𝐌_i-1 (as in the right of Figure <ref>) (Lines <ref>- <ref>).
  
  The training process is repeated until 𝐌_0 is obtained. To obtain 𝐌_i-1 from 𝐌_i the mapping information of 𝐆_i-1 is used, where M_i[u] = M_i-1[v] iff u ∈ V_i is a super node of v ∈ V_i-1.
  
  
  
  
  
  provides support for large-scale graphs for which the memory footprint of the training exceeds the device memory. Even for practical sizes, e.g., |V| = 128M and d = 128, the number of entries in the matrix is approximately 16G. With double precision, one needs to have 128GB memory on the device to store the entire matrix. 
  
  For each G_i, initially checks if both G_i and 𝐌_i can fit in the GPU (Line <ref>). If so, it proceeds by copying G_i and the projection of 𝐌_i to the GPU and carrying out the embedding of G_i in a single step (Lines  <ref>-<ref>). Otherwise, it generates positive samples in CPU and carries out the embedding of G_i by copying the relative portions of the samples and 𝐌_i and training the graph in batches (Line  <ref>).
  
  
  Using multilevel coarsening arises an interesting problem; let e be the total number of epochs one performs on all levels. With a naive approach, one can distribute the epochs evenly to each level. However, when more epochs are reserved for G_is in the lower levels, the process will be faster. To add, the corresponding embedding matrices will have a significant impact on the overall process as they are projected to further levels. On the other hand, when more epochs are reserved for the higher levels, i.e., larger G_is, the embedding is expected to be more fine-tuned. So the question is how to distribute the epoch budget to the levels. Based on our preliminary experiments, employs a mixed strategy; a portion, p, of the epochs are distributed uniformly and the remaining (e × (1-p)) epochs are distributed geometrically. That is, training at level i uses e_i = e/D + e'_i epochs where e'_i is half of e'_i+1. 
  
   The value p is called the smoothing ratio and is left as a configurable parameter for the user to establish an interplay between the performance and accuracy.=-1
  
  Another parameter that has a significant impact on embedding quality is the learning rate. 
  For multilevel embedding, a question that arises is how to set the learning strategy for each level. 
  In , we use the same initial learning rate for each level, i.e., for the training of each 𝐌_i and decrease it after each epoch. 
  That is, the learning rate for epoch j at the ith level is equal to lr ×max(1 - j/e_i, 10^-4). 
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Graph embedding in 
  
  
  
  
  implements a GPU parallel, lock-free learning step for the embedding algorithm. As mentioned above, following VERSE, we use an SGD-based optimization process for training. 
  
  
  As shown in <cit.>, a lock-free implementation of SGD, which does not take the race conditions, i.e., simultaneous updates, into account, does not have a significant impact on the learning quality of the task on multi-core processors. However, our preliminary experiments show that on a GPU, where millions of threads are being executed in parallel, such race conditions significantly deteriorate the quality of the embedding. For , we follow a slightly more restricted implementation which is still not race-free. 
  
  
  To reduce the impact of race conditions, we synchronize the epochs and ensure that no two epochs are processed at the same time. For an epoch to train 𝐌_i, traverses V_i in parallel by assigning a single vertex to a single GPU-warp. For a single (source) vertex, multiple positive and negative samples are processed one after another and updates are performed as in Algorithm <ref>. With this implementation and vertex-to-warp assignment, a vertex v ∈ V_i cannot be a source vertex for two concurrent updates but it may be (positively or negatively) sampled by another vertex while the warp processing v is still active. Similarly, v can be sampled by two different source vertices at the same time. Hence, the reads/writes on 𝐌_i[v] are not completely race free. However, according to our experiments, synchronizing the epochs and samples for the same source vertex is enough to robustly perform the embedding process. 
  
  As shown in Algorithm <ref>, the positive and negative samples are generated in the GPU. For a source vertex v ∈ V_i, the positive sample u ∈ V_i is chosen from Γ_G_i(v). 
  Negative samples, on the other hand, are drawn from a noise distribution as mentioned in Section <ref> which we model as a uniform random distribution over V_i. 
  
  After each sampling, the corresponding update is performed by using the procedure in Algorithm <ref>. 
  
  
  
  During the updates for a source vertex src, the threads in the corresponding warp perform (1 + n_s) × d accesses to 𝐌_i[src]. For large values of n_s and d, performing this many global memory accesses dramatically decreases the performance. To mitigate this, before processing src, copies 𝐌_i[src] from global to shared memory. Then for all positive and negative samples, the reads and writes for the source are performed on the shared memory. Finally, 𝐌_i[src] is copied back to global memory. On the other hand, for the sampled u vertices, 𝐌_i[u] is always kept in global memory since each entry is read and written only once. To perform coalesced accesses on 𝐌_i[u], the reads and writes are performed in a round-robin fashion. That is 𝐌_i[u][j + (32 × k)] is accessed by thread j at the kth access where 32 is the number of threads within a warp.  
  
  
  
    §.§.§ Embedding for small dimensions
   
  
  Assuming a warp contains 32 threads, when d ≤ 16 dimensions are used for embedding, a single source vertex does not keep all the threads in a warp busy. In this case, 32 - d warp threads remain idle which yields to the under-utilization of the device. To tackle this problem, we integrate a specialized implementation for small dimension embedding. We set the number of threads responsible for a source vertex as the smallest multiple of 8 larger than or equal to d, i.e., 8 or 16. Hence, depending on d, we can assign 2 or 4 vertices to a single warp. 
  
  
  
   §.§ Graph coarsening
  
  
  
  
  
  
  
  
  
  
  
   
  
  employs a fast algorithm to keep the structural information within the coarsed graphs while maximizing the coarsening efficiency and effectiveness. Coarsening efficiency at the ith level is measured by the rate of shrinking defined as (|V_i - 1| - |V_i|)/|V_i - 1|. On the other hand, the effectiveness is measured in terms of its embedding quality compared to other possible coarsenings of the same graph embedded with the same parameters. 
  We adapt an agglomerative coarsening approach, , which generates vertex clusters in a way similar to the one used in <cit.>. Given G_i = (V_i, E_i), the vertices in V_i are processed one by one. If v is not marked, it is marked, and mapped to a cluster, i.e., a new vertex in V_i+1 and its edges are processed. If an edge (v,u) ∈ E_i, where u is not marked, u is added to v's cluster. Then, all of the vertices in v's cluster are shrunk into a super vertex v_sup∈ G_i + 1.=-1
  
  
  preserves both the first- and second-order proximites <cit.> in a graph. The former measures the pairwise connection between vertices, and the latter represents the similarity between vertices' neighborhoods. It achieves that by collapsing vertices that belong to the same neighborhood around a local, hub vertex.
  However, if this process is handled carelessly, two, giant hub vertices can be merged. We observed that this degrades the effectiveness and efficiency of the coarsening. The effectiveness degrades since the structural equivalence is not preserved in the lower levels of the coarsening, where most of the vertices are represented by a small set of super vertices. Furthermore having a small set of giant supers inhibits the graph from being coarsened further, resulting in an insufficient efficiency. To mitigate this, a new condition for matching is introduced to the algorithm, where u ∈ V_i can not be put into the cluster of v ∈ V_i if |Γ_G_i(u)| and |Γ_G_i(v)| are both larger than |E_i|/|V_i|. Consequently, assuming that the hub vertices will have a higher degree than the density of G_i, two of them can no longer be in the same cluster. Our preliminary experiments show that this simple rule has a significant effect on both the efficiency and the effectiveness of the coarsening.
  
  As mentioned above, when a vertex is marked and added to a cluster, its edges are not processed further and it does not contribute to the coarsening. Performing the coarsening with an arbitrary ordering may degrade the efficiency since large vertices can be locked by the vertices with small neighborhoods. Hence, when an edge (u, v) ∈ E_i is used for coarsening for a hub-vertex v ∈ V_i, to maximize efficiency, we prefer u ∈ V_i to be inserted in to the cluster of origin v. To provide this, an ordering is procured by sorting the vertices with respect to their neighborhood size and this ordering is used during coarsening. This ensures processing vertices with a higher degree before the vertices with smaller neighborhoods and this results in a substantial increase in the coarsening efficiency.
  
  The details of the coarsening phase are given in Algorithm <ref>. The algorithm takes an uncoarsened graph G = G_0 and returns the set of coarsened graphs 𝒢 along with the mapping information to be used to project the embedding matrices ℳ. 𝒢 and ℳ are initialized as {G_0} and empty set, respectively. Starting from i = 0, the coarsening continues until a graph G_i+1 with less than threshold vertices is generated. As mentioned above, first the vertices in G_i are sorted with respect to their neighborhood sizes. Then the coarsening is performed and a smaller G_i+1 is generated. We also store the mapping information map_i used to shrink G_i to G_i+1. This will be used later to project the embedding matrix M_i+1 obtained for G_i+1 to initialize the matrix M_i for G_i. 
  
  To add, threshold = 100 is used for all the experiments in the paper which is the default value for . 
  
  
  
  
  
    §.§.§ Complexity analysis:
   
  All the algorithms, coarsening and embedding, use the Compressed Sparse Row (CSR) graph data structure. In CSR, an array, adj holds the neighbors of every vertex in the graph consecutively. It is a list of all the neighbors of vertex 0, followed by all the neighbors of vertex 1, and so on. Another array, xadj, holds the starting indices of each vertex's neighbors in adj, with the last index being the number of edges in the graph. In other words, the neighbors of vertex i are stored in the array adj from adj[xadj[i]] until adj[xadj[i+1]]. 
  
  has three stages; sorting (line <ref>), mapping (lines <ref>–<ref>) and coarsening (line <ref>). A counting sort is implemented for the first stage with a time complexity of 𝒪(|V| + |E|). For mapping, the algorithm traverses all the edges in the graph. This has a time complexity of 𝒪(|V| + |E|). Finally, coarsening the graph requires sorting the vertices with respect to their mappings and going through all the vertices and their edges within the CSR, which also has a time complexity of 𝒪(|V| + |E|).=-1
  
  
  
    §.§.§ Parallelization:
   
  When the embedding is performed on the CPU, as the literature suggests, embedding dominates the total execution time. However, with fast embedding as in , this is not the case. Thus, we parallelize the coarsening on the CPU.
  
  In a parallel coarsening with τ threads, one can simply traverse V_i in parallel and perform the mapping with no synchronization. However, this creates race conditions and inconsistent coarsenings.
  
  
  To avoid race conditions, we use a lock per each entry of map_i. To update map_i[v] and map_i[u] as in lines <ref> and <ref>, the thread first tries to lock map_i[v] and map_i[u], respectively. If the lock is obtained, the process continues. Otherwise, the thread skips the current candidate and continues with the next vertex. One caveat is the update on the counter cluster. Hence, instead of using a separate variable for super vertex ids, the parallel version uses the hub-vertex id for mapping. That is map_i[v] is set to v unlike line <ref> of the sequential algorithm. Note that with this implementation, map_i does not provide a mapping to actual vertex IDs in G_i+1. This can be fixed in 𝒪(|V|) time via sequential traversals of the map_i array, which first detect/count the vertices that has map_i[v] = v and reset the map_i values for all.   
  
  
  
  
  
  
  
  The parallel coarsened graph construction is not straightforward. After the mapping, the degrees of the (super) vertices in G_i+1 are not yet known. To alleviate that, we allocate a private E^j_i+1  region in the memory to each thread t_j, 1 ≤ j ≤τ. These threads create the edge lists of the new vertices on these private regions which are then merged on a different location of size |E_i+1|. To do that first a sequential scan operation is performed to find the region in E_i+1 for each thread. Then the private information is copied to E_i+1.  
  
  An important problem that needs to be addressed for all the steps above is load imbalance. Since the degree distribution on the original graph can be skewed and becomes more skewed for the coarsened graphs, a static vertex-to-thread assignment can reduce the performance. Hence uses a dynamic scheduling strategy, which uses small batch sizes for all the steps above.
  
  
  
   §.§ Handling large graphs
  
  One of the strongest points of is the ability to quickly generate a high-quality embedding of a large-scale graph on a single GPU. Here we present the techniques used when the graph and the embedding matrix do not fit in the GPU memory. The base algorithm for requires storing d × |V| and (|V| + 1) + |E| entries for the matrix and the graph, respectively, for a graph G =(V,E). For large graphs with millions of vertices, a single GPU is unable to store this data.=-1
  
  
  
  We mitigate the bottleneck of storing 𝐌_i on the GPU by using a partitioning schema similar to <cit.>, in which 𝐌_i is partitioned and embedding is carried out on the sub-matrices instead of the entirety of 𝐌_i. Formally, we partition V_i into K_i disjoint subsets of vertices 𝒱_i = {V_i^0, V_i^1, …, V_i^K_i-1}. 
  
  Let 𝒫_i = {𝐌_i^0, 𝐌_i^1, …, 𝐌_i^K_i-1} be the sub-matrices of 𝐌_i corresponding to the vertex sets in 𝒱_i.
  With partitioning, embedding G_i becomes the process of moving the sub-matrices in 𝒫_i to the GPU, carrying out the training on these parts, and switching them out for the next sub-matrices, and so on.
  
  In order to carry out the embedding correctly, all possible negative samples, chosen from  V_i× V_i, must be able to be processed. To do this, handles the embedding in rotations. During a rotation there will always be a point in time when 𝐌_i^j and 𝐌_i^k are together in the GPU for all 0 ≤ j < k < K. When this happens, B positive and B × n_s negative samples are chosen from V_i^j (or V_i^k) for every vertex in V_i^k (or V_i^j).
  
  This way, in each rotation, we run a total of at most B × K_i updates on every vertex, which makes a full rotation (almost) equivalent to B × K_i epochs. We use the term almost since a vertex v ∈ V_i^j may not have a neighbor in V_i^k.
  
  In this case, no positive updates are performed for v whenV_i^j and V_i^k are on the device. Hence, running e_i/B × K_i rotations is (almost) equivalent to running e_i epochs for the embedding. We set B = 5 as the default value in and experiment with different values to see its impact on embedding quality and performance.=-1
  
  
  
  Although partitioning the embedding matrix solves the first memory bottleneck, storing G_i on the GPU can still be problematic since this will leave less space for the sub-matrices. This is why we opt not to store G_i on the GPU. The positive samples are chosen on the host and only when required, they are sent to the GPU. 
  
  Negative samples, on the other hand, are still generated on the GPU: the kernel for the parts (V_i^j, V_i^k) draws the negative samples for vertices in V_i^j randomly from V_i^k and vice versa for vertices in V_i^k.=-1
  
  
  
  
  
  
    §.§.§ Minimizing the data movements:
  
  Since the embedding is performed on pairs of sub-matrices, the rotation order in which we process the pairs is important to minimize the data movement operations. We follow an order resembling the inside-out order proposed in <cit.> which
  
  formally defines the part pairs as (V_i^a_0, V_i^b_0), (V_i^a_1, V_i^b_1), ⋯, (V_i^a_ℓ, V_i^b_ℓ) where ℓ = K_i(K_i+1)/2 and
  
      (V_i^a_j, V_i^b_j) = 
      
      (V_i^0, V_i^0)    j = 0
      (V_i^a_j-1, V_i^b_j-1+1)     j > 0  and  a_j-1 > b_j-1
      
      (V_i^a_j-1+1, V_i^0)     a_j-1 = b_j-1
  
  
  
  
  
    §.§.§ Choosing the sub-matrices and samples to be stored:
   
  
  Let P_GPU be the number of sub-matrices stored on the GPU at a time. Since we require every sub-matrix pair to exist on the GPU together during a single rotation, the smallest acceptable value is 2. However, P_GPU = 2  means that there will be time instances where all the kernels processing the current sub-matrices finish and a new kernel cannot start until a new sub-matrix is copied to the GPU. This leaves the GPU idle during the copy operation. On the other hand, using P_GPU > 2 occupies more space but allows an overlap of data transfers with kernel executions. For instance, assume 𝐌_i^1, 𝐌_i^2 and 𝐌_i^4 are on GPU and the three upcoming kernels are (V_i^4, V_i^1), (V_i^4,V_i^2) and (V_i^4,V_i^3). The first two kernels are dispatched and after the first finishes, while the second is running, 𝐌_i^1 is replaced with 𝐌_i^3, thus hiding the latency. 
  
  A large P_GPU increases the amount of overlap. However, it also consumes more space on the GPU and increases K_i, i.e., the number of sets in 𝒱_i. This leads to a rotation containing more kernels, i.e., pairs to be processed. 
  
  For this reason, we set P_GPU = 3 to both hide the latency and occupy less GPU memory.=-1
  
  Since we do not keep the large graphs on GPU memory and draw positive samples on the CPU, these samples must also be transferred to the GPU for each kernel. However, if all these samples are transferred at once, similar to above, K_i increases and the performance decreases. To solve this issue, we only keep samples for S_GPU pairs on the GPU and dynamically replace a pool once it is consumed. We set S_GPU = 4 as we've experimentally found it to be an adequate value for all our graphs. 
  
  
  
  
  
  
  
  
  
  
    §.§.§ Implementation details:
  
  Embedding a large graph G_i requires an orchestrated execution of multiple tasks on the host and the device. coordinates the following tasks by single host thread:=-1 
  
  
    * The main thread dispatches the embedding kernels to the GPU, as well as moves the sub-matrices forth and back between the host and GPU. Multiple GPU streams are used to allow for multiple kernel dispatches at once to maximize the utilization.
  
    * The thread performs (positive) sampling into pools. When necessary, creates a team of threads to generate samples for a single sample pool. Once a pool is 
  
  ready to be sent to the GPU, it is kept in a buffer.
  
    * dispatches the sample pools to the GPU. Once a sample pool is used up on the GPU, and becomes free to be overwritten, dispatches a ready pool to the GPU. 
  
  
  
  Coordination among these threads is provided through condition variables. 
  A high-level overview of the large graph embedding process is shown in Figure <ref> and its pseudocode is given in Algorithm <ref>.=-1
  
  In Line <ref> of Algorithm <ref>, we compute K_i and the number of rotations e`. Line <ref> initializes an array GPUState of size P_GPU which keeps track of embedding sub-matrices currently stored on the GPU. GPUState[j] = k means that bin j on the GPU currently holds M_i^k where an entry -1 indicates that the bin j is empty. At the beginning, the first P_GPU sub-matrices are copied to the GPU (lines <ref>–<ref>) by calling the function (j, k) which copies 𝐌_i^j out of the GPU, replaces it with 𝐌_i^k, and returns the new GPUState. Afterwards, the and threads are started and the main thread starts to perform the embedding rotations. Lines <ref>–<ref> run the main embedding loop. 
  
  When the sub-matrices 𝐌_i^m, 𝐌_i^s and the sample pool S_i^m, s are ready, the function (𝐌_i^m, 𝐌_i^s, n_s, lr) runs the embedding kernel on the sub-matrices pair (𝐌_i^j, 𝐌_j^k) using n_s negative samples and learning rate lr. Finally, the function (GPUState, a, b) will determine the next sub-matrix to be switched into the GPU after running (𝐌_i^j, 𝐌_i^k, n_s, lr) given GPUState.
  
  
  
  
  
  § EXPERIMENTS
  
  
  
  We will first explain the ML pipeline we used to evaluate embeddings and to compute AUCROC (Area Under the Receiver Operating Characteristics) scores. Then the data-sets will be summarized and the state-of-the-art tools used to evaluate the performance of will be listed. Lastly, the results will be given. 
  
  
  
  
   §.§ Evaluation with link-prediction pipeline
  
  
  We evaluate the embedding quality of , , , and with link prediction, which is one of the most common ML tasks that  embedding algorithms are evaluated by <cit.>.
  
  First, the input graph G is split into train and test sub-graphs as G_train = (V_train, E_train) and G_test = (V_test, E_test) respectively. G_train contains 80% of the edges of G, where G_test contains the remaining 20%. Then, we remove all the isolated vertices from G_train and also all (u,v) edges from G_test, where u or v is not in G_train.
  
  This guarantees that V_test⊆ V_train.
  
  Next, we execute the tools to generate embeddings of G_train.
  
  
  Finally, we employ a Logistic Regression model which uses the embeddings generated in the previous step to predict the existence of edges in G_test. For medium scale graphs, we used the LogisticRegression module from scikit-learn. However, logistic regression becomes too expensive for large-scale graphs. Thus, for such graphs, we use the SGDClassifier module from scikit-learn with a logistic regression solver. 
  
  
  
  For the prediction pipeline, we create two matrices, 𝐑_train, and 𝐑_test. 
  Each vector 𝐑_train[i] represents either a positive or a negative sample, and we obtain these vectors by doing element-wise multiplication of two vectors from 𝐌_0, corresponding to two vertices in the graph.  
  𝐑_train includes all the edges in  G_train as positive samples. 
  Moreover, we generate |E_train| number of negative samples from (V_train× V_train) ∖ E_train and add them as vectors to 𝐑_train to make a balanced training set for the logistic regression classifier. 
  In addition to d values obtained via element-wise multiplications, for 𝐑_train, a label representing a positive or negative sample is concatenated to the end of the vector. Hence, the length of each 𝐑_train vector is d+1. 
  
  We create 𝐑_test in a similar fashion by using G_test instead of G_train as the source of the samples. We first train the logistic regression model using 𝐑_train, and then test the validity of the model with 𝐑_test. Finally, we report the AUCROC score obtained from the test set <cit.>.=-1
  
  
  
  
  
  
   §.§ Datasets used for the experiments
  
  We use various graphs in the evaluation process to cover many structural variations and to evaluate the tools in terms of performance and quality as fairly and thoroughly as possible. The graphs differ in terms of their origin, the number of vertices, and density. The properties of these graphs are given in Table <ref>. The medium-scale graphs, with less than 10M vertices, and large-scale ones are separated in the table.
  
  
  
  
  
  
   §.§ State-of-the art tools used for evaluation
  
  To evaluate the performance and the quality of , we use the results of the following state-of-the-art tools as a baseline.
  
  
  
  
  : is a recent multi-core graph embedding tool <cit.>. It can employ different vertex-similarity measures including PPR, adjacency lists, and SimRank. We use the PPR similarity measure and α = 0.85 as recommended by the authors. For , we set the epoch number to e = 600, 1000, and 1400, use a learning rate of lr = 0.0025 and report the best AUCROC. Larger learning rates produce worse results. 
  
  
  : is a state-of-the-art, fast multi-GPU graph embedding tool. However, according to <cit.>, the algorithm cannot embed graphs with |V| > 12,000,000 on a single GPU. We use the default values for the hyperparameters as recommended by the authors and LINE is chosen as the base embedding method. As for the number of epochs, we use two settings; a fast setting with e = 600 epochs, and a slow setting with e = 1000 epochs.
  
  
  : performs embedding by coarsening a graph into multiple levels similar to  <cit.>. It trains the smallest level and refines the embeddings up the coarsening levels. However, unlike , uses a neural network to project the coarsened graph embedding to that of the original graph. We use the following parameters for the model: DeepWalk as a base embedding method, MD-GCN as a refinement method, 8 levels of coarsening, and a learning rate of lr = 0.001. As for the epochs used during training, we note that does not allow the number of epochs to be configured. Hence, that decision was left to the model itself. 
  
  
  
  
  
  
  
  
  For , we use three configurations: fast, normal and slow, with parameters given in Table <ref>. The configurations differ in terms of the number of epochs, smoothing ratio and learning rate. Compared to -slow, -fast uses a smaller p and a larger lr to compensate for the lesser number of epochs on the original graph with faster learning. Furthermore, we include in the experiments a version of which does not perform coarsening. This configuration spends all of the epochs on the original graph. In addition to these differences, for medium-scale graphs, a larger number of epochs is used for each configuration compared to large-scale graphs. 
  
  For the experiments with and , we define a single epoch as sampling |E| target vertices. We do so to match the definition of an epoch given by  <cit.> for the fairness of the experiments. With this definition, using fewer epochs for large-scale graphs makes more sense since a single epoch implies billions of samples and updates for such graphs, whereas for medium-scale graphs this number is in the order of tens of millions.
  
  
  
  
  All the experiments run on a single machine with 2 sockets, each with 8 Intel
  
  E5-2620 v4 CPU cores running at 2.10GHz with two hyper-threads per core (32 logical cores in total), and 198GB RAM. To avoid the effects of hyper-threading, we only use 16 threads for parallel executions. The GPU experiments use a single Titan X Pascal GPU with 12GB of memory. 
  
  The server has Ubuntu 4.4.0-159 as the operating system. All the CPU codes are compiled with gcc 7.3.0 with -O3 as the optimization parameter. For CPU parallelization, OpenMP multithreading is used in general. Only for large graphs, a hybrid implementation with OpenMP and C++11 threading is employed. For GPU implementations and compilation, nvcc with CUDA 10.1 and optimization flag -O3 are used. The GPUs are connected to the server via PCIe 3.0 x16. For GPU implementations, all the relevant data structures are stored on the device memory, unified memory is not used.
  
  
  
   §.§ Experiments on coarsening performance
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Table <ref> provides the properties of the coarsenings obtained from the sequential and parallel coarsening algorithms with τ = 32 threads. 
  As the results show, parallel coarsening reaches a similar number of levels and the graphs at the last-level are of similar sizes. 
  Hence, there is a negligible difference regarding the quality of graphs generated by the two algorithms. Only for soc-sinaweibo, there is a one-level difference for which |V_D - 1| also has a difference of 142 vertices.=-1
  
  With a similar coarsening quality, the parallel algorithm is 5–10× faster compared to the sequential counterpart. As described in Section <ref>, the 
  
  time complexity is 𝒪(|V| + |E|) and in practice, |E| dominates the workload. Although there are other parameters, the variation in the speedups is in concordance with the variation in the number of edges. For instance, soc-sinaweibo only has 200M edges and yields the smallest speedup value of 5.8×. On the other hand, the largest speedup 10.5× is obtained for com-friendster, which is the largest in our data-set with 1.8B edges.=-1 
  
  
  
  
  
  
  
    §.§.§  vs 
  
  
  In Table <ref>, we show a brief comparison of MILE and with 16 threads on the graph com-orkut. Since MILE does not have a stopping criterion for coarsening, we used the same amount of coarsening levels for both algorithms. While coarsening a graph of 3 million vertices and 100 million edges, is 264 times faster than MILE. Moreover, is a lot more efficient regarding the number of vertices obtained at each level. For instance, in 8 levels shrinks com-orkut to only 230 vertices, while MILE shrinks it to 12062 vertices. This is important since the training time is affected by the number of vertices at each level.=-1
  
  
  
  
  
  
  
  
   §.§ Experiments on handling large graphs
  
  
  
  Figure <ref> shows the effect of adjusting the batch size B for large-graph embedding. The top figure provides the execution time and the bottom figure provides the AUCROC scores. We can see the trade-off between performance and quality while increasing B. The execution time decreases since the number of embedding rounds is reduced. 
  
  However, the quality also decreases since 
  
  increasing B results in increasing the number of updates being carried out on a subset of the graph in isolation from the rest of the embedding process. 
  
  To be as efficient as possible and not to decrease the accuracy significantly, we use B = 5 as the default value for . 
  
  
  
  
   §.§ Experiments on embedding quality
  
  
  Tables <ref> and <ref> provide the execution times and AUCROC scores of the tools evaluated in this work on medium-scale and large-scale graphs, respectively. The parameters for the tools are given in Section <ref> (see Table <ref> for configurations).
  
  The first observation is that coarsening does not have a significant negative effect on the quality of the embedding. While performs worse on some medium-scale graphs with coarsening, for others, the scores are approximately the same or even better. Since all of the epochs are reserved for the original graph in the non-coarsened version, we expect the results of this configuration to be more fine-tuned. However, the data shows that training in the coarser levels may have a more prominent effect than training more on the original graph.=-1
  
  
  
  In addition to configurations, we use , , and on medium-scale graphs. The results are given in Table <ref>.
  
  
  To evaluate the runtime performance of the tools, we use the execution time of as the baseline and present the speedups by each tool. From these experiments, we have the following observations:
  
      
      
    * -fast is an ultra-fast solution that produces very accurate embeddings at a fraction of the time compared to all the systems under evaluation. It can achieve a speedup over of up to three orders of magnitude and an average of 600× with a maximum loss in AUCROC of 2% and an average loss of 1.16%. When compared to , it is superior in terms of AUCROC in six out of eight graphs while being at least two orders of magnitude faster. As for , we see that, at an average loss in AUCROC of 0.54%, can achieve an average speedup of 23.44×. 
      
    * -normal demonstrates the speed/quality trade-off of and its flexibility. Switching from -fast to -normal results in an average AUCROC increase of 0.76% while only reducing the speed on average by a factor of 3×.=-1
      
    * The flexibility is demonstrated further by -slow whose accuracy becomes close to the best tool for every graph. Compared to , this configuration has an average loss of 0.24% in AUCROC, but still has an average speedup of 79.24×.=-1   
      
    * To compare with the state-of-the-art GPU implementation , we used the best AUCROCs of the tools. For 4/8 graphs, configurations produce better AUCROC compared to configurations. The values are similar; on average, achieves 0.16% higher AUCROCs than . However, is 5.2× faster than on average.
  
  
  
  
    §.§.§ Large-scale graphs
  
  
  The results of the experiments for large-scale graphs can be seen in Table <ref>. We observe the following: 
  
  
  results are not reported since, for all the large-scale graphs, the executable runs out of CPU memory on our machine. We find that on hyperlink2012 is reported to achieve 94.3% link-prediction AUCROC after an embedding for 5.36 hours using four Tesla P100s GPUs <cit.>. -normal achieves an AUCROC of 97.20% after an embedding taking only 0.2 hours using a single Titan X GPU (26.8× speedup). It is also reported that takes 20.3 hours on com-friendster <cit.> where -normal requires only 0.76 hours (26.7× speedup). 
  
  
  cannot embed hyperlink2012 and soc-sinaweibo before the 12 hour timeout. Furthermore, the executable failed to run on the other two graphs due to insufficient memory.
  
  
  timed out for 3 of the 4 graphs, as shown in Table <ref>. However, it performed the embedding on soc-sinaweibo  successfully. Compared to -slow, it scores a 0.52% higher AUCROC. On the other hand, -slow achieves a 26× speedup over . 
  
  
  
  
  
  
  
  
  
  
   §.§ Experiments on smaller dimensions
  
  
  We analyzed the performance of when multiple vertices are assigned to a single warp, where d is small. 
  The results on com-orkut and soc-LiveJournal are given in Table <ref>.
  Without small-dimension technique (SM), takes approximately the same time for d = 8, 16 and 32 where 4× and 2× less work is performed for d = 8 and 16. With SM, we observe 2.63× and 1.84× speedups for d = 8 and 16, respectively. Moreover, for soc-LiveJournal, we obtain 2.70× and 1.85× speedups for d = 8 and d = 16, respectively. As expected, with or without SM, d = 32 timings are almost the same.=-1 
  
  
  
  
  
  
   §.§ Speedup breakdown
  
  
  
  For a more detailed analysis of performance improvements, we run intermediate versions of and report the speedup over 16-thread CPU implementation. The experiments are conducted with six graphs; two large-scale graphs (com-friendster, and hyperlink2012), and four medium-scale graphs. We did not run the GPU implementations that are not using coarsening on the large-scale graphs as they take a long amount of time. The results are presented in Figure <ref>.
  
  
  
  The first version is the Naive GPU implementation that results in an average slowdown of 3.3×. The Optimized GPU version leverages architecture-specific optimizations to reduce memory access overhead. That is, global memory is organized to have coalesced accesses, and shared memory is utilized to reduce the number of global memory accesses. This version is 5.4× faster than the 16-thread one. These two versions do not use coarsening.=-1 
  
  The next version employs Sequential Coarsening, as well as all the GPU optimizations from the previous one. This version scores an average speedup of 45× over the CPU version, while maintaining the embedding quality as shown in Table <ref>. This is due to the cumulative nature of the updates on the coarsened graphs, where a single update on a super vertex is propagated to all the vertices it contains.=-1  
  
  The Parallel Coarsening version, which is the final , further improves the performance. As discussed in Section <ref>, the performance difference between Parallel Coarsening and Sequential Coarsening is expected to be more for larger graphs. For instance, on com-friendster, sequential and parallel coarsening phases take 2468.52 and 235.38 seconds, respectively (Table <ref>). On the same graph, the total run-time of -normal is 2720.82 seconds (Table <ref>). In other words, parallel coarsening results in an 80% improvement on performance.
  
  
  
  § RELATED WORK
  
  There are various approaches proposed for embedding graphs. A survey on these techniques can be found in <cit.>. For instance, in matrix factorization based embedding, the matrix representing the relationship between vertices in the graph is factorized <cit.>. Another approach in the literature is the sampling-based embedding <cit.> in which the samples are drawn from the graph and used to train a single-layer neural network. Different embedding algorithms use different sampling strategies, most notable of which are random walks <cit.>. 
  
  Each algorithm in the literature tries to capture the structural and role/class information based on a similarity measure and/or a sampling strategy. There also have been attempts <cit.> for a more generalized embedding process. 
  
  Graph embedding is an expensive task and hard to perform without high-performance hardware. 
  Graph coarsening <cit.>, as well as distributed system approaches <cit.> have been previously used to tackle this issue. 
  
  However, these methods do not fully utilize a specialized, now ubiquitous piece of hardware, GPU. The literature does not usually focus on the performance of coarsening since, on CPUs, embedding is slow and coarsening time is negligible. However, this is not the same for GPUs. Furthermore, even state-of-the-art embedding tools are using coarsening schemes that cannot shrink the graphs well.=-1
  
  There have been attempts to make the embeddings faster by utilizing GPUs; generates samples on the host device and performs the embedding on the GPU <cit.>. However, when the total GPU memory is not capable of storing the embedding matrix, the tool cannot perform the embedding. To the best of our knowledge, 
  there is no sampling-based graph embedding tool in the literature which is designed for memory restricted but powerful GPUs. tries to utilize the power of GPUs by applying coarsening and using a judiciously orchestrated CPU-GPU parallelism.=-1
  
  
  
  § CONCLUSION AND FUTURE WORK
   
  
  In this paper, we introduce a high-quality, fast graph embedding approach that utilizes CPU and GPU at the same time. The tool can embed any graph by employing a  partitioning schema along with dynamic on-CPU sampling. In addition to this, we provide optimization techniques to minimize GPU idling and maximize GPU utilization during embedding. Furthermore, a parallel coarsening algorithm, which outperforms the state-of-the-art coarsening techniques both in terms of efficiency and speed, is proposed. We fine-tuned the coarsening approach to produce high-quality embeddings at a fraction of the time spent by the state-of-the-art. Our experiments demonstrate the effectiveness of on a wide variety of graphs.
  
  In the future, we are planning to make publicly available with easy-to-use interfaces for widely used software such as Matlab and scikit-learn. Furthermore, we will extend our work for other ML tasks such as classification and anomaly detection.=-1
  
  
  
  
  
  
  
  
  
  ACM-Reference-Format
  
  
  