HBMax: Optimizing Memory Efficiency for Parallel Influence Maximization on Multicore Architectures

Let $G=(V,E)$ be a graph with $n$ vertices and $m$ edges. Given $G$ and a stochastic diffusion process, the influence maximization problem is one of identifying a set $S$ of top $k$ vertices in $G$ (called “seeds”) that maximizes the expected influence spread ($\mathbb{E}[I(S)]$) in $G$, measured by the number of activated vertices. Kempe et al. (Kempe et al., 2003) have shown that $\mathbb{E}[I(S)]$ is a submodular function of $S$ for two simple but powerful diffusion models: the independent cascade model (IC) and the linear threshold model (LT). Submodular functions have the property of diminishing marginal gains and leads to efficient approximation algorithms (Nemhauser et al., 1978). Leveraging the submodularity framework, they propose a greedy hill-climbing algorithm that offers an approximation guarantee of $1-1/e-\epsilon$ (Kempe et al., 2003).

The IC and LT diffusion models are broadly studied in the literature. The IC model comes from the physics of interacting particles. In this model, each newly activated vertex has a single chance to activate its neighbors. Assuming a directed graph $G=(V,E,\omega)$ as an example, each edge $(u,v)\in E$ is assigned with a probability $p(u,v)$ to trigger the activation of $v$ from $u$. When the vertex $u$ gets activated at time $t$, then at time $t+1$, the vertex $v$ gets activated with the probability $p(u,v)$.

In contrast, the LT model captures mass behaviours. Each vertex has a threshold modeling their resistance to adopting the mass behavior and each edge $(u,v)$ has a weight $w$ representing the capacity of $u$ to influence $v$. At time $t$, a vertex $v$ gets activated if the sum of the weights on the incoming edges from its already active neighbors exceeds its threshold. It is worth noticing that the edge weights $\omega$ are provided by input graphs, but they are not corresponding to the probability of activation. Thus they only affect the LT model instead of the IC model.

The LT diffusion model tends to produce very small RRR sets and more than $99\%$ of them end up with less than 10 vertices (Minutoli et al., 2020a). On the contrary, the sizes of RRR sets from the IC model is less skewed and tend to be of larger size. From a utility standpoint, the IC model is more generally applicable to abstract a range of diffusion processes. It can be viewed as a special case of the Susceptible-Infected-Recovered (SIR) models (Minutoli et al., 2020b) used in epidemics. However, its large memory footprint (as shown in Figure 1) limits its scalability, so we focus on the IC diffusion model in this paper.

Reverse Influence Sampling. Following the seminal work (Kempe et al., 2003), Borgs et al. (2014) proposed an approximation algorithm based on the idea of Reverse Influence Sampling (RIS). Their scheme attempts to find which are the most likely causes of activation for each vertex in the graph. The algorithm randomly samples vertices $v$ and simulates the diffusion model in reverse collecting sets of vertices which may cause the activation of $v$. The seed sets $S$ is later decided by solving a maximum coverage problem over the sets collected through RIS. This approach provides the same approximation guarantee as the greedy hill-climbing algorithm (Kempe et al., 2003).

The RIS approach is the fundamental building block of the IMM algorithm from Tang et al. (2015), and its state-of-the-art parallel implementation Ripples (Minutoli et al., 2019) is the starting point for our work. In the following, we give important definitions and a high-level description of the fundamental building blocks of the IMM and Ripples algorithm. We direct the reader to the original work for a more exhaustive presentation.

In the following discussion, we use RRR sets to refer to the collection of random reverse reachable sets. We use rr or sample to refer to one RRR. Since each RRR is a set of vertices, we refer the cardinality of an RRR as the size of RRR.

