1 Introduction↩︎

Belief Propagation (BP) is a widely used approximate inference algorithm for probabilistic graphical models (PGMs) [1]. In these models, random variables are represented as nodes, and their dependencies are captured by the graph structure. The primary task of inference is to estimate the posterior marginal distributions of these variables given observed evidence. Exact inference is often computationally intractable, making approximate algorithms a more practical alternative. BP is a prominent approximate inference algorithm [2] that operates by iteratively passing messages between neighboring nodes until convergence or for a specified number of iterations. Due to its simplicity and efficiency, BP has been successfully applied in diverse fields including computer vision [3], channel coding [4], protein folding [5], and, notably, program analysis [6].

In recent years, the combination of program analysis and probabilistic reasoning has enabled the application of BP to program analysis [6]–[15]. For instance, BINGO [6] converts the Datalog program analysis into a directed acyclic graph, transforms the graph into a Bayesian network by attaching probabilities to logic rules, and applies BP to infer the likelihood of bug reports. SmartFL [9], [14] models program semantics as a factor graph by utilizing information from static analysis results and dynamic execution traces, and employs BP to infer the probabilities of program elements being faulty.

Program analysis on large-scale software often yields large PGMs, leading to high computational costs for BP. Graphics Processing Units (GPUs), optimized for massive parallel operations, are a natural choice for reducing the inference cost of BP. This potential has motivated prior GPU-based BP approaches [16]–[21].

However, existing GPU-based BP approaches suffer from two key limitations when applied to program analysis. First, they lack support for flexible message update strategies required in program analysis. In the existing BP literature, researchers have devoted considerable attention to update strategies as a means of improving BP performance. The most straightforward scheme is synchronous BP, where all messages are updated simultaneously [2]. In contrast, asynchronous schemes introduce a degree of sequential updating, such as tree-based updates [22] and greedy message selection [18], [23]–[25]. While synchronous BP exposes high levels of parallelism and can approximate the correct posteriors well when it converges, asynchronous BP has been shown, both empirically and theoretically, to converge more reliably when synchronous BP fails to converge [23]. Program analysis applications often require different update strategies depending on usage scenarios, yet existing GPU-based implementations support only limited update strategies. For example, widely used tools such as PGMax [21] support only synchronous BP, restricting their applicability in program analysis.

Second, BP with local structures [26] has not yet been explored in GPU-based implementations, although it is theoretically efficient. In program analysis, probabilistic constraints often exhibit inherent local structures, where many variable assignments share the same probability value. A typical example is the probabilistic Horn constraints in the form of \(p:x_1\land\cdots\land x_n\rightarrow y\), which is ubiquitous in modeling program semantics. While traditional BP requires an exponential number of operations (\(O(n\cdot2^n)\)) to marginalize over \(n\) variables in such a constraint, BP with local structures can reduce this complexity to linear time (\(O(n)\)) by grouping assignments with identical probabilities and extracting common terms. However, existing GPU-based approaches typically rely on dense tabular representations that fail to capture these computational redundancies. Consequently, due to these limitations in generality and algorithmic efficiency, CPU-based approaches [26], [27] still remain the preferred choice for program analysis applications.

We present FastLBP, a GPU-accelerated BP framework designed for program analysis. First, to accommodate diverse update strategies required in program analysis, we propose a unified representation for user-specified update strategies. We further develop a dependency analysis algorithm that identifies groups of messages that can be updated safely in parallel. By partitioning edges into ordered groups, FastLBP preserves the semantics of user-specified update strategies while enabling parallel execution. Second, to address the inefficiency caused by traditional BP, we incorporate BP with local structures into our GPU design. We organize messages and factor configurations in a compact layout to enable efficient memory access. We further design a GPU execution scheme that maps each message computation to a thread, where each thread evaluates a sequence of sub-messages. This design reduces the computational complexity by efficiently implementing BP with local structures.

We conduct extensive experiments to evaluate the efficiency and accuracy of our algorithm. We implement FastLBP in a prototype system and compare it with state-of-the-art approaches [21], [26]. We evaluate FastLBP using two representative approaches: SmartFL, a probabilistic fault localization approach, and BINGO, a Bayesian program analysis approach, to demonstrate the effectiveness of our approach. Experimental results show that: (1) compared to the state-of-the-art CPU-based approach [26], our approach achieves average speedups of \(17.42\times\) and \(2.82\times\) on SmartFL and BINGO, respectively, while maintaining the accuracy. Compared to the state-of-the-art GPU-based approach [21], our approach achieves an average speedup of \(6.14\times\) on SmartFL. (2) FastLBP supports update strategies used in BINGO that are not supported by existing GPU-based approaches, highlighting its advantage in enabling flexible update strategies.

In summary, our contributions are as follows:

• We present FastLBP, a GPU-accelerated BP framework for program analysis, enabling flexible update strategies and efficient inference on large-scale PGMs.

• We propose a unified representation for user-specified update strategies and a dependency analysis algorithm for effective parallelization.

• We design a GPU implementation of BP with local structures by assigning each GPU thread to compute a single message and organizing the memory in a compact way.

• We implement a high-performance prototype and conduct extensive experiments to demonstrate the flexibility and efficiency of our approach.

2 Motivating Example↩︎

In this section, we illustrate the key ideas of our approach using an example generated by BINGO [6], which performs BP on PGMs derived from analysis rules.

2.1 Support for Update Strategies↩︎