IMM and Ripples. Although the reverse reachable scheme of Borgs et al. greatly improves the efficiency of MC simulations, it could overestimate the required number of MC simulations and waste a lot of computation. The IMM algorithm proposed by Tang et al. (2015) adopts a two-phase design to algorithmically determine the sampling effort by leveraging a martingale strategy¹¹1The martingale strategy is a betting strategy in which a person doubles the bet every time they lose knowing that they must win at some point. It first samples a small number of RRR sets to calculate the achieved influence and then doubles the number of samples until the influence reaches an estimated lower bound. . This strategy greatly improves the performance so that IMM algorithm can analyze large graphs with millions of vertices.

Specifically, the two-phase design works as follows. In the Sampling phase, the algorithm will produce $\theta$ RRRs starting from random vertices in the input graph and simulating the diffusion model (IC or LT) in reverse. In the case of the IC and LT model, the task of generating a RRR set resembles a randomized breadth-first search (BFS) where only a subset of the neighbors of a vertex enter the next frontier. To estimate the required sampling effort ($\theta$) to achieve the approximation guarantee (controlled through the parameter $\epsilon$), the algorithm uses two important results. Borgs et al. showed that the fraction of RRR covered during the seed selection process is an unbiased estimator of the influence function while Tang et al. were able to prove a lower bound on the sampling effort using an estimation of the influence function. In Equation (1), Tang et al. shows the algorithm starts with a guess on $\theta$ and doubles the sampling effort at every iteration until the exit condition derived from the lower lower bound is satisfied and the final value of $\theta$ is computed. Here $n$ is the number of vertices, $k$ is the number of seeds, $\epsilon$ is the error factor, and $\theta$ is the required number of sampling. The larger $n$ and $k$ are, the larger $\theta$ will be. On the other hand, small error factor $\epsilon$ will increase $\theta$ non-linearly. This process resembles to the martingale betting strategy.

The Seed Selection algorithm is based on the greedy maximum coverage algorithm (Vazirani, 2001). The method iterates $k$ times over the RRR sets to select the vertex $v$ appearing most frequently. At every iteration $i$, the RRR sets covered by one of the $i$ seeds already selected are ignored. Minutoli et al. (2019, 2020a) devise efficient parallel schemes that perform the counting by either updating the state of the previous iteration or recount from scratch when more profitable.

Ripples (Ripples, 2022) is a state-of-the-art parallel software framework that provides fast and scalable implementation for the IM problem. According to (Minutoli et al., 2019, 2020a), its CPU version provides a speedup of 580$\times$ over the best sequential baseline using 1024 nodes, and its CPU+GPU version provides a speedup of 760$\times$ over the best sequential baseline.

Huffman coding is a classical data compression technique (Huffman, 1952). It assigns variable-length codes to encode target characters based on their relative frequency, which gets better reduction when the data has skewed distribution (Roth and Van Horn, 1993). The Huffman codes are prefix-free and are typically created as a binary tree with the encoded characters stored at the leaves. There are works that adopt Huffman coding to reduce the memory footprint. For example, Suontausta and Tian (Suontausta and Tian, 2003) applied Huffman coding for efficiently storing decision tree parameters to minimize the memory footprint. Ficara et al. (Ficara et al., 2010) adopted Huffman coding to improve counting bloom filters in terms of fast access and limited memory consumption.

Bitmap is a mapping from some domain (e.g., a range of integers) to bits (Foley et al., 1996). It not only reduces the data size but also allows efficient direct operations such as binary logic operations (Martínez-Bazan et al., 2012). Thus, many efficient bitmap compression schemes have been extensively studied in the database systems, such as BBC, WAH, EWAH, and Roaring (Wang et al., 2017). However, none of these work exploited data analytics on the Huffman or bitmap encoded data directly.

Due to the growing sizes of graphs, researchers have been studying compression techniques for graphs. For example, Randall et al. (Randall et al., 2002) tried to compress the links of web graphs by leveraging the locality of web graphs. Ligra (Shun and Blelloch, 2013) and SlimGraph (Besta et al., 2019) focused on providing frameworks that can facilitate graph analytics with compression scheme. The Ligra framework mainly utilizes the vertex degrees and graph density to select a scheme to map vertices or edges to integer-arrays or bit-vectors, while SlimGraph focuses more on utilizing statistics of local parts in the input graph to map vertices and edges to higher hierarchies. Thus, its compression kernels can preserve some critical graph properties. However, the use of these frameworks needs to re-program the graph applications with their specific syntax and semantics. This is a non-trivial effort for most of the graph algorithms and applications in the literature.