BP is a message-passing algorithm widely used for performing inference on factor graphs. 1 shows a subgraph of a factor graph generated by BINGO. A factor graph is a bipartite graph consisting of two types of nodes, variables and factors. For example, the variables are \(V = \{v_1, v_2, v_3, v_4, v_5\}\) and the factors are \(A = \{a_1, a_2, a_3, a_4, a_5\}\) in 1. BP operates by iteratively passing messages along edges between variables and factors. A factor-to-variable message is computed depending on incoming variable-to-factor messages, while a variable-to-factor message is computed based on incoming factor-to-variable messages. For each message update, the incoming messages are those from neighboring nodes except the target node. For example, message from \(a_4\) to \(v_1\), denoted as \(\mu_{a_4\to v_1}\), depends on messages from \(v_2\) and \(v_4\) to \(a_4\), denoted as \(\mu_{v_2\to a_4}\) and \(\mu_{v_4\to a_4}\). Moreover, \(\mu_{v_2\to a_4}\) depends on \(\mu_{a_2\to v_2}\), and \(\mu_{v_4\to a_4}\) depends on \(\mu_{a_5\to v_4}\).

Update strategies play a critical role in the performance and convergence behavior of BP. They define how message updates are scheduled. For instance, synchronous BP requires all messages to be updated in parallel, whereas asynchronous BP enforces a degree of sequentialism in message updates. Although synchronous BP is well suited for parallel execution, asynchronous BP has been shown to converge more reliably [23], which is important for many program analysis applications. For example, BINGO ranks alarms according to the inferred probabilities. Asynchronous BP often produces more accurate probabilities than synchronous BP, thereby reducing the effort required for further inspection of bug reports. At the same time, synchronous BP remains useful because it runs faster and produces good results when it converges.

However, existing approaches are unable to support flexible update strategies. One line of GPU-based work designs update strategies for specific domains, such as computer vision [17] and channel code decoding [28]. PGMax [21] is a Python package for running BP efficiently on both CPUs and GPUs and supports arbitrary graphs, but it only supports synchronous BP. Another line of work studies the convergence behavior of different update strategies in sequential and parallel scenarios [18], [23]–[25], but lacks a general mechanism to capture user-specified dependencies among message updates. Our work differs from with these works by providing a systematic way to capture updating dependencies, thereby enabling a wider range of update strategies for GPU parallel computation.

We illustrate a complex update strategy that is not supported by existing GPU-based approaches. To balance convergence stability and computational efficiency, BINGO adopts a round-robin asynchronous update strategy, where messages are updated sequentially in a predefined order during each iteration. When updating a message, the most recently updated neighboring messages are used whenever available. Although this sequential schedule improves convergence behavior, it introduces strict data dependencies between message updates. These dependencies arise from both the user-specified order and the message update rules of BP. As a result, updating multiple messages simultaneously may violate the semantics of the asynchronous schedule.

Consider the factor graph in 1. Suppose the factor-to-variable messages are updated in the following order during an iteration: \[\begin{align} \mu_{a_1\to v_1},\; \mu_{a_2\to v_2},\; \mu_{a_3\to v_1},\; \mu_{a_3\to v_3},\; \mu_{a_4\to v_1},\\ \mu_{a_4\to v_2},\; \mu_{a_4\to v_4},\; \mu_{a_5\to v_3},\; \mu_{a_5\to v_4},\; \mu_{a_5\to v_5}. \end{align}\] Let \(\mu_{a\to v}^{i}\) denote the message from factor \(a\) to variable \(v\) at iteration \(i\). In this example, \(\mu_{a_4\to v_1}^{i + 1}\) depends on \(\mu_{a_2\to v_2}^{i + 1}\) because the latter is updated earlier in the schedule, and on \(\mu_{a_5\to v_4}^i\) because that message has not yet been updated. Such ordering constraints create dependencies that limit naive parallel execution. However, not all messages are mutually dependent. Two messages can be safely updated in parallel if neither requires the updated value of the other. In this example, \(\mu_{a_2\to v_2}^{i+1}\) is independent of \(\mu_{a_1\to v_1}^{i+1}\), and \(\mu_{a_3\to v_1}^{i+1}\) is independent of both. Therefore, these messages can be updated concurrently without violating the asynchronous semantics. However, existing GPU-based approaches fail to capture such dependencies.

To enable parallel execution while preserving the semantics of different user-specified update strategies, our approach first provides a unified representation of update strategies, namely a preordered set \((E, \lesssim)\) over the edge set \(E\) of the factor graph. Intuitively, specifying \(e_1\lesssim e_2\) means that, in a single iteration, the update of the message on edge \(e_1\) must happen no later than that on edge \(e_2\). We then perform dependency analysis that captures both the update strategy constraints and BP computation dependencies. Based on this analysis, messages can be partitioned into ordered groups that can be updated in parallel.

In this example, our approach represents the update strategy as a total order induced by the preorder according to the updating sequence: \[\begin{align} &(a_1, v_1)\prec (a_2, v_2)\prec (a_3, v_1)\prec (a_3, v_3)\prec (a_4, v_1)\prec\\ &(a_4, v_2)\prec (a_4, v_4)\prec (a_5, v_3)\prec (a_5, v_4)\prec (a_5, v_5). \end{align}\] Our approach then performs a topological traversal on the Hasse diagram of this order and partitions the messages into groups. Specifically, for \((a, v)\prec (a', v')\), our approach ensures that if \(\mu_{a'\to v'}\) depends on \(\mu_{a\to v}\) in the BP computation, then \(\mu_{a\to v}\) must be computed strictly earlier than \(\mu_{a'\to v'}\); otherwise, the two updates can be performed in parallel. First, we can update \(\mu_{a_1\to v_1}\), \(\mu_{a_2\to v_2}\), and \(\mu_{a_3\to v_1}\) in parallel, as we discussed earlier in this section. Next, we place \(\mu_{a_3\to v_3}\) into a new group because it depends on \(\mu_{a_1\to v_1}\). We then add \(\mu_{a_4\to v_1}\), \(\mu_{a_4\to v_2}\), and \(\mu_{a_4\to v_4}\) to the same group one by one, because these four messages are mutually independent. But we place \(\mu_{a_5\to v_3}\) into another group because it depends on \(\mu_{a_4\to v_4}\). Finally, \(\mu_{a_5\to v_4}\) and \(\mu_{a_5\to v_5}\) can be placed in the same group as \(\mu_{a_5\to v_3}\). As a result, the messages can be partitioned into three groups: \[\begin{gather} \{\mu_{a_1\to v_1},\; \mu_{a_2\to v_2},\; \mu_{a_3\to v_1}\}, \\ \{\mu_{a_3\to v_3},\; \mu_{a_4\to v_1},\; \mu_{a_4\to v_2},\; \mu_{a_4\to v_4}\}, \\ \{\mu_{a_5\to v_3},\; \mu_{a_5\to v_4},\; \mu_{a_5\to v_5}\}. \end{gather}\] Updating each group in parallel, while processing the groups sequentially, preserves the semantics of the original asynchronous schedule and exposes substantial parallelism.

2.2 Implementing BP with Local Structures↩︎

Many probabilistic constraints arising in program analysis exhibit local structures [26]. Local structures occur when multiple assignments of variables share the same probability value in a factor. Instead of representing the factor as a full table over all assignments, these repeated patterns can be grouped and computed more efficiently. A common example in program analysis is the probabilistic Horn constraint \(p : x_1 \land x_2 \land \cdots \land x_n \rightarrow y\). Such a constraint states that if all premises hold, the conclusion holds with probability \(p\). When the constraint is encoded as a factor in a factor graph, many assignments of the variables share identical probabilities. Traditional BP evaluates factor-to-variable messages by enumerating all assignments of neighboring variables, resulting in an exponential complexity of \(O(n\cdot2^n)\). BP with local structures exploits the repeated patterns in the factor to group assignments with identical probabilities, reducing the computation to linear time \(O(n)\).

To leverage this advantage on GPUs, we extend the standard message-level parallelization of BP to support BP with local structures. We assign each GPU thread to compute a single factor-to-variable message. By providing each thread with the corresponding edge, factor, and variable identifiers, we can efficiently access the required messages and factor configurations. 2 illustrates the computation process using the message group: \[\{\mu_{a_3\to v_3},\; \mu_{a_4\to v_1},\; \mu_{a_4\to v_2},\; \mu_{a_4\to v_4}\}\text{.}\] In BP with local structures, a factor-to-variable message is computed as the sum of several sub-messages. Each sub-message corresponds to a group of assignments that share the same probability value. It is either zero or computed as a constant multiplied by the product of incoming variable-to-factor messages. For instance, for the constraints generated by BINGO, messages are defined over binary variables (domain \(\{0, 1\}\)). Each entry is computed as the sum of five sub-messages, resulting in ten sub-messages in total.

Another problem is how to represent factor graphs efficiently, as those generated by program analysis applications are typically large and sparse. To achieve memory-efficient storage and efficient access, we adopt the Compressed Sparse Row (CSR) [29] format to represent messages. In CSR, a sparse matrix is stored as a flat row-major array. To compute a factor-to-variable message \(\mu_{a\to v}\), we require all incoming variable-to-factor messages to factor \(a\). Therefore, the variable-to-factor message matrix is organized such that rows are indexed by factor identifiers, grouping messages associated with the same factor contiguously. Conversely, to compute a variable-to-factor message \(\mu_{v\to a}\), we require all incoming factor-to-variable messages to variable \(v\), and thus the factor-to-variable message matrix is indexed by variable identifiers. For example, the first row (corresponding to \(v_1\)) contains \(\{\mu_{a_1\to v_1}, \mu_{a_3\to v_1}, \mu_{a_4\to v_1}\}\), and the second row (corresponding to \(v_2\)) contains \(\{\mu_{a_2\to v_2}, \mu_{a_4\to v_2}\}\). Additionally, factor configurations are stored in a similar manner as variable-to-factor messages. Unlike standard BP, which sums over all factor values, BP with local structures computes sub-messages based on an indicator specifying whether a variable is fixed or treated as a wildcard. Instead of using a dense tabular representation, we store these configurations in a flat array. 3 illustrates the memory layout of the factor graph in 1 and how \(\mu_{a\to v}\) accesses the required messages and factor configurations.

3 Preliminaries↩︎

In this section, we introduce the background of our approach, including formal definitions of factor graphs and belief propagation, as well as a brief overview of GPU programming.

3.1 Factor Graph↩︎

A factor graph is a special type of PGM. Common PGMs such as Bayesian networks and Markov random fields can be converted into factor graphs, enabling a unified approach to inference. A factor graph is a bipartite graph representing the factorization of joint probability distribution. Given a set of random variables \(\mathbf{X} = \{X_1, X_2, \dots, X_n\}\), suppose their joint probability distribution can be factorized into a product of functions: \[P(\mathbf{X}) = \frac{1}{\mathcal{Z}} \prod_{a \in A} f_a(\mathbf{X}_a),\] where \(\mathbf{X}_a \subseteq \mathbf{X}\) is a subset of \(\mathbf{X}\), and \(\mathcal{Z}\) is the normalization constant. The corresponding factor graph \(G = (V, A, E)\) consists of variable nodes \(V = \{v_1, v_2, \dots, v_n\}\), factor nodes \(A = \{a_1, a_2, \dots, a_m\}\), and edges \(E \subseteq V \times A\). Each node \(v_i\) corresponds to the random variable \(X_i\) (or \(X_{v_i}\) for notational convenience) and each node \(a_j\) corresponds to the function \(f_{a_j}\). An edge \((v_i, a_j) \in E\) exists if and only if \(X_i \in \mathbf{X}_{a_j}\).