Moreover, there are a few works that investigate data analytics directly on compressed data without decompression. For example, Zhang et al. proposed an approach to perform document analytics (word count, inverted index, and sequence count) directly on compressed textual data on CPUs (Zhang et al., 2021b) and GPUs (Zhang et al., 2021a; Pan et al., 2021; Liu et al., 2022). Moreover, Zhang et al. developed a new storage engine, called CompressDB, which can support data processing for databases without decompression (Zhang et al., 2022). Furthermore, Macko et al. (Macko et al., 2015) proposed LLAMA that performs graph storage and analysis on the compressed sparse row (CSR) representation and achieves performance gain on graph benchmarks (i.e., PageRank, BFS, and triangle counting) due to in-memory execution. In addition, Mofrad et al. (Mofrad et al., 2019) proposed a compression technique specifically designed for matrix-vector operations based on the compressed sparse column (CSC) representation and leverage this compression to accelerate distributed graph benchmarks (e.g., PageRank, single source shortest path, and BFS). However, no work has been done on using compression to improve memory efficiency and accelerate a real-world diffusion-based graph application such as influence maximization.

The key component of the Monte Carlo (MC) diffusion process (Tang et al., 2015; Minutoli et al., 2019) is a probabilistic Breadth First Search. The intermediate collections of vertices that may cause activation of the BFS root are saved in the form of RRR sets. We benchmark six real-world large graphs used in the literature and study the memory usages for the MC diffusion process and generating/saving the intermediate RRR sets. As shown in Figure 1, the space attributed to storing of the intermediate RRR sets dominate, consuming between 87% to 99% of the memory footprint.

To further understand the memory footprint usage patterns in Figure 1, we studied the distribution of the RRR sets by their sizes. Figure 2 shows the sizes of RRR sets, i.e., their vertex counts, for different graphs. The figure illustrates the shapes of RRR sets’ distribution based on their sizes: (1) In the top row, the two graphs (i.e., DBLP and YouTube) have long tails, and their shapes correspond to a power-law distribution. (2) In the middle row, the shapes of the two graphs (i.e., Skitters and Orkut) show close to a linear decay without any long tails. (3) In the bottom row, the shapes of the two graphs (i.e., Pokec and LiveJournal) show a very uniform two-step (flat-head) distribution. The overall memory footprint of the RRR sets corresponds to the area under the curves in the respective plots. This well explains why the last two graphs spends much more memory in saving the RRR sets. Further, these observations raise a critical question to the design of optimizations for memory footprint reduction: What are the characteristics of the various RRR sets of the above three different categories? (see Section 3.2)

The shapes in Figure 2 roughly categorize the RRR sets’ distributions into two types: the first four graphs (i.e., DBLP, YouTube, Skitter and Orkut) have skew-distributed distribution while the last two graphs (i.e., Pokec and LiveJournal) have a uniform two-step distribution. We use two quantities, i.e., skewness and density, to characterize the shapes of RRR sets’ distributions.

The skewness score ($\mathcal{S}$) measures how distributions are asymmetric about their mean: When $\mathcal{S}<0$, the distribution has a longer left tail; When $\mathcal{S}>0$, the distribution’s right tail is longer. Equation (2) shows the calculation of skewness score for the size of RRR sets $X=(X_{1},....,X_{\theta})$, where each $X_{i}$ is a sample of the number of visited vertices starting from a seed vertex in a MC simulation. Thus all $X_{i}$’s are at least $1$ because there will be no empty resultant RRR set (it will at least contain the starting seed vertex); $\theta$ is the total number of sampled RRR sets; $\bar{X}$ is the average RRR size, and $s$ is their standard deviation.²²2The standard deviation $s\neq 0$ because $X_{i}$’s are the sizes of RRR sets, which are not all equal. Given the elements in Equation (2) are all non-zero, the skewness score and density will not overflow by dividing zeros for real-world social networks. In Table 1, we show the skewness scores of the first four graphs (i.e., DBLP, YouTube, Skitter and Orkut) are all positive, which help to distinguish them from the last two graphs (i.e., Pokec and LiveJournal) whose skewness scores are negative.