3.2 Belief Propagation↩︎

Belief Propagation (BP) is a message-passing algorithm for performing inference on factor graphs. BP iteratively exchanges messages between variables and factors along the edges of a factor graph and computes marginal distributions of variables with these messages. For each edge \((v, a) \in E\), two types of messages are defined, a message \(\mu_{v\to a}\) from \(v\) to \(a\), and a message \(\mu_{a \to v}\) from \(a\) to \(v\). Both types of messages, \(\mu_{v\to a}, \mu_{a \to v}: \text{Range}(X_v)\to\mathbb{R}\), are real-valued functions whose domain is the set of values that can be taken by the random variable associated with \(v\).

The message sent from a variable \(v\) to a neighboring factor \(a\) is defined as the product of the incoming factor-to-variable messages: \[\label{eq:message-vtoa} \mu_{v \to a}(x_v) = \prod_{a^* \in N_V(v) \setminus \{a\}} \mu_{a^* \to v}(x_v),\tag{1}\] where \(N_V(v)\) denotes the set of neighboring factors of \(v\).

The message sent from a factor \(a\) to a neighboring variable \(v\) is defined as the product of the factor function with the incoming variable-to-factor messages: \[\label{eq:message-atov} \mu_{a \to v}(x_v) = \sum_{\mathbf{x}_a \setminus x_v} f_a(\mathbf{x}_a) \prod_{v^* \in N_A(a) \setminus \{v\}} \mu_{v^* \to a}(x_{v^*}),\tag{2}\] where \(N_A(a)\) denotes the set of neighboring variables of \(a\), and \(\sum_{\mathbf{x}_a \setminus x_v}\) denotes summation over all possible values of \(\mathbf{x}_a\) while keeping \(x_v\) fixed.

After convergence or a certain number of iterations, the marginal probability of a variable \(v\) is computed as \[\label{eq:marginal} P(X_v = x_v) = \frac{\prod_{a \in N_V(v)} \mu_{a \to v}(x_v)}{\sum_{x_v} \prod_{a \in N_V(v)} \mu_{a \to v}(x_v)}.\tag{3}\]

If the factor graph is a tree, BP yields exact marginal distributions after a finite number of iterations. For graphs with cycles, BP is not guaranteed to converge but often provides good approximate solutions in practice.

3.3 BP with Local Structures↩︎

The main computational bottleneck of BP lies in the factor-to-variable message \(\mu_{a\to v}\) in 2 , which requires enumeration over an exponential number of factor values. Local structures in a constraint can be exploited to simplify this computation by grouping assignments that share the same factor value. The computation of \(\mu_{a\to v}\) is then reduced to the summation of multiple sub-messages based on the local structure of the factor function. Formally, \[\label{eq:wu-atov} \mu_{a \rightarrow v}(x_v) = \sum_{k=1}^{l} \sum_{\mathbf{x}_a \backslash x_v} t_k(\mathbf{x}_a) \prod_{v^* \in N_A(a) \setminus \{v\}} \mu_{v^* \to a}(x_{v^*})\text{,}\tag{4}\] where each \(t_k\) specifies a weight together with indicators indicating whether a variable is constrained to a fixed value or treated as a wildcard. According to 4 , each message is computed using 4, where \(v\) is the target variable index, \(x_v\) is the target value, \(p[l]\) denotes the weights associated with the factor, \(b = |N_A(a)|\) is the number of variables adjacent to \(a\), and \(l\) is the number of sub-messages. The indicator arrays are \(indicator[b][l]\), and \(\mu[b][r]\) stores messages from \(N_A(a)\) to \(a\). As shown in 4, the complexity is reduced from \(O(b \cdot r^b)\) to \(O(l \cdot b \cdot r)\), where \(b\) can be large.

3.4 GPU Programming↩︎

We introduce GPU programming using Compute Unified Device Architecture (CUDA) as a representative example. Modern GPUs consist of multiple streaming multiprocessors (SMs), each of which contains many compute units for massive parallel execution. The functions executed on GPU are referred to as kernels. A kernel is launched to run multiple threads in parallel. These threads are grouped into thread blocks, while these blocks together form a grid. Each thread block is scheduled to run on a single SM. Within a thread block, threads are divided into groups of 32 threads known as warps. Threads within a warp execute instructions in a Single-Instruction Multiple-Threads (SIMT) paradigm.

4 FastLBP Framework↩︎

In this section, we first decribe the overall workflow of FastLBP (4.1). We then present two key designs of FastLBP: (1) the update strategy representation with the dependency analysis algorithm (4.2), and (2) the implementation of BP with local structures on GPUs (4.3).

4.1 Overall Workflow↩︎

5 shows the heterogeneous workflow of FastLBP. Given a factor graph derived from probabilistic program analysis and a user-defined update strategy, FastLBP ultimately produces marginal probabilities of variables as results. The workflow consists of the following four major components.

Initialization. FastLBP first takes as input a factor graph, typically generated by probabilistic program analysis frameworks such as BINGO or SmartFL. Based on the graph structure and factor functions, FastLBP allocates GPU memory for messages using a CSR-based layout. This representation enables compact storage and efficient access for sparse message structures, avoiding the prohibitive memory cost of dense adjacency matrices on GPUs.

Dependency Analysis. FastLBP requires the user to specify an update strategy. FastLBP represents the update strategy as a preordered set \((E, \lesssim)\) over graph edges. This abstraction uniformly captures different update strategies. It then analyzes the dependencies implied by this preordered set and BP computation process. The result is a partitioning of edges into ordered groups, each representing a set of factor-to-variable messages that can be updated simultaneously.

Kernel Launch. Based on the dependency analysis results, FastLBP organizes message updates into batches corresponding to the computed edge groups. For each group, a GPU kernel is launched, where each thread is assigned to compute a single message. The required factor, variable, and edge identifiers are passed to each thread to enable efficient data access. This batching strategy reduces kernel launch overhead and ensures that each kernel exposes sufficient parallelism for efficient GPU execution.

Message Passing. Finally, FastLBP performs iterative message passing on the GPU. In each iteration, kernels are launched following the dependency order, and messages are updated accordingly using the message-level parallel implementation. The process continues until convergence or a user-specified number of iterations, after which marginal probabilities are computed and returned to the analysis framework.

4.2 Representation and Dependency Analysis↩︎

In this section, we first briefly introduce the representation of update strategies, then formalize its definition and semantics, and finally present the dependency analysis algorithm.

4.2.1 Update Strategies↩︎

An update strategy in FastLBP is represented as a preordered set \((E, \lesssim)\) over graph edges. 1 describes several typical update strategies. Users can specify the strategy name, and FastLBP encodes it together with the edge set into a unified intermediate representation, i.e., a preordered set. Additionally, users may define a custom preordered set directly as the update strategy.

4.2.2 Representation Formalization↩︎

To perform dependency analysis in a uniform manner, we introduce a unified representation of update dependencies. Our key insight is that update strategies can be uniformly characterized as preorders over the edge set of a factor graph. This abstraction allows different scheduling semantics to be handled within a single formal framework.

We first characterize the dependencies induced by different update strategies. The problem lies in whether a message depends on messages from the current iteration or the previous iteration. We focus on factor-to-variable updates, since they are directly used to compute marginal probabilities. We assume that each factor-to-variable message is updated once per iteration. When computing a factor-to-variable message, all incoming variable-to-factor messages must first be computed. These variable-to-factor messages, in turn, depend on the incoming factor-to-variable messages. By treating the computation of variable-to-factor messages as an intermediate step, the essential dependency lies in whether the required factor-to-variable messages come from the current or the previous iteration. Formally, by substituting Equation (1 ) into Equation (2 ) and introducing superscripts to indicate iteration indices, we obtain: \[\label{eq:message-dependence} \begin{align} &\mu_{a\to v}^{i+1}(x_v) = \\ &\sum_{\mathbf{x}_a\setminus x_v} f_a(\mathbf{x}_a) \prod_{v^* \in N_A(a)\setminus \{v\}} \prod_{a^* \in N_V(v^*)\setminus \{a\}} \mu_{a^*\to v^*}^{i+\delta((a, v), (a^*, v^*))}(x_{v^*})\text{,} \end{align}\tag{5}\] where \(\delta: E \times E \to \{0, 1\}\) indicates whether a dependency comes from the current iteration (\(\delta=1\)) or the previous iteration (\(\delta=0\)).

This expression shows that the update of a message on edge \((a,v)\) depends only on a specific subset of edges determined by the graph structure. To make this explicit, we define the edge dependency neighborhood.

The key to defining a legitimate update strategy is specifying a reasonable \(\delta\). We now give the definition of intermediate representation used in FastLBP for update strategies, along with the corresponding definition of \(\delta\).

Intuitively, specifying that \(e_1 \lesssim e_2\) in the preorder means that, during one iteration, the message on edge \(e_2\) cannot be updated before the message on edge \(e_1\). Specifically, \(e_1 \prec e_2\) means the message on \(e_1\) must be updated before the message on \(e_2\) if \(e_2\) depends on \(e_1\) (\(\delta(e_2, e_1) = 1\)).

The following Theorem 1 shows that the function \(\delta\) defined in Definition 2 is legitimate.

Theorem 1 shows that the update strategy defined by Definition 2 contains no cyclic dependencies. That is, within a single iteration, starting from any message, although it may depend on other messages computed in the current iteration, and those messages in turn may depend on yet others from the same iteration, this chain of dependencies must eventually terminate at messages whose values are taken from the previous iteration. This termination is guaranteed by the finiteness of the edge set \(E\). Therefore, the current iteration’s results can always be computed starting from the messages of the previous iteration, which confirms the well-definedness of Definition 2.

Definition 2 captures a sufficiently broad subset of valid update strategies. For example, PARALL can be represented by a preorder that all edges are in the same equivalence class, SEQFIX corresponds to a total order over all edges, and TOPO is expressed as a topological order.

4.2.3 Dependency Analysis Algorithm.↩︎

We now present the dependency analysis algorithm based on the intermediate representation of the update strategy.