The density ($\mathcal{D}$) measures the proportion of non-zero elements that are required to represent the sampled RRR sets if they are stored in matrix format. Although density is not correspond to the shape of distributions, it can help to select the data structure to store RRR sets. It is sufficient to use 32-bit unsigned integers to represent each vertex in the sampled RRR sets in this study. Thus we can use $3.12\%$ to be the threshold. When the density $\mathcal{D}\leq\frac{1}{32}=3.12\%$, it is more efficient to use a sparse representation of RRR sets, i.e., to explicitly store each vertex. On the other hand, when the density $\mathcal{D}>3.12\%$, it is more efficient to use a dense bitmap coding to store RRR sets. Equation (2) shows the calculation of density for the size of RRR sets $X$ (Similarly, $\theta$ is the total number of RRR sets; $X_{i}$ is the number of vertices in each RRR set; $n$ is the number of vertices in graph $G$). In Table 1, we show the last two graphs (Pokec and LiveJournal) have much higher densities ($\geq 50\%$) than the rest four graphs (DBLP, YouTube, Skitter, Orkut).

The different shapes of RRR sets’ distribution also characterize the influence of individual vertices which can help us to find data localities. To verify this intuition, we compare the selected seeds of each graph by varying the random starting vertex for MC samplings. Table 2 shows the Rank-Biased Overlap (RBO) scores (Webber et al., 2010) ($1.0$ means highly overlapping rank positions and $0.0$ means no overlapping). For the graphs with skew-distributions, the same top-1 influential vertex is consistently selected from different random starting vertices. However, for the graphs with flat-head distributions, many vertices achieve the maximum influence, so the RBO score is zero.

From the first four graphs (DBLP, YouTube, Skitter, Orkut), the sampled RRR sets are close to power-law or linear decay distributions. Under this situation, there are only a small number vertices that influence a lot of other vertices; while many vertices only influence one or two other vertices. On the other hand, for the last two graphs (Pokec and LiveJournal), the flat-head shows that many vertices can influence many other vertices. In other words, many vertices are equally influential, which makes the flat-head distributed RRRs lack the data locality in the skew-distributed RRR sets. Thus, they need different optimization strategies to reduce their memory footprints.

In summary, the memory required to store intermediate RRR sets is one to two order of magnitude larger than the memory footprint of the diffusion process. We observe two types of graphs with regards to their different behaviors in diffusion process and hence the data localities. Considering the observed skewed-distributed RRR sets, it is promising to leverage variable-length encoding to reduce the memory footprint of RRR sets for those graphs. On the ther hand, using bitmap coding can reduce the memory footprints for the flat-head distributed RRR sets because of their high density of non-zero elements.

The overview of the “compress-to-compute” workflow is shown in Figure 3. It contains three main phases: warm-up, sample-and-encode, and decode-and-select. Specifically, we first characterize the skewness score and the density from the distribution of the RRR sets in the warm-up phase. Based on that, we then choose an efficient compression technique to encode the intermediate RRR sets and query on the encoded RRR sets.