In 6, \(\operatorname{TopologicalSort}\) on line [line:dependence-toposort] refers to a topological sort performed on the Hasse diagram of the partial order on the set of equivalence class. Edges in the same equivalence class can be in any order. After the topological sort, the algorithm greedily groups the edges into batches. When processing the edge \(e_i\), if it depends on any edge already in \(NextSet\), then the current \(NextSet\) is output as a batch, and a new batch is started. Finally, the algorithm outputs a sequence of edge sets \(S\). Updating factor-to-variable messages according to this sequence yields the expected result under the specified update strategy. The overall time complexity of the dependency analysis is \(O(|E|^2)\) in the worst case, as the algorithm performs a quadratic scan in the number of edges. Since this analysis is executed only once prior to message iterations, its cost is negligible compared to the overall inference procedure.

Theorem 2 establishes the correctness of the algorithm, guaranteeing that all message dependencies for a given batch are correctly sourced, either from previously computed batches in the current iteration or from the previous iteration, as required.

4.3 Implementation of BP with Local Structures↩︎

To leverage the computational benefits of BP with local structures on GPUs, we develop a message-level parallel implementation that maps batches of message updates to GPU execution.

4.3.1 Message-Level Parallelization↩︎

We adopt a message-level parallelization strategy, where each GPU thread is responsible for computing a single message on an edge. For each edge \((a, v) \in E\), a thread computes either a variable-to-factor message \(\mu_{v\to a}\) or a factor-to-variable message \(\mu_{a\to v}\). To compute a batch of factor-to-variable messages obtained from dependency analysis, we organize the computation into two phases. This separation avoids read-after-write dependencies and ensures that all messages can be computed fully in parallel based on the precomputed messages. In the first phase, we compute all required variable-to-factor messages \(\mu_{v\to a}\) in a kernel launch. According to 1 , each \(\mu_{v\to a}\) is computed as the product of incoming factor-to-variable messages and can be computed independently. We assign one thread per edge to compute these messages in parallel. In the second phase, we compute the factor-to-variable messages \(\mu_{a\to v}\) using BP with local structures. As shown in 4, each message is evaluated as the sum of multiple sub-messages. Each thread iterates over the sub-messages, checks the indicator conditions, and accumulates the results accordingly. This design exposes massive parallelism across edges while maintaining good data locality within each thread.

4.3.2 Memory Organization↩︎

Efficient memory organization is critical for achieving high performance on GPUs, especially for large-scale factor graphs generated by program analysis. In FastLBP, we adopt a memory layout that ensures both compact storage and efficient memory access. We store all messages in global GPU memory using a CSR format. Specifically, we maintain two separate CSR arrays for variable-to-factor messages and factor-to-variable messages, respectively. Each CSR representation consists of two key auxiliary arrays for indexing. One is the row offset array that represents the starting position of each row. The other is the column index array that contains the column indices of the corresponding value. To facilitate efficient access patterns, we align the CSR rows with the computation direction. When computing a factor-to-variable message \(\mu_{a \to v}\), all required variable-to-factor messages associated with factor \(a\) are stored contiguously by using factor identifiers as row indices. Conversely, when computing \(\mu_{v \to a}\), we organize messages using variable identifiers as row indices. This layout allows each thread to access all necessary inputs via sequential memory reads, reducing indexing overhead.

For BP with local structures, instead of storing factor functions as dense tables, we store them as a flattened array of sub-message indicators to effective support BP with local structures. We represent the indicator array as \(indicator[N][L]\), where \(N = |E|\) is the number of edges and \(L = \max_{a \in A} l_a\) is the maximum number of terms among all factors. For efficient GPU storage and access, we organize this array in a CSR-like format using factor identifiers as row indices, so that all indicators associated with the same factor are stored contiguously. For each edge \(e = (a, v)\), \(indicator[e]\) is an integer array of length \(L\). For \(i \leq l_a\), \(indicator[e][i]\) follows the original semantics in 3.3. That is, \(indicator[e][i] = c\) if \(v\) is constrained to value \(c\) in the \(i\)-th term, and \(indicator[e][i] = -1\) otherwise. For \(i>l_a\), we pad \(indicator[e][i]\) with a sentinel value (e.g., -2) indicating an invalid term. This representation avoids storing factor tables in exponential size and enables each thread to perform BP with local structures, significantly reducing memory footprint and improving computational efficiency.

4.3.3 Kernel Launch↩︎

To execute message updates efficiently on GPUs, we organize computation into batches derived from dependency analysis. Each batch consists of messages that can be computed independently without synchronization. For each batch, we launch a GPU kernel where each thread is assigned a single message computation. Each thread is provided with the corresponding factor, variable, and edge identifier. The edge identifier is computed using the CSR representation and determines the location for storing the output message. These identifiers allow each thread to efficiently locate input messages and write results without contention. Batching not only reduces kernel launch overhead but also improves efficiency by exposing massive parallelism.

5 Experimental Evaluation↩︎

In this section, we aim to answer the following research questions:

1. How efficient is FastLBP compared with state-of-the-art approaches?

2. How important is supporting flexible update strategies for program analysis?

We describe our experimental setup in 5.1, and answer the above questions in 5.2 and 5.3, respectively.

5.1 Experimental Setup↩︎

We implement FastLBP in C++ and CUDA, reusing the frontend of libDAI [27] for reading factor graphs. We conduct all experiments on a Linux machine with 2.40 GHz Intel Xeon processors, 256 GB RAM, and NVIDIA GeForce RTX 4090 GPUs.

5.1.1 Benchmarks↩︎

We use the same benchmarks as in [26]. Specifically, we evaluate our approach on two probabilistic program analysis applications: SmartFL [9] and BINGO [6]. Both systems construct large factor graphs derived from program analysis and perform probabilistic inference using BP.

For SmartFL, we use the Defects4J benchmark [30]. After filtering invalid cases where the factor graphs cannot be constructed by SmartFL, we obtain 194 factor graphs. On average, each benchmark contains about 102K lines of code, and the generated probabilistic model contains approximately 258K variables and 812K edges.

For BINGO, we use the benchmark consisting of eight Java programs collected from previous work [31], [32] and the DaCapo benchmark suite [33]. BINGO performs a flow- and context-sensitive datarace analysis [34] on the benchmark. The average program size is about 81K lines of code. The corresponding factor graphs generated by BINGO contain on average 62K variables and 159K edges.

5.1.2 Baselines↩︎

We compare FastLBP against two baselines: [26] and PGMax [21]. [26] implemented their approach on SmartFL and libDAI, evaluating it on SmartFL and BINGO, respectively. Both implementations are single-threaded. We adopt these implementations as the state-of-the-art CPU-based baselines. PGMax [21] is a Python/JAX [35] package that implements BP for discrete factor graphs with arbitrary shapes. It only supports synchronous BP with a fixed number of iterations. PGMax can run on GPUs since JAX supports GPU acceleration. We adopt PGMax as the state-of-the-art GPU-based baseline in our evaluation.

5.2 Efficiency↩︎

5.2.1 Efficiency on SmartFL↩︎

To evaluate the efficiency of FastLBP, we compare it with [26] and PGMax [21] on SmartFL. Since the original setting of SmartFL adopts synchronous BP, all approaches use synchronous BP in this experiment to ensure a fair comparison. The number of message-passing iterations is set to 500, at which point BP has converged on most factor graphs (190/194).

Overall, FastLBP significantly outperforms both baselines. 7 shows the inference time against the number of variables in the factor graphs. Compared to [26], FastLBP achieves a \(17.42\times\) speedup on average (geometric mean). The advantage becomes more pronounced as the size of the factor graph increases. For small graphs (around \(10^2\) to \(10^3\) variables), the performance differences between them are relatively small due to the limited amount of parallelism available. However, as the number of variables grows beyond \(10^4\), FastLBP exhibits significantly lower inference time than [26]. FastLBP also outperforms PGMax, achieving an average speedup of \(6.14\times\). This improvement mainly comes from FastLBP’s efficient implementation of BP with local structures. We note that Chart-15 is not shown in 7 because PGMax runs out of memory on this factor graph. This happens because PGMax adopts a tabular representation of factors, while Chart-15 contains factors connected to many variables.

We also check the correctness of FastLBP on SmartFL. As shown in 2, FastLBP is mathematically equivalent to sequential implementations, but minor floating-point errors might still occur due to implentation details. We quantify the differences in marginal probabilities using relative error. Across all test cases, FastLBP has relative errors less than \(10^{-8}\) compared to [26], confirming its correctness.

5.2.2 Efficiency on BINGO↩︎

We compare FastLBP with [26] on BINGO. BINGO performs probabilistic inference with user feedback as evidence. In each iteration, it queries the user about the alarm with the highest likelihood, then updates the posterior distribution based on the feedback. Efficient probabilistic inference is therefore crucial for quickly identifying true alarms. We evaluate the overall efficiency of inpecting true alarms. We adopt the asynchronous BP introduced in 2.1, which is also the update strategy originally used by BINGO. Empirically, BP does not converge in most factor graphs generated by BINGO, making asynchronous update strategy more suitable for this workload. We do not compare PGMax in this experiment because it only supports synchronous BP. We follow the original stopping criterion, limiting the number of message-passing iterations to between 500 and 1,000 when BP does not converge.

8 shows the number of true alarms inspected over the accumulated inference time for each benchmark. On most benchmarks, FastLBP identifies all true alarms in significantly less time. The only exception is weblech, where the factor graph is relatively small (with 696 variables). Under the constraints of the asynchronous update strategy, FastLBP cannot achieve high parallelism on such small graphs. The benefits of GPU parallelism do not outweigh the overhead introduced by dependency analysis and memory transfers. As the graph size increases, the advantage of FastLBP becomes increasingly evident. On hedc, FastLBP is slightly faster than [26]. On the remaining benchmarks (luindex, jspider, avrora, xalan, sunflow, and ftp), FastLBP runs several times faster. Although FastLBP is marginally slower on small graphs, this has negligible practical impact because the total BP time on such graphs is inherently small (less than 3 seconds on weblech). Thus, the performance difference does not significantly affect user experience. In contrast, users are primarily concerned with execution time on large-scale graphs, where the total inference time can extend to several hours or even exceed one day. It is precisely in these large-scale scenarios that FastLBP excels, significantly accelerating BP and substantially reducing user waiting time.

FastLBP achieves an average speedup of \(2.82\times\) in total inference time compared to [26]. The relatively moderate overall speedup is mainly due to two reasons. First, the asynchronous update strategy required by BINGO limits the amount of parallelism that can be exploited by GPU execution. Second, several benchmarks contain relatively small factor graphs. When excluding benchmarks where [26] finishes within 15 minutes (i.e., weblech, hedc, luindex, and jspider), the average speedup increases to \(5.12\times\), with a maximum speedup of \(10.46\times\). These results demonstrate that FastLBP effectively supports asynchronous update strategies required by program analysis applications.

5.3 The Importance of Supporting Flexible Update Strategies↩︎