More specifically, we introduce a block-based sampling-and-encoding approach and an integrated selection approach. To achieve a better scalability on multicore architectures, we design two schemes - (1) Huffmax: We use Huffman coding to encode skew-distributed RRR sets, and leverage the data locality knowledge acquired from the warm-up phase to accomplish the selection without fully decoding the compressed RRR sets. (2) Bitmax: We use bitmap coding to encode flat-head distributed RRR sets. The corresponding selection can be accomplished directly on the encoded data by efficient bit operations. These two schemes are complementary to each other so that the proposed workflow does not rely on the data locality of skewed distributions. It can work with RRR sets that have non-skewed distributed sizes such as normal distributions or uniform distributions (their skewness scores are zero). In Figure 2, we show the flat-headed distributions of Pokec and Livejournal are two special cases that are close to uniform distributions. Furthermore, we use a heuristic approach to optimize the parallel reduction operation to find the vertex with the highest frequency. We use Table 3 to list the notation used in the paper.

Next, we describe the proposed two encoding methods.

Lastly, we describe two selection approaches for efficiently querying in Huffmax  and Bitmax, respectively. The selection iterates two core computations for $k$ times to select the most influential seeds. The first computation is to locate and remove the RRR sets which contain the current most frequent vertex $u^{*}$. The second computation is to reconstruct the frequency table $\hat{h}$ after every RRR set that contains the previous $u^{*}$ is removed. After that, a new $u^{*}$ is selected based on the updated frequency table $\hat{h}$.

In this section, we evaluate the performance of our proposed HBMax. We first evaluate on the local workstation to compare the memory usage and time-to-solution of HBMax with Ripples because it has larger memory (376 GBs). The experiment results are measured with 16 cores. For a fair comparison, we repeat the experiments on each tested graph for five times and take the averages to compare memory usage and time-to-solution of HBMax with Ripples. We show the averaged experiment results in Table 6, 7 and 8. We then use Bridges-2 to study the scalability of HBMax because its compute node has more CPU cores (i.e., 64 cores per CPU). We show the scalability results in Figure 5 and 6.

We now present the strong scaling result of HBMax with up to 64 CPU cores, as shown in Figure 6. We break down HBMax into four operations: (1) sampling the RRR sets with MC diffusion process, (2) encoding the RRR sets with Huffmax, (3) encoding the RRR sets with Bitmax method, and (4) selecting the most influential seeds. We analyze their impacts to the overall scalability separately.

First, Figure 5 shows that the parallel sampling step is the dominant part on all evaluated graphs (takes an average of $83.30\%$ of the total time) and has a strong scalability due to the independent nature of each sampling operation; thus, the scalability of HBMax is affected by the encoding step and selection step.

Second, for Huffmax, we note that building the Huffman codebook does not affect the overall scalability. This is because we only build up the codebook once during the warm-up phase, which takes an average of $2.19\%$ of the total time. Also, the Huffman encoding takes an average of $8.77\%$ of the total time as it only includes a small number of synchronizations to build frequency table. Thus, Huffmax’s encoding shows a strong scalability up to 64 cores.

Third, for all the large graphs (Pokec, LiveJournal, Arabic, Twitter7), the encoding of Bitmax shows a strong scalability up to 64 cores due to its independent nature. The bitmap encoding takes an average of $8.19\%$ of the total time.

Fourth, the selection operation in both Huffmax and Bitmax does not affect the overall scalability. First, our optimized reduction with parallel merge plays an important role, as shown in Figure 4, to minimize the synchronization between threads. Without this optimization, the selection operation would become the bottleneck due to the original non-scalable parallel reduction. Second, both the selection operations of Huffmax and Bitmax are parallelized considering the NUMA effect so that each thread can maintain its local frequency table by only accessing the locally encoded data. Note that selection of Bitmax scales better than Huffmax because it does not have a decoding step.

Finally, we note that for the two small graphs DBLP and YouTube, Ripples does not scale well. This is because of the highly imbalanced workload between vertices in the seed-selection step, which can be also reflected in its high skewness score. While HBMax achieves better scalability (can scale up to 32 cores) thanks to the proposed parallel merge and the consideration of NUMA effect. The performance starts to degrade from 64 cores due to the relatively small workload. For the other six larger graphs, HBMax can scale up to 64 cores and achieve an average speedup of 12.98$\times$ on 64 cores.