To investigate the importance of supporting flexible update strategies in program analysis, we evaluate how different BP update strategies affect the quality of inference results on BINGO. In particular, we compare asynchronous BP used by FastLBP and [26] with synchronous BP implemented in PGMax. We do not evaluate on SmartFL because synchronous BP already converges reliably on SmartFL. We use the default stopping criterion of BINGO for FastLBP and [26], which stops at 500 to 1,000 iterations based on the convergence rate. Since PGMax runs BP for a fixed number of iterations, we set the number of iterations to 1,000 to match the upper bound of BINGO’s stopping criterion.

Following BINGO, we employ three metrics to quantify the quality of inference results: inversion count, Rank-100%-T, and Rank-90%-T. An inversion is defined as a pair of alarms where a false alarm is ranked higher than a true alarm in the user’s inspection sequence. Formally, let \(l_1, \ldots, l_n \in \{0, 1\}\) denote the ground-truth labels of the alarms provided by the user. Then, the inversion count is calculated as: \(\text{Inversion}(l_1, \dots, l_n)=\sum_{i=1}^n (1 - l_i) \sum_{j=i+1}^n l_j\). Rank-100%-T and Rank-90%-T represent the numbers of interactions to identify all and 90% of the true alarms. Lower values indicate that users can prioritize and inspect true alarms more rapidly. Consequently, higher accuracy in BP results correlates with improved alarm ranking and reduced inspection effort.

2 shows the results of the three approaches. Compared with [26], FastLBP shows only negligible differences in all three metrics. Theoretically, FastLBP should produce identical results to the sequential implementation under the same update strategy. The small discrepancies observed in practice are primarily due to floating-point errors between the two implementations. This leads to differences in the marginal probabilities, thus affecting the final alarm rankings. These results confirm that FastLBP preserves the inference quality of the sequential implementation while enabling GPU acceleration.

In contrast, PGMax exhibits significantly worse results across all metrics. The inversion count, Rank-100%-T, and Rank-90%-T increase by \(256.14\%\), \(186.68\%\), and \(128.48\%\), respectively. This degradation significantly affects user experience because users must inspect many more false alarms before discovering true bugs. The main reason is that PGMax only supports synchronous BP. In program analysis frameworks such as BINGO, BP often fails to converge and tends to exhibit more oscillating behaviors under synchronous BP than asynchronous BP.

These results highlight the importance of supporting flexible update strategies for program analysis. In particular, asynchronous BP is essential for workloads where synchronous BP leads to unstable inference results. By supporting both synchronous and asynchronous update strategies, FastLBP can adapt to different workloads while maintaining both efficiency and accuracy.

6 Related Work↩︎

Our work is related to research on (1) GPU acceleration for BP, (2) parallel and distributed BP on CPUs, and (3) probabilistic program analysis using BP. We discuss the related work below.

GPU acceleration for BP. GPU acceleration of BP has been explored in several domain-specific applications, such as computer vision [16], [17], [36] and channel coding [19], [20], [28]. These approaches design GPU implementations tailored to the characteristics of PGMs in specific domains. In contrast, our work focuses on accelerating BP in the context of program analysis, where the resulting PGMs are significantly larger in scale and more complex in graph structure. Beyond domain-specific applications, [21] implement BP using JAX [35], enabling GPU acceleration for general PGMs. Compared to their work, our approach supports more sophisticated update strategies and a more efficient BP algorithm. Randomized Belief Propagation [18] randomly selects a subset of messages to update in each iteration, effectively balancing convergence behavior and parallel efficiency. In contrast, our approach supports a broader range of update strategies and can be easily adapted to incorporate their update strategy.

Parallel and distributed BP on CPUs. Synchronous BP is straightforward to parallelize. To improve convergence, prior work on CPU-based parallel BP has focused on frontier-based update strategies [24], [25]. They design update strategies that use the concept of residuals [23] and update the most influential messages during each iteration for both parallelism and convergence behavior. While effective on CPUs, these approaches incur significant overhead on GPUs [18], limiting their applicability in GPU settings. In contrast, we propose a framework that enables flexible update strategies through a unified representation and dependency analysis, which can also be naturally extended to CPU execution.

Probabilistic program analysis using BP. Recent years have seen increasing interest in applying probabilistic reasoning to program analysis. In particular, Bayesian program analysis transforms static analysis derivation graphs into PGMs and reasons about the likelihood that a reported alarm corresponds to a real bug. Prior work in this area can be classified in three directions. The first line of work focuses on incorporating diverse forms of posterior information to refine the alarm ranking. Posterior information includes user feedback on the alarms [6], differences between the derivations of alarms [7], test execution results [8], and informal information in the program [12]. The second line explores how to refine Bayesian networks to mitigate false generalization on posterior information. [10] proposed to learn new rules and probabilities. [11] and [13] aimed at learning to select a suitable abstraction. The last line applies probabilistic program analysis to augment other software engineering tasks, such as fault localization [9], [14] and fuzzing [15]. These approaches typically model program semantics as PGMs and perform BP for marginal inference. However, the resulting models are often large and complex, making BP a major performance bottleneck. Our work can accelerate these systems by providing an efficient and accurate inference library.

7 Conclusion↩︎

We present FastLBP, a GPU-accelerated Belief Propagation framework for program analysis. By introducing a unified representation for update strategies and a dependency analysis algorithm, FastLBP enables parallel execution for flexible update strategies. Additionally, our memory-efficient GPU implementation of BP with local structures effectively handles logical constraints in program analysis. Our evaluation on SmartFL and BINGO demonstrates that FastLBP significantly outperforms both CPU-based and GPU-based approaches in efficiency without sacrificing accuracy.

Data Availability Statement↩︎

The artifact for this paper is available at Zenodo [37]. It includes source codes, scripts, and data required to reproduce the results presented in 5.

References↩︎