# Performing Bayesian Network Inference Using Amortized Region Approximation with Graph Factorization

Authors: Peng Lin; Changsheng Dou; Nannan Gu; Zhiyuan Shi; Lili Ma
Journal: International Journal of Intelligent Systems
 (2023)

---

## Abstract

Exact inference for large, directed graphical models, also known as Bayesian networks (BNs), can be intractable as the space complexity grows exponentially in the tree-width of the model. Approximate inference, such as generalized belief propagation (GBP), is used instead. GBP treats inference as the Bethe/Kikuchi energy function optimization problem. The solution is found using iterative message passing, which is inefficient and convergent problematic. Recent progress on amortized technique for GBP is an attractive alternative solution that can optimize the Bethe/Kikuchi energy function using (deep) neural networks, requiring no message passing. Despite being efficient, the amortized technique for GBP is applied to undirected graphical models with specific structures and factors, with no guarantee of the approximation quality. This is because the energy function to be amortized is defined by a region (or factor) graph that is ad hoc and difficult to construct to ensure sensible approximations. This paper proposes a new amortized GBP algorithm applied to BN for efficient inference. The proposed algorithm is composed of the following: (i) a new pairwise conversion (PWC) algorithm that converts all the conditional probability distributions in the BN into pairwise factors to facilitate efficient region graph constructions; (ii) following PWC, an improved loop structured region graph (LSRG) algorithm was derived to generate a valid region graph satisfying desired regional properties; and (iii) the energy function defined by the proposed PWC-LSRG region graph can be directly amortized using (deep) neural networks to ensure sensible approximations. Empirical studies show that the proposed amortized PWC-LSRG algorithm is of practical use and significantly improves convergence and efficiency compared to conventional algorithms.


---

## Body

## 1. Introduction and Motivation

Directed graphical models, also known as Bayesian networks (BNs) [1], encode causal and dependent information among random variables. BN is considered a nature-inspired graphical tool for human probabilistic reasoning [2]. In recent years, applications driven by BN modeling have increased rapidly in medical [3], finance [4], biology [5], cybersecurity [6], and forensics and law [7]. Meanwhile, the BN scale also grows in the number of variables and the complexity of the model, making exact probabilistic inference methods, such as the junction tree algorithm [8], intractable. One must rely on approximate methods instead.One important class of approximate inference is variational inference (VI) [9, 10], designed explicitly for BNs. An alternative approach, called generalized belief propagation (GBP) [11], is a variant of VI originally applied to undirected graphical models. VI and GBP are both optimization-based methods, which turn the original global inference problem into an optimization problem subject to local terms with relaxed space complexity constraints. VI is an analytical solution and provides an effective evidence lower bound (ELBO) for the objective function. In contrast, GBP cannot yield an ELBO but is empirically more accurate than VI for marginal inference tasks. Besides, GBP is convenient for automated inference when all continuous variables are discretized [12] into discrete variables.However, the accuracy of GBP depends on constructing a region graph from the model, which is ad hoc and difficult for users. Representative region graph algorithms to alleviate this problem include iterative joint graph propagation (IJGP) [13], loop structured region graph (LSRG) [14], and triplet region construction (TRC) [15]. The TRC is a hybrid algorithm using both structural and factor information of the model and is more accurate than IJGP and LSRG under bounded cluster size. However, TRC is designed explicitly for BNs, so not directly applicable to Markov networks. Instead, LSRG and IJGP use structural information only but applied to general forms of factors, such as conditional random fields encoding both directed and undirected factors [1].The space complexity for VI or GBP is low. However, running VI or GBP is inefficient for large and none-conjugate exponential family models. This is because, under these circumstances, they feature an iterative VI [16] or GBP [11, 17] message-passing procedure, in which computing gradients can be inconvenient [18] and sensitive to numerical issues [15].An attractive alternative to improving efficiency is using (deep) neural network (DNN) methods to perform approximate inference. The DNN methods for graphical models can be divided into variational inference-based DNN (VI-DNN) or belief propagation-based DNN (BP-DNN). Much VI-DNN work involves learning a directed graphical model with structured data modeling [19–21], in which latent variables are involved. Hence, Monte Carlo sampling is used to estimate posterior distributions. The VI-DNN inference is carried out by mapping the ELBO function to a loss function of the neural network. Parameters of interest are then amortized [22] by the neural network parameters. The original ELBO optimization is amortized by training a (deep) neural network, which outputs the ELBO estimation of the interests. So the VI-DNN inference is also known as amortized variational inference. The drawback of VI-DNN is that the estimation variance caused by sampling can be high [23].The BP-DNN is designed for undirected graphical models that employ the (deep) neural network in two ways. The first approach uses a graph neural network (GNN) [24–26], directly mapping the message-passing nodes used in GBP into neural network nodes. The simulation is an analogy to the message-passing procedure. The second approach is amortizing [27–29] the energy function, namely Bethe or Kikuchi function [11] defined by the BP or GBP, using reparametrized neural network loss function and parameters. The Kikuchi function generalizes the Bethe function so the amortized Bethe/Kikuchi are generalized to amortized Kikuchi. Compared to the GNN methods, the amortized technique is more efficient because it directly minimizes the objective function without using message passing. As a result, the (deep) neural network output the predicted marginal belief and an estimated “none-ELBO” upon convergence. However, unlike the VI-DNN, the BP-DNN does not ensure the estimated “none-ELBO” is of any bound to the ground truth but compares more favorably to the VI-DNN on accuracy [28]. Besides, the BP-DNN is more flexible than VI-DNN on parameter settings, because, for example, VI-DNN restricts the factorization of the latent variable distributions into simple forms, such as mean-field factorization [9, 18].Motivated by both the accuracy and efficiency gain obtained by BP-DNN compared to discussed methods, we propose to employ amortized Kikuchi approximation rather than amortized variational inference for BN inference problems. However, despite promising, the amortized Bethe or Kikuchi are only applied to undirected graphical models, such as the Markov random fields (MRF), with pairwise and singleton factors [28, 29] and specific model structures. This is because general MRFs are assumed to be decomposable into pairwise MRFs. Moreover, the accuracy of amortizing depends on how the Kikuchi energy function is defined.Therefore, using the amortized Kikuchi approximation for BN inference involves resolving three critical issues: (i) reducing the size of the conditional probability distribution (CPD) in the BN to a minimal level so the region graph algorithm can be used conveniently. This is because the region size is a lower bound by the CPD (factor) size [15]. (ii) Improve the region graph construction algorithm to ensure sensible approximations for BNs. (iii) DNN is sensitive to parameter tuning, so the justification of the optimal parameter setting is also crucial.This paper assumes all discrete BNs given continuous variables can be discretized [12]. We propose a new amortized Kikuchi approximate inference algorithm for BN inference tasks. The novelty of this work is fourfold:(1)
The amortized Bethe or Kikuchi approximation is normally applied to the undirected graphical model with specific model structures. We propose to extend the amortized Kikuchi approximation for use in BN inference.(2)
We propose a new pairwise conversion (PWC) algorithm applicable to conditional probability distributions in BNs, that can transform the BN into a unified pairwise model representation. The converted BN provides significant flexibility, allowing us to build a region graph using the maximum region size from two to higher orders. Furthermore, this algorithm can be used independently.(3)
We propose a new region graph algorithm, which improves the loop structured region graph algorithm (LSRG) by incorporating node ordering information in a BN, resulting in more accurate region graph approximations than conventional algorithms. This algorithm can also be used independently.(4)
Combining 1–3, we obtain an amortized PWC-LSRG algorithm that significantly improved efficiency and convergence compared to competing inference algorithms.This paper is an extended version of the conference paper [30] with significant new technical improvements over the amortized PWC-LSRG algorithm and more empirical validations. The paper is organized as follows: Section 2 introduces the related BN inference methods. Section 3 develops our proposed amortized PWC-LSRG algorithm. Finally, Section 4 carries out the experiments, and Section 5 concludes the paper.


## 2. BNs and Inference Methods

### 2.1. BN and Exact Inference

In this section, we briefly review the BN notations that will be used in the paper and describe why it is sufficient to solve the BN inference problem using a complete BN structure. Following that, we discuss the relation between cluster size and the size of the conditional probability distributions because cluster size is a measure of the space complexity of both exact and approximate inference.A BN is a directed acyclic graph (DAG)G=<V,E>, with nodes V=1,…,m representing random variables and a set of edges E=s,t:s,t∈V. For each i∈V, Xi is a random variable.Two variables,Xi and Xj, are conditional independent (CI) given K⊆V\i,j, noted Xi⊥XjXK: PXi=xiXj=xj,XK=xK=PXi=xiXK=xk. The key property of a BN is the ordered Markov property which defines an ordering πV for nodes such that a node only depends on its immediate parents, Xi⊥XprediXpai, where pai are the parents of the node and i and predi are the predecessors of i. A conditional probability distribution (CPD) is associated with each node, depending on its parent nodes, if there are any. The BN represents the joint distribution, p of the random variables X1, . . . , Xm as the product of their associated CPDs, as shown as follows:(1)pX1,…,XmG=∏i=1mpXiXpai,where each term pXiXpai is a CPD. The number of variables involved in a CPD is the CPD size.In the worst case, there are no CI assumptions in the BN; the BN graph is a complete DAG (or complete BN) withm nodes, i.e., every pair of nodes is connected by a directed edge. Performing inference on such a BN graph represents the worst case space complexity OKm if each node has K states and is usually intractable for exact methods. We say a BN is a sparse BN if it encodes CI information.Theoretically, any BN model can be represented by a complete BN (by removing CI assumptions). Hence, any algorithm that performs efficient inference for complete BN models will also be efficient for arbitrary BN models. Therefore, it is sufficient to find an efficient approximate inference algorithm for complete BN models.Each CPD functionpXiXpai in a complete DAG G is a function of the variables, it contains pXiXpai∼fiX1:i, as shown in Figure 1(a). In Figure 1(b), the BN contains CI information, so the number of connected edges is less than 1(a). However, in Figure 1(a), BN suffices to represent the BN in Figure 1(b) by setting the removed edges with uniform distribution.Figure 1
(a) A 5-dimensional complete BNG with the CPDs’ list aside each node; (b) a 5-dimensional sparse BN.
(a)(b)Theorem 1.
Any sparse BNG of m nodes can be converted to an equivalent m dimensional complete BN G′, with the equivalence pG=pG′.Proof.
We can always set CPD for the edge inG′ that is missed in G, a uniform distribution. Then, following the chain rule (1), such uniform CPDs will not change the joint distribution of the BN.
Theorem1 ensures that we can always convert a sparse BN to its complete BN format, ensuring we obtain a unified representation for all BNs.
To perform exact inference for BNs, we use the Junction tree (JT) algorithm [8]. The JT construction involves moralization, triangulation, and node elimination processes. The moralization converts the BN G to an undirected graph (also known as Markov network) by connecting an edge Xi,Xj if there is an edge between Xi and Xj in BN G, or if Xi and Xj are parents of the same child node. The added edge between the parents that share the same child node is called a moral edge. The undirected graph is called triangulated if every cycle of length strictly greater than 3 possesses a chord.
After performing these steps, we can obtain a maximal spanning tree as a JT if the triangulation is optimal. Each vertex in the JT is a cluster containing a set of BN variables. The space complexity of the JT is measured by the maximum cluster size −1, called the tree-width of a BN. The maximum cluster size controls the space complexity and is lower bounded by the maximum CPD size. This means the maximum cluster size is always larger than the CPD size. For a complete BN, as shown in Figure1(a), the maximum CPD size always equals the maximum cluster size (under JT) so there is no space for efficiency improvement by running approximate inference for complete BNs, unless the maximum CPD size is reduced.
Therefore, when we rely on approximate inference, we first need to reduce the CPD size of the BN to a minimum level so that the cluster size can be reduced in approximate methods. Then, we can further reduce the maximum cluster size using region graph-based methods. We address the CPD size reduction in section3.1.


### 2.2. Generalized Belief Propagation and Region Graph Construction

This section reviews generalized belief propagation, the associated Bethe/Kikuchi energy functions, and essential region graph construction methods.


#### 2.2.1. BFE/KFE Functions

In the early work, the notation “Belief Propagation” (BP) refers to Bayes updating rules for directed acyclical graph models [31]. However, we restrict the BP algorithm as described in Yedidia (2005) [11], which is a formulation of using factor graph representation for BNs.For a BN with joint distributionpX1,…,Xm, a factor graph representation is a bipartite graph connecting the factor nodes a (a∈Φ, with Φ be all factors/CPDs in the BN, and the variable nodes i (i=1,…,m). A factor node a has a factor φaXa, in which Xa are variables in the factor φa, and the length of Xa equals the factor size.The joint distributionpX1,…,Xm can then be expressed as follows:(2)pX1,…,Xm=1Z∏a=1ΦφaXa,where Z is the normalization constant or partition function. Using Boltzmann’s law [11], we can formularize (2) as the probability of a nonphysical system at state x, pX=x=1/ZT∙e−Ex/T. Here, T is the temperature and EX=x=−∑a=1ΦlnφaXa is the energy term of the system at state x. The partition function ZT is calculated as ZT=∑x∈Se−Ex/T, where S is the space of all possible states of the system.WhenZT cannot be computed exactly, we can use a belief bX to construct an approximation:(3)FBFEb=Ub−Hb=−lnZ+Dbp,where FBFEb is the Bethe free energy (BFE) [11], Ub=∑x∈SbxEx and Hb=−∑x∈Sbxlnbx are variational energy and variational entropy, respectively. When the Kullback–Leibler distance Dbp=0, the BFE achieves its minimal and recovers the partition function Z. Therefore, minimizing Fbetheb approximates −lnZ when Dbp≠0. The beliefs contain limited information exchange in the factor graph representation since the factor nodes are connected to singleton variable nodes, with only singleton node interactions captured.Compared to BP, GBP is more general as it allows higher-ordered node interactions during the approximation. This is achieved by using a region graph instead of a factor graph. For an undirected graphical model with joint distributionpX1,…,Xm, (2) is rewritten as follows:(4)pX1,…,Xm=1Z∏r∈RϕrXr,where R is a collection of subsets of variables V; we call each r∈R a region. Associated with each region, r is a non-negative kernel function ϕrXr. A region graph GR is a directed graph in which each vertex v∈GR corresponds to a region r. We say vp is an ancestor of vc if vc⊂vp, and we say vp is a parent of vc if vp⟶vc is a directed edge.Unlike BP which optimizes the BFE function, GBP optimizes the Kikuchi cluster free energy (KFE) [11] function. However, the BFE can be viewed as a special case of the KFE. Optimizing the KFE involves defining a factored energy function regarding entropies over all regions as an approximation of FKFE, and enforcing a set of regions to obey consistency on the local polytope (rather than global marginal polytope). Specifically, the key components of the GBP are as follows:(1)
Define a collection ofR-pseudo-beliefs brXr,r∈R;(2)
Construct a constrained minimization problem of the form:(5)prXr≅br∗Xr≔argminbrXr∈LRFKFEbrXr,where(6)FKFEbrXr≔∑r∈Rcr∑XrbrXrErXr+∑XrbrXrlnbrXr,and the Lagrangian term LR in (5) ensures belief consistency between parent-to-child regions, and the value in the belief is normalized:(7)LR≔brXr,r∈R:∀r,u∈Rs.t.u⊂r,∑Xr\ubrXr=buXu∀u∈R,∑XubrXu=1.The termEr in (6) represents the energy associated with each region r. cr=1−∑r^∈Ancestorrcr^ is a counting number of each region, given cr^ is degrees of freedom for the region r^.Optimizing (6) with constraints (7) can be achieved by using GBP message passing or its convergent alternative concave-convex procedure (CCCP) [17]. The benefit of using the Kikuchi approximation is that we can use a collection of bounded-sized regions in equations (5–7) to avoid exponentially sized clusters in exact methods, such as JT. The JT can also be viewed as the special case of Kikuchi approximation when the GBP cluster size increases to tree − width + 1 [32], in which the partition function Z exactly recovers the log probability of the model logpe [1] (also called the probability of evidence).


#### 2.2.2. Region Graph Algorithms

The cluster variation method (CVM) [11] is typically used to construct a region graph, ensuring each node is counted exactly once. The CVM is a multilevel region graph with the first level regions, called the outer regions, absorbing all factors in the factor graph. The subsequent levels of regions are inner regions, which intersect the previous level regions. All outer regions use a default counting number 1, and all inner regions’ counting numbers are calculated using the counting number definition. Figure 2 compares a factor graph Gf and a CVM region graph GR for Figure 1(b) BN. The CVM region graph encodes higher-ordered interactions among variables than the factor graph. So, region graph-based methods can be more accurate than factor graph-based methods.Figure 2
(a) Factor graphGf and (b) CVM region graph GR of Figure 1(b) BN G.
(a)(b)However, CVM does not define the outer regions and only ensures constructing a valid region graph when the outer regions are given. The outer region selection algorithms can be grouped into structural, factor, and hybrid information-based algorithms. The loop structured region graph (LSRG) [14] represents the structural information-based approach to selecting outer regions solely based on the model structure. The triplet region construction (TRC) [15] is the latest hybrid outer region selection algorithm using the model’s structural and factor information. TRC is designed explicitly for BNs (not MNs) in which all CPDs are converted into triplets. TRC is more accurate than LSRG but less flexible than LSRG to apply to general forms of factors. For example, a conditional random field [1] is a partial BN encoding both CPDs and Markov factors. Also, if we are keen to use outer regions of size two, LSRG is more convenient. We omit the factor information-based region selection algorithms because they are significantly inefficient compared to LSRG and TRC.The merit of LSRG is that it satisfies two important region properties: maxent-entropy normal and tree robustness.The maxent-normal property requires that if a constrained region-based free energy approximation is valid, the corresponding region-based entropy achieves its maximum when all the beliefs are uniform. The CVM algorithm can satisfy this property when the outer regions are determined.The tree-robustness property requires that a structured region graph can be reduced so that inference on any spanning tree ofG is exact. Therefore, the LSRG selects cycle faces from an undirected graph model G as outer regions by ensuring any spanning tree of G runs exactly under these selected outer regions. However, this property is difficult to satisfy for general structured models. So, the LSRG [14] only ensures the tree-robustness property for two typical graph structures: the planar MN and the complete MN, in which each edge of the MN encodes a pairwise factor.Definition 1.
A tree-robustness region graph [14] for a complete MN is a CVM region graph with the outer regions defined by the following: (i) choose a root node i of the MN. (ii) Then, from the node i, create a spanning tree of the MN, with the off-tree edges j,k. (iii) Create i,j,k as the outer regions.
The region graph for a complete MN obtained by Definition1 satisfies the tree-robustness property. When the subgraphs of MN G are mutually singly connected planar and complete MNs, the LSRG can also ensure the tree-robustness by constructing outer regions for each subgraph. For more general MNs, if the structure is too complex to admit the tree-robustness property, the LSRG will sacrifice finding a sequence of outer regions that is tree exact for many spanning trees of G instead of all spanning trees. The construction is, in turn, satisfying the partial tree-robustness.
The LSRG is promising but has several limitations when applied to BNs. First, BN is a directed graph model, so choosing which node as the root node is undefined by LSRG. The cycle size found by LSRG is not equal, presenting difficulties in controlling the space complexity. We will address these issues in section3.2.


## 2.1. BN and Exact Inference

In this section, we briefly review the BN notations that will be used in the paper and describe why it is sufficient to solve the BN inference problem using a complete BN structure. Following that, we discuss the relation between cluster size and the size of the conditional probability distributions because cluster size is a measure of the space complexity of both exact and approximate inference.A BN is a directed acyclic graph (DAG)G=<V,E>, with nodes V=1,…,m representing random variables and a set of edges E=s,t:s,t∈V. For each i∈V, Xi is a random variable.Two variables,Xi and Xj, are conditional independent (CI) given K⊆V\i,j, noted Xi⊥XjXK: PXi=xiXj=xj,XK=xK=PXi=xiXK=xk. The key property of a BN is the ordered Markov property which defines an ordering πV for nodes such that a node only depends on its immediate parents, Xi⊥XprediXpai, where pai are the parents of the node and i and predi are the predecessors of i. A conditional probability distribution (CPD) is associated with each node, depending on its parent nodes, if there are any. The BN represents the joint distribution, p of the random variables X1, . . . , Xm as the product of their associated CPDs, as shown as follows:(1)pX1,…,XmG=∏i=1mpXiXpai,where each term pXiXpai is a CPD. The number of variables involved in a CPD is the CPD size.In the worst case, there are no CI assumptions in the BN; the BN graph is a complete DAG (or complete BN) withm nodes, i.e., every pair of nodes is connected by a directed edge. Performing inference on such a BN graph represents the worst case space complexity OKm if each node has K states and is usually intractable for exact methods. We say a BN is a sparse BN if it encodes CI information.Theoretically, any BN model can be represented by a complete BN (by removing CI assumptions). Hence, any algorithm that performs efficient inference for complete BN models will also be efficient for arbitrary BN models. Therefore, it is sufficient to find an efficient approximate inference algorithm for complete BN models.Each CPD functionpXiXpai in a complete DAG G is a function of the variables, it contains pXiXpai∼fiX1:i, as shown in Figure 1(a). In Figure 1(b), the BN contains CI information, so the number of connected edges is less than 1(a). However, in Figure 1(a), BN suffices to represent the BN in Figure 1(b) by setting the removed edges with uniform distribution.Figure 1
(a) A 5-dimensional complete BNG with the CPDs’ list aside each node; (b) a 5-dimensional sparse BN.
(a)(b)Theorem 1.
Any sparse BNG of m nodes can be converted to an equivalent m dimensional complete BN G′, with the equivalence pG=pG′.Proof.
We can always set CPD for the edge inG′ that is missed in G, a uniform distribution. Then, following the chain rule (1), such uniform CPDs will not change the joint distribution of the BN.
Theorem1 ensures that we can always convert a sparse BN to its complete BN format, ensuring we obtain a unified representation for all BNs.
To perform exact inference for BNs, we use the Junction tree (JT) algorithm [8]. The JT construction involves moralization, triangulation, and node elimination processes. The moralization converts the BN G to an undirected graph (also known as Markov network) by connecting an edge Xi,Xj if there is an edge between Xi and Xj in BN G, or if Xi and Xj are parents of the same child node. The added edge between the parents that share the same child node is called a moral edge. The undirected graph is called triangulated if every cycle of length strictly greater than 3 possesses a chord.
After performing these steps, we can obtain a maximal spanning tree as a JT if the triangulation is optimal. Each vertex in the JT is a cluster containing a set of BN variables. The space complexity of the JT is measured by the maximum cluster size −1, called the tree-width of a BN. The maximum cluster size controls the space complexity and is lower bounded by the maximum CPD size. This means the maximum cluster size is always larger than the CPD size. For a complete BN, as shown in Figure1(a), the maximum CPD size always equals the maximum cluster size (under JT) so there is no space for efficiency improvement by running approximate inference for complete BNs, unless the maximum CPD size is reduced.
Therefore, when we rely on approximate inference, we first need to reduce the CPD size of the BN to a minimum level so that the cluster size can be reduced in approximate methods. Then, we can further reduce the maximum cluster size using region graph-based methods. We address the CPD size reduction in section3.1.


## 2.2. Generalized Belief Propagation and Region Graph Construction

This section reviews generalized belief propagation, the associated Bethe/Kikuchi energy functions, and essential region graph construction methods.


### 2.2.1. BFE/KFE Functions

In the early work, the notation “Belief Propagation” (BP) refers to Bayes updating rules for directed acyclical graph models [31]. However, we restrict the BP algorithm as described in Yedidia (2005) [11], which is a formulation of using factor graph representation for BNs.For a BN with joint distributionpX1,…,Xm, a factor graph representation is a bipartite graph connecting the factor nodes a (a∈Φ, with Φ be all factors/CPDs in the BN, and the variable nodes i (i=1,…,m). A factor node a has a factor φaXa, in which Xa are variables in the factor φa, and the length of Xa equals the factor size.The joint distributionpX1,…,Xm can then be expressed as follows:(2)pX1,…,Xm=1Z∏a=1ΦφaXa,where Z is the normalization constant or partition function. Using Boltzmann’s law [11], we can formularize (2) as the probability of a nonphysical system at state x, pX=x=1/ZT∙e−Ex/T. Here, T is the temperature and EX=x=−∑a=1ΦlnφaXa is the energy term of the system at state x. The partition function ZT is calculated as ZT=∑x∈Se−Ex/T, where S is the space of all possible states of the system.WhenZT cannot be computed exactly, we can use a belief bX to construct an approximation:(3)FBFEb=Ub−Hb=−lnZ+Dbp,where FBFEb is the Bethe free energy (BFE) [11], Ub=∑x∈SbxEx and Hb=−∑x∈Sbxlnbx are variational energy and variational entropy, respectively. When the Kullback–Leibler distance Dbp=0, the BFE achieves its minimal and recovers the partition function Z. Therefore, minimizing Fbetheb approximates −lnZ when Dbp≠0. The beliefs contain limited information exchange in the factor graph representation since the factor nodes are connected to singleton variable nodes, with only singleton node interactions captured.Compared to BP, GBP is more general as it allows higher-ordered node interactions during the approximation. This is achieved by using a region graph instead of a factor graph. For an undirected graphical model with joint distributionpX1,…,Xm, (2) is rewritten as follows:(4)pX1,…,Xm=1Z∏r∈RϕrXr,where R is a collection of subsets of variables V; we call each r∈R a region. Associated with each region, r is a non-negative kernel function ϕrXr. A region graph GR is a directed graph in which each vertex v∈GR corresponds to a region r. We say vp is an ancestor of vc if vc⊂vp, and we say vp is a parent of vc if vp⟶vc is a directed edge.Unlike BP which optimizes the BFE function, GBP optimizes the Kikuchi cluster free energy (KFE) [11] function. However, the BFE can be viewed as a special case of the KFE. Optimizing the KFE involves defining a factored energy function regarding entropies over all regions as an approximation of FKFE, and enforcing a set of regions to obey consistency on the local polytope (rather than global marginal polytope). Specifically, the key components of the GBP are as follows:(1)
Define a collection ofR-pseudo-beliefs brXr,r∈R;(2)
Construct a constrained minimization problem of the form:(5)prXr≅br∗Xr≔argminbrXr∈LRFKFEbrXr,where(6)FKFEbrXr≔∑r∈Rcr∑XrbrXrErXr+∑XrbrXrlnbrXr,and the Lagrangian term LR in (5) ensures belief consistency between parent-to-child regions, and the value in the belief is normalized:(7)LR≔brXr,r∈R:∀r,u∈Rs.t.u⊂r,∑Xr\ubrXr=buXu∀u∈R,∑XubrXu=1.The termEr in (6) represents the energy associated with each region r. cr=1−∑r^∈Ancestorrcr^ is a counting number of each region, given cr^ is degrees of freedom for the region r^.Optimizing (6) with constraints (7) can be achieved by using GBP message passing or its convergent alternative concave-convex procedure (CCCP) [17]. The benefit of using the Kikuchi approximation is that we can use a collection of bounded-sized regions in equations (5–7) to avoid exponentially sized clusters in exact methods, such as JT. The JT can also be viewed as the special case of Kikuchi approximation when the GBP cluster size increases to tree − width + 1 [32], in which the partition function Z exactly recovers the log probability of the model logpe [1] (also called the probability of evidence).


### 2.2.2. Region Graph Algorithms

The cluster variation method (CVM) [11] is typically used to construct a region graph, ensuring each node is counted exactly once. The CVM is a multilevel region graph with the first level regions, called the outer regions, absorbing all factors in the factor graph. The subsequent levels of regions are inner regions, which intersect the previous level regions. All outer regions use a default counting number 1, and all inner regions’ counting numbers are calculated using the counting number definition. Figure 2 compares a factor graph Gf and a CVM region graph GR for Figure 1(b) BN. The CVM region graph encodes higher-ordered interactions among variables than the factor graph. So, region graph-based methods can be more accurate than factor graph-based methods.Figure 2
(a) Factor graphGf and (b) CVM region graph GR of Figure 1(b) BN G.
(a)(b)However, CVM does not define the outer regions and only ensures constructing a valid region graph when the outer regions are given. The outer region selection algorithms can be grouped into structural, factor, and hybrid information-based algorithms. The loop structured region graph (LSRG) [14] represents the structural information-based approach to selecting outer regions solely based on the model structure. The triplet region construction (TRC) [15] is the latest hybrid outer region selection algorithm using the model’s structural and factor information. TRC is designed explicitly for BNs (not MNs) in which all CPDs are converted into triplets. TRC is more accurate than LSRG but less flexible than LSRG to apply to general forms of factors. For example, a conditional random field [1] is a partial BN encoding both CPDs and Markov factors. Also, if we are keen to use outer regions of size two, LSRG is more convenient. We omit the factor information-based region selection algorithms because they are significantly inefficient compared to LSRG and TRC.The merit of LSRG is that it satisfies two important region properties: maxent-entropy normal and tree robustness.The maxent-normal property requires that if a constrained region-based free energy approximation is valid, the corresponding region-based entropy achieves its maximum when all the beliefs are uniform. The CVM algorithm can satisfy this property when the outer regions are determined.The tree-robustness property requires that a structured region graph can be reduced so that inference on any spanning tree ofG is exact. Therefore, the LSRG selects cycle faces from an undirected graph model G as outer regions by ensuring any spanning tree of G runs exactly under these selected outer regions. However, this property is difficult to satisfy for general structured models. So, the LSRG [14] only ensures the tree-robustness property for two typical graph structures: the planar MN and the complete MN, in which each edge of the MN encodes a pairwise factor.Definition 1.
A tree-robustness region graph [14] for a complete MN is a CVM region graph with the outer regions defined by the following: (i) choose a root node i of the MN. (ii) Then, from the node i, create a spanning tree of the MN, with the off-tree edges j,k. (iii) Create i,j,k as the outer regions.
The region graph for a complete MN obtained by Definition1 satisfies the tree-robustness property. When the subgraphs of MN G are mutually singly connected planar and complete MNs, the LSRG can also ensure the tree-robustness by constructing outer regions for each subgraph. For more general MNs, if the structure is too complex to admit the tree-robustness property, the LSRG will sacrifice finding a sequence of outer regions that is tree exact for many spanning trees of G instead of all spanning trees. The construction is, in turn, satisfying the partial tree-robustness.
The LSRG is promising but has several limitations when applied to BNs. First, BN is a directed graph model, so choosing which node as the root node is undefined by LSRG. The cycle size found by LSRG is not equal, presenting difficulties in controlling the space complexity. We will address these issues in section3.2.


## 2.2.1. BFE/KFE Functions

In the early work, the notation “Belief Propagation” (BP) refers to Bayes updating rules for directed acyclical graph models [31]. However, we restrict the BP algorithm as described in Yedidia (2005) [11], which is a formulation of using factor graph representation for BNs.For a BN with joint distributionpX1,…,Xm, a factor graph representation is a bipartite graph connecting the factor nodes a (a∈Φ, with Φ be all factors/CPDs in the BN, and the variable nodes i (i=1,…,m). A factor node a has a factor φaXa, in which Xa are variables in the factor φa, and the length of Xa equals the factor size.The joint distributionpX1,…,Xm can then be expressed as follows:(2)pX1,…,Xm=1Z∏a=1ΦφaXa,where Z is the normalization constant or partition function. Using Boltzmann’s law [11], we can formularize (2) as the probability of a nonphysical system at state x, pX=x=1/ZT∙e−Ex/T. Here, T is the temperature and EX=x=−∑a=1ΦlnφaXa is the energy term of the system at state x. The partition function ZT is calculated as ZT=∑x∈Se−Ex/T, where S is the space of all possible states of the system.WhenZT cannot be computed exactly, we can use a belief bX to construct an approximation:(3)FBFEb=Ub−Hb=−lnZ+Dbp,where FBFEb is the Bethe free energy (BFE) [11], Ub=∑x∈SbxEx and Hb=−∑x∈Sbxlnbx are variational energy and variational entropy, respectively. When the Kullback–Leibler distance Dbp=0, the BFE achieves its minimal and recovers the partition function Z. Therefore, minimizing Fbetheb approximates −lnZ when Dbp≠0. The beliefs contain limited information exchange in the factor graph representation since the factor nodes are connected to singleton variable nodes, with only singleton node interactions captured.Compared to BP, GBP is more general as it allows higher-ordered node interactions during the approximation. This is achieved by using a region graph instead of a factor graph. For an undirected graphical model with joint distributionpX1,…,Xm, (2) is rewritten as follows:(4)pX1,…,Xm=1Z∏r∈RϕrXr,where R is a collection of subsets of variables V; we call each r∈R a region. Associated with each region, r is a non-negative kernel function ϕrXr. A region graph GR is a directed graph in which each vertex v∈GR corresponds to a region r. We say vp is an ancestor of vc if vc⊂vp, and we say vp is a parent of vc if vp⟶vc is a directed edge.Unlike BP which optimizes the BFE function, GBP optimizes the Kikuchi cluster free energy (KFE) [11] function. However, the BFE can be viewed as a special case of the KFE. Optimizing the KFE involves defining a factored energy function regarding entropies over all regions as an approximation of FKFE, and enforcing a set of regions to obey consistency on the local polytope (rather than global marginal polytope). Specifically, the key components of the GBP are as follows:(1)
Define a collection ofR-pseudo-beliefs brXr,r∈R;(2)
Construct a constrained minimization problem of the form:(5)prXr≅br∗Xr≔argminbrXr∈LRFKFEbrXr,where(6)FKFEbrXr≔∑r∈Rcr∑XrbrXrErXr+∑XrbrXrlnbrXr,and the Lagrangian term LR in (5) ensures belief consistency between parent-to-child regions, and the value in the belief is normalized:(7)LR≔brXr,r∈R:∀r,u∈Rs.t.u⊂r,∑Xr\ubrXr=buXu∀u∈R,∑XubrXu=1.The termEr in (6) represents the energy associated with each region r. cr=1−∑r^∈Ancestorrcr^ is a counting number of each region, given cr^ is degrees of freedom for the region r^.Optimizing (6) with constraints (7) can be achieved by using GBP message passing or its convergent alternative concave-convex procedure (CCCP) [17]. The benefit of using the Kikuchi approximation is that we can use a collection of bounded-sized regions in equations (5–7) to avoid exponentially sized clusters in exact methods, such as JT. The JT can also be viewed as the special case of Kikuchi approximation when the GBP cluster size increases to tree − width + 1 [32], in which the partition function Z exactly recovers the log probability of the model logpe [1] (also called the probability of evidence).


## 2.2.2. Region Graph Algorithms

The cluster variation method (CVM) [11] is typically used to construct a region graph, ensuring each node is counted exactly once. The CVM is a multilevel region graph with the first level regions, called the outer regions, absorbing all factors in the factor graph. The subsequent levels of regions are inner regions, which intersect the previous level regions. All outer regions use a default counting number 1, and all inner regions’ counting numbers are calculated using the counting number definition. Figure 2 compares a factor graph Gf and a CVM region graph GR for Figure 1(b) BN. The CVM region graph encodes higher-ordered interactions among variables than the factor graph. So, region graph-based methods can be more accurate than factor graph-based methods.Figure 2
(a) Factor graphGf and (b) CVM region graph GR of Figure 1(b) BN G.
(a)(b)However, CVM does not define the outer regions and only ensures constructing a valid region graph when the outer regions are given. The outer region selection algorithms can be grouped into structural, factor, and hybrid information-based algorithms. The loop structured region graph (LSRG) [14] represents the structural information-based approach to selecting outer regions solely based on the model structure. The triplet region construction (TRC) [15] is the latest hybrid outer region selection algorithm using the model’s structural and factor information. TRC is designed explicitly for BNs (not MNs) in which all CPDs are converted into triplets. TRC is more accurate than LSRG but less flexible than LSRG to apply to general forms of factors. For example, a conditional random field [1] is a partial BN encoding both CPDs and Markov factors. Also, if we are keen to use outer regions of size two, LSRG is more convenient. We omit the factor information-based region selection algorithms because they are significantly inefficient compared to LSRG and TRC.The merit of LSRG is that it satisfies two important region properties: maxent-entropy normal and tree robustness.The maxent-normal property requires that if a constrained region-based free energy approximation is valid, the corresponding region-based entropy achieves its maximum when all the beliefs are uniform. The CVM algorithm can satisfy this property when the outer regions are determined.The tree-robustness property requires that a structured region graph can be reduced so that inference on any spanning tree ofG is exact. Therefore, the LSRG selects cycle faces from an undirected graph model G as outer regions by ensuring any spanning tree of G runs exactly under these selected outer regions. However, this property is difficult to satisfy for general structured models. So, the LSRG [14] only ensures the tree-robustness property for two typical graph structures: the planar MN and the complete MN, in which each edge of the MN encodes a pairwise factor.Definition 1.
A tree-robustness region graph [14] for a complete MN is a CVM region graph with the outer regions defined by the following: (i) choose a root node i of the MN. (ii) Then, from the node i, create a spanning tree of the MN, with the off-tree edges j,k. (iii) Create i,j,k as the outer regions.
The region graph for a complete MN obtained by Definition1 satisfies the tree-robustness property. When the subgraphs of MN G are mutually singly connected planar and complete MNs, the LSRG can also ensure the tree-robustness by constructing outer regions for each subgraph. For more general MNs, if the structure is too complex to admit the tree-robustness property, the LSRG will sacrifice finding a sequence of outer regions that is tree exact for many spanning trees of G instead of all spanning trees. The construction is, in turn, satisfying the partial tree-robustness.
The LSRG is promising but has several limitations when applied to BNs. First, BN is a directed graph model, so choosing which node as the root node is undefined by LSRG. The cycle size found by LSRG is not equal, presenting difficulties in controlling the space complexity. We will address these issues in section3.2.


## 3. Amortized BN Inference

This section develops our proposed, amortized PWC-LSRG algorithm. Section3.1 derives the pairwise conversion algorithm for BN CPDs. Section 3.2 presents an improved LSRG algorithm bespoke for BN inference. Finally, section 3.3 depicts the amortized PWC-LSRG algorithm.


### 3.1. Pairwise Decomposition of CPDs

As the CPD size lower bounds cluster size, ideally, the CPD size should be minimized thus the cluster size can also be reduced. The worst but also the most general case is the complete BN model in which all edges encode pairwise CPDs, so the outer region size can be reduced to two (which is the minimal size for an outer region). By Theorem1, any BN (including the partial BN) is a subset of the associated complete BN model. Suppose we can efficiently decompose a BN into all pairwise CPD representations and resolve its associated complete BN inference using region-based methods. In that case, we will resolve general BN inference problems efficiently.Reducing the CPD size can be achieved from either the factor level or the CPD level. Representing methods are binary factorization (BF) [15, 33] and MRF decomposition [34]. The BF algorithm reduces the CPD size by introducing intermediate variables to ensure each node connects no more than two parent nodes. The MRF decomposition [34] can decompose the general MN factors into pairwise factors by introducing an auxiliary node, with the new pairwise factors defined as the fraction power of the original factor. This method is inconvenient to apply to BN when the CPD size is large.Our proposed pairwise decomposition method is a two-step hybrid algorithm. We use BF to reduce the CPD size to three initially. Then, for each triplet CPD, we decompose it into three pairwise CPDs and two associated hidden pairwise factors. The process is exact. The joint distribution of the BN is then expressed as the product of all pairwise factors.We binary factorize the BNG by adding intermediate nodes and obtaining an equivalent BF model G′, in which every node has at the most two parent nodes. Equivalence means that the BF model G′ can always rebuild the CPDs in the BN G through intermediate nodes. It has been proved in [15] that the BF process is exact for discrete BNs, so any ordered joint distributions defined on the same variables between G and G′ are the same. The benefit of using the BF algorithm is to ensure the CPD size is no larger than three.Figure3 presents the BF algorithm for factorizing a 4-dimensional complete discrete BN G into its BF model G′. In general, the CPD for the added intermediate node E has the form:(8)pE=eijX=xp,Y=yq,=1,ifp=iandq=j,0,otherwise,where X and Y are parents of E and xp, (p=1,…,m1) and yq (q=1,…,m2) are node states, and eij (i=1,…,m1, j=1,…,m2) is an m1×m2 CPD matrix for E. Following (8) we can always obtain a binary factorized BN, with a maximum CPD size of three.Figure 3
(a) A 4-dimensional complete discrete BNG with all binary variables; (b) a binary factorized BN G′ of (a). The singleton marginal distributions are obtained using the JT algorithm in AgenaRisk [35].
(a)(b)In general, for a triplet CPDpX3X1,X2 with the associated possible pairwise CPDs pX1Xpa1 and pX2Xpa2, we can introduce an intermediate node Z with the settings in (9) to decompose pX3X1,X2 into pairwise CPDs. The scope of Z is the Cartesian product: SZ=SX1×SX2×SX3. Therefore, Z=zj (j=1,…,SZ) takes the value of all assignments of X1,X2,X3. We have the following:(9)pZ=φ123X1,X2,X3,pX1=x1Z=zj=1,ifx1⊂assignmentj,pX2=x2Z=zj=1,ifx2⊂assignmentj,pX3=x3Z=zj=1,ifx3⊂assignmentj,where φ123X1,X2,X3 is the function of pX3X1,X2. Due to the complementary property of the CPD, other entries not in (9) are set to zero. Because (9) adds new connections Z⟶X1 and Z⟶X2 for constructing the new CPDs, the original CPDs pX1Xpa1 and pX2Xpa2 cannot be explicitly represented in the BN. So, they are set as hidden pairwise factors (as opposed to CPDs) in the inference phase. Note that we cannot multiply the CPD pX1Xpa1 into pX1Z as that breaks the pairwise factor constraint.Figure4 shows an example using (9). The CPD pCA,B in Figure 4(a)G does not associate with pairwise CPDs for the parent nodes A and B, so we simply introduce a new node Z0 in 4(b)G^ following (9). Node C is then disconnected to A and B and connects to Z0 instead. Similarly, we need to introduce a new node Z1 to decompose the CPD pFD,E in 4(b). The CPD pFD,E is associated with two pairwise CPDs pDB and pEC, which needs to be set as hidden pairwise factors, shown as dashed edges in Figure 4(b). After exact inference, both 4(a)G and 4(b)G^ produced identical marginal results.Figure 4
(a) A binary factorized BNG; (b) pairwise converted BN G^ of G. The dashed edges represent hidden CPDs.
(a)(b)Proposition 1.
The pairwise conversion (PWC) (9) is an exact procedure to convert a binary factorized BN G into its equivalent pairwise representation BN G^.Proof.
The CPDpX1=x1Z=zj=1 in (9) can be viewed as an indicator function, i.e., taking the relevant values in φ123X1,X2,X3. So, multiplying all the CPDs in (9) and integrating Z out, we end up with φ123X1,X2,X3, equaling pX3X1,X2. The hidden pairwise factors are original BN CPDs. Therefore, our pairwise converted BN G^ does not change the chain rule (1) for the binary factorized BN G. Hence, pG^=pG.


### 3.2. LSRG Algorithm with Node Ordering Information

We can apply LSRG to the moral graph of a pairwise converted BNG^. However, the LSRG does not incorporate the node ordering information of the BN, i.e., ignoring the ordered Markov property that results in less optimal region graph constructions for BNs. To improve that, we introduce a node ordering π using the BN G^.Definition 2.
The node orderingπ is defined as π≔X1⟶,…,⟶Xm for the m variables in the pairwise BN G^, where for any Xj that is, the successor to Xi, Xj∉Xpai.
In complete BNs, the node orderingπ is fixed in the model. As introduced in section 2.2, Definition 1, the tree-robustness region graph can be obtained conveniently for a complete MN. And a complete MN has the same structure as a complete BN’s moral graph. So, the prescription of Definition 1 also works for complete BNs. For sparse BNs, we can use Definition 2 and Theorem 1, adding directed edges to the BN, obtaining an equivalent complete BN, and using Definition 1 to build a tree-robustness region graph.
For example, following Definition2, we define a node ordering for Figure 4(b)G^ with π≔Z0⟶A⟶B⟶C⟶Z1⟶D⟶E⟶F to obtain an 8-dimensional complete BN G^π as shown in Figure 5(a). The pairwise CPDs that are not contained in G^ encode uniform distributions. We can choose Z0 as the root node for LSRG to generate the outer regions. The resulting outer region size is three regardless of how large the tree-width of the BN is.
Although we can always use the complete BN representationG^π to obtain a tree-robustness region graph for a sparse BN G^, we are left with a dramatically increased number of outer regions along with the model dimensionality. This is because most of the outer regions are added unnecessarily. Also, all the outer regions contain the root node, resulting in possible numerical issues [15] associated with the region graph.
We will therefore generate the LSRG outer regions directly from the moral graph ofG^ instead of G^π, as shown in Figure 5(b). The resulting region graph will only ensure partial tree-robustness rather than tree-robustness. Following that, we obtain a significant reduction in outer regions while retaining sensible approximations. For example, the number of LSRG outer regions will be reduced from 21 to 6 when using the Figure 5(b) model rather than Figure 5(a) model.
After performing the PWC algorithm, all BN nodes will be connected to no more than two nodes. So, the subgraphs of the BN will be of simple forms, such as planar graphs and/or four-dimensional complete DAGs. It is impossible to contain a complete DAG larger than four dimensions in the subgraph because it violates the PWC algorithm.
If we bound the maximum region size to two, we can use the node pairs connecting an edge in the moral graph as outer regions. In this case, the Kikuchi approximation equals the Bethe approximation. If we bound the maximum region size to three, we will obtain LSRG outer regions with triplet cycles and cycles of length>3. The outer regions found by LSRG are not equal sized. So, the length >3 outer regions will need further decompositions. This is one limitation that the LSRG does not control the outer region sizes.
Our improved LSRG algorithm contains two more steps: (i) add chordal edges to cycles with length>3 to make it chordal. (ii) The chordal edges are added in a sequence following the node ordering π.
For example, in Figure5(b), there is a length four cycle face Z0,B3,Z1,C4 found by LSRG. Following the node ordering π, a chordal edge being added is B3,C4 rather than Z0,Z1 because Z1 is ordered after B3 and C4. The chordal edge connecting B3 and C4 has higher priority than the edge connecting Z0 and Z1. This is one instance of the ordered Markov property that the BN holds for conditional independence Z1⊥Z0pB3,C4. So, the set of outer regions is Z0,B3,C4, and B3,C4,Z1 is more favorable than Z0,B3,Z1 and Z0,C4,Z1. It can be verified that with or without using Z0,B3,Z1 and Z0,C4,Z1 as outer regions, there is no accuracy difference in the approximation. Hence, one significant difference between using a region-based approximation for BN and MN models is that the node ordering information of a BN cannot be ignored.
Given all CPDs are pairwise and singletons for the Figure4(b) BN (with the moral graph shown in Figure 5(b)), the joint distribution of the BN is rewritten as pZ0,…,F7=1/Z∏i,j∈EφijXi,Xj∏i∈VφiXi,i,j=0…7. The improved LSRG region graph for this BN is shown in Table 1.
In Table1, the region counting numbers are listed aside from each region. The inner regions are obtained using the CVM algorithm. There are no third-level regions because their counting numbers are zeros. It can be verified that this region graph satisfies the tree-robustness property because the subgraphs are all planar graphs.
Our improved LSRG algorithm follows directly with the PWC algorithm. So, the proposed PWC-LSRG algorithm is formally presented below in Algorithm1:
We omit the outer region size two scenarios in Algorithm1 because it is straightforward to select all node pairs connecting an edge in the moral graph as the outer regions. The associated energy function is the Bethe function.
The energy function defined by the PWC-LSRG region graph can be either optimized by the usual GBP/CCCP message-passing algorithms or amortized by training a (deep) neural network to receive performance gains. We show the amortized technique in the next section.Figure 5
(a) Complete BNG^π for Figure 4(b) BN G^, with bold dashed edges, a spanning tree of G^π. Each node is indexed with a number. Note that all the CPDs in G^π are pairwise and singleton factors. (b) Moral graphs of Figure 4(b)G^ with the dashed edge added as a chordal edge.
(a)(b)Table 1
Improved LSRG region graph for Figure4(b) BN.Algorithm 1: PWC-LSRG algorithm.
Input: a discrete BNG with CPD for each variable Xi(i=1,…,m)Output: region graphGR(1)
Perform BF (8) to convert G to its binary factorized BN G1, with variables Xj(j=1,…,m1);(2)
Initialize a BNG^=G1;(3)
for each CPDpXjXpaj in G^:(4)
ifpXjXpaj=3:(5)
add an intermediate nodeZ with the CPD defined in (9) and reconnect all nodes;(6)
convertpXpajXpapaj as hidden pairwise factors in G^;(7)
end if;(8)
end for;(9)
define a node orderingπ for all variables Xk(k=1,…,m2) in G^;(10)
convertG^ to its moralised graph MG^;(11)
find tree-robust subgraphsH of MG^; //construct tree-robust cycles(12)
Add simple cyclesB in H as outer regions to GR;(13)
Mark all edges inH as used;(14)
Mark all nodes inH as visited;(15)
for each unused edgee=s,t from a visited node s: //construct partial tree-robust cycles(16)
ift is unvisited:(17)
add an eare1 from s through edge e to some visited node u;(18)
find the shortest pathe2 from s to u through used edges;(19)
add cyclee1∪e2 as outer region to GR;(20)
end if;(21)
end for;(22)
for each size>3 outer region r(r∈R) in GR:(23)
add chordal edges to the cycle corresponding tor in MG^ by following the node ordering π to make the cycle chordal;(24)
replacer with its triplet decompositions r1,r2,…, in GR;(25)
end for;(26)
construct then-level region graph GR=CVMGR and assign all the pairwise/singleton factors into relevant outer regions;(27)
returnGR;


### 3.3. Amortizing

Amortizing [27–29] refers to minimizing the BFE or KFE objectives using (deep) neural networks (DNNs) and their associated loss functions. The parameters to optimize in the BFE/KFE are amortized by optimizing the neural network parameters. The input of the DNN are CPDs of the BN, and the loss function is defined by the BFE/KFE together with a penalty term. The output of the DNN is the predicted marginal beliefs composing the BEF/KFE objective function. The inference is carried out by training the DNN using forward and backward propagations to reconfigure the weights and bias parameters until convergence. The key to the DNN optimization is to define the penalty terms derived from the factor graph (for the BFE) or the region graph (for the KFE).It has been shown [11, 29] that the KFE is superior to the BFE, provided that the region graph is well defined, no matter whether the objective function is minimized by GBP message passing or amortized by neural networks. The BFE is also viewed as a special case of the KFE. So, we omit the BFE and its amortizing formulas [28] for brevity. Amortizing the KFE objective function involves the following steps:(1)
Given ann-level region graph GR, predict the region beliefs for outer regions brr∈Rlevel1 by the following:(10)brxr;ω=softmaxfGR,r;ω,where fGR,r;ω is a vector of scores for all variable configurations in region r, produced by a neural network with parameters ω. br is then predicted by using a SoftMax function of the neural network.(2)
Compute the inner region beliefsbuu∈Rlevel2:n by the following:(11)buxu;ω=1Rpau∑rp∈Rpau∑xp\xubpxp;ω,where Rpau denotes all the parent regions of u and ∙ is the length of the set. Equation (11) means that the inner region belief is the average over the same subset of all the parent region beliefs.(3)
The constrained KFE objective (as shown in (5)) is then minimized as a loss function:(12)minωFKFEbrxr+γ∑r∑rp∈Rpardbr,∑xp\xrbpxp;ω,where γ is a tuning parameter regularizing the belief consistency, and d∙,∙ is a L2 distance measuring the belief mismatch between the current and the estimated belief difference for determining convergence. Suppose the number of penalty terms in d∙,∙ is F, γ is closely related to F.The input of the DNN are pairwise factors converted from the BN CPDs, and the region graph constraints are used as penalties. The amortized technique can be viewed as consuming a region graph and outputting the beliefs as predictions.Specifically, each nodeXi has undergone a transformer [36] layer (the use of a transformer instead of conventional methods, such as LSTM or direct feedforward NN, is that we aim to explicitly model attention mechanism, both for extended use in learning tasks and gaining improved accuracy, i.e., [28] to conventional methods.) to obtain a hidden node representation hi, which is associated with a learnable embedding vector ei and shared with all layers. The outer regions ri,j,k∈R are of size three and their beliefs are obtained by applying an affine layer followed by a SoftMax function, i.e., bi,j,kxi,xj,xk=softmaxWhi;hj;hk+b. The inner region beliefs are calculated by (11). Consistency of the region’s beliefs is carried out by marginalizing irrelevant variables and comparing the difference between the shared variables. The DNN is continuously trained until the region beliefs’ mismatch is reduced to a certain threshold. The gradient of the loss function (12) is computed automatically by the torch [37] build-in functionality.Given all the input factors are pairwise and the proposed PWC-LSRG algorithm generates the region graph, we can directly apply the amortized inference framework (10)–(12) to the BN inference. The resulting amortized PWC-LSRG algorithm is presented in Algorithm 2.Algorithm 2: Amortized PWC-LSRG algorithm for BN inference.
Input: a Discrete BN G with CPD for each variable Xi(i=1,…,m)Output: marginal beliefs of G(1)
perform PWC-LSRG algorithm toG obtaining a pairwise BN G^ with variables Xj(j=1,…,m1), state dimension N, factors ϕk(k=1,…,m2), and its associated region graph GR with outer regions rs(s=1,…,m3);(2)
transform [36] Xj into NN node representation hj and initialize NN unary (m1×RNandouterregionm3×RN∙N∙N) embedding parameters;//attention layer(3)
mlp⟵ linear residual layers;//multilayer perception(4)
assign relevant factorsϕk to each outer region rs;(5)
assign uniform factors to each inner region;(6)
while region consistent error>ϵdo;(7)
obtain region logitslrs by mlp;(8)
obtain region beliefbr⟵softmaxlr;(9)
compute KFE (5) after forward propagation;(10)
compute region consistent error by comparing the parent and child region beliefs using (11);(11)
compute loss function (12) and perform backward propagation to update all parameters;(12)
end while;(13)
compute outer region marginal beliefs using (10);(14)
return the marginals ofXi in G by marginalizing irrelevant variables from the beliefs in #13Algorithm2 presents the algorithm when the outer region size is bounded to three, which is a balance between accuracy and efficiency. If we bound the outer region size to four, the region graph will extend to four levels, with all outer region sizes enforced to four. In our amortized framework, the DNN tensors must be equal sized thus outer regions of size <4 must be merged into size four. This process will introduce extra overheads. We can also lower the outer region size to two and change the #02 step of the algorithm to be m2×RN∙N embedding parameters. The outer regions are then composed of all node pairs connecting an edge in the moral graph of G^. The approximation is equivalent to amortized Bethe approximation, which is more efficient than using size three outer regions.The computational cost of the associated DNN is determined by both the number of nodes,m1 in BN G^, and the number of outer regions m3 generated by the PWC-LSRG. Each node is associated with a fixed state dimension N (i.e., 100). Given the number of transformer layers L is fixed, the node embeddings are of size L∙N∙m1. The region embeddings are then of size 3∙N∙m3, as it is a concatenation (in this case, triplets) of the node embeddings. The embeddings are then undergone the MLP and softmax layers. So, large BN with high tree width is likely more computationally intensive than low tree-width BNs, if there are more outer regions; hence, more penalty terms (m3) are produced. However, the tree-width does not necessarily determine the number of penalty terms because it only measures the length of the largest cluster. The number of penalty terms is determined by the cycles found by the PWC-LSRG algorithm (in Algorithm 1). We show tests in section 4.3 for empirical details.The flowchart of the amortized PWC-LSRG algorithm is shown in Figure6.Figure 6
Flowchart of the amortized PWC-LSRG algorithm.


## 3.1. Pairwise Decomposition of CPDs

As the CPD size lower bounds cluster size, ideally, the CPD size should be minimized thus the cluster size can also be reduced. The worst but also the most general case is the complete BN model in which all edges encode pairwise CPDs, so the outer region size can be reduced to two (which is the minimal size for an outer region). By Theorem1, any BN (including the partial BN) is a subset of the associated complete BN model. Suppose we can efficiently decompose a BN into all pairwise CPD representations and resolve its associated complete BN inference using region-based methods. In that case, we will resolve general BN inference problems efficiently.Reducing the CPD size can be achieved from either the factor level or the CPD level. Representing methods are binary factorization (BF) [15, 33] and MRF decomposition [34]. The BF algorithm reduces the CPD size by introducing intermediate variables to ensure each node connects no more than two parent nodes. The MRF decomposition [34] can decompose the general MN factors into pairwise factors by introducing an auxiliary node, with the new pairwise factors defined as the fraction power of the original factor. This method is inconvenient to apply to BN when the CPD size is large.Our proposed pairwise decomposition method is a two-step hybrid algorithm. We use BF to reduce the CPD size to three initially. Then, for each triplet CPD, we decompose it into three pairwise CPDs and two associated hidden pairwise factors. The process is exact. The joint distribution of the BN is then expressed as the product of all pairwise factors.We binary factorize the BNG by adding intermediate nodes and obtaining an equivalent BF model G′, in which every node has at the most two parent nodes. Equivalence means that the BF model G′ can always rebuild the CPDs in the BN G through intermediate nodes. It has been proved in [15] that the BF process is exact for discrete BNs, so any ordered joint distributions defined on the same variables between G and G′ are the same. The benefit of using the BF algorithm is to ensure the CPD size is no larger than three.Figure3 presents the BF algorithm for factorizing a 4-dimensional complete discrete BN G into its BF model G′. In general, the CPD for the added intermediate node E has the form:(8)pE=eijX=xp,Y=yq,=1,ifp=iandq=j,0,otherwise,where X and Y are parents of E and xp, (p=1,…,m1) and yq (q=1,…,m2) are node states, and eij (i=1,…,m1, j=1,…,m2) is an m1×m2 CPD matrix for E. Following (8) we can always obtain a binary factorized BN, with a maximum CPD size of three.Figure 3
(a) A 4-dimensional complete discrete BNG with all binary variables; (b) a binary factorized BN G′ of (a). The singleton marginal distributions are obtained using the JT algorithm in AgenaRisk [35].
(a)(b)In general, for a triplet CPDpX3X1,X2 with the associated possible pairwise CPDs pX1Xpa1 and pX2Xpa2, we can introduce an intermediate node Z with the settings in (9) to decompose pX3X1,X2 into pairwise CPDs. The scope of Z is the Cartesian product: SZ=SX1×SX2×SX3. Therefore, Z=zj (j=1,…,SZ) takes the value of all assignments of X1,X2,X3. We have the following:(9)pZ=φ123X1,X2,X3,pX1=x1Z=zj=1,ifx1⊂assignmentj,pX2=x2Z=zj=1,ifx2⊂assignmentj,pX3=x3Z=zj=1,ifx3⊂assignmentj,where φ123X1,X2,X3 is the function of pX3X1,X2. Due to the complementary property of the CPD, other entries not in (9) are set to zero. Because (9) adds new connections Z⟶X1 and Z⟶X2 for constructing the new CPDs, the original CPDs pX1Xpa1 and pX2Xpa2 cannot be explicitly represented in the BN. So, they are set as hidden pairwise factors (as opposed to CPDs) in the inference phase. Note that we cannot multiply the CPD pX1Xpa1 into pX1Z as that breaks the pairwise factor constraint.Figure4 shows an example using (9). The CPD pCA,B in Figure 4(a)G does not associate with pairwise CPDs for the parent nodes A and B, so we simply introduce a new node Z0 in 4(b)G^ following (9). Node C is then disconnected to A and B and connects to Z0 instead. Similarly, we need to introduce a new node Z1 to decompose the CPD pFD,E in 4(b). The CPD pFD,E is associated with two pairwise CPDs pDB and pEC, which needs to be set as hidden pairwise factors, shown as dashed edges in Figure 4(b). After exact inference, both 4(a)G and 4(b)G^ produced identical marginal results.Figure 4
(a) A binary factorized BNG; (b) pairwise converted BN G^ of G. The dashed edges represent hidden CPDs.
(a)(b)Proposition 1.
The pairwise conversion (PWC) (9) is an exact procedure to convert a binary factorized BN G into its equivalent pairwise representation BN G^.Proof.
The CPDpX1=x1Z=zj=1 in (9) can be viewed as an indicator function, i.e., taking the relevant values in φ123X1,X2,X3. So, multiplying all the CPDs in (9) and integrating Z out, we end up with φ123X1,X2,X3, equaling pX3X1,X2. The hidden pairwise factors are original BN CPDs. Therefore, our pairwise converted BN G^ does not change the chain rule (1) for the binary factorized BN G. Hence, pG^=pG.


## 3.2. LSRG Algorithm with Node Ordering Information

We can apply LSRG to the moral graph of a pairwise converted BNG^. However, the LSRG does not incorporate the node ordering information of the BN, i.e., ignoring the ordered Markov property that results in less optimal region graph constructions for BNs. To improve that, we introduce a node ordering π using the BN G^.Definition 2.
The node orderingπ is defined as π≔X1⟶,…,⟶Xm for the m variables in the pairwise BN G^, where for any Xj that is, the successor to Xi, Xj∉Xpai.
In complete BNs, the node orderingπ is fixed in the model. As introduced in section 2.2, Definition 1, the tree-robustness region graph can be obtained conveniently for a complete MN. And a complete MN has the same structure as a complete BN’s moral graph. So, the prescription of Definition 1 also works for complete BNs. For sparse BNs, we can use Definition 2 and Theorem 1, adding directed edges to the BN, obtaining an equivalent complete BN, and using Definition 1 to build a tree-robustness region graph.
For example, following Definition2, we define a node ordering for Figure 4(b)G^ with π≔Z0⟶A⟶B⟶C⟶Z1⟶D⟶E⟶F to obtain an 8-dimensional complete BN G^π as shown in Figure 5(a). The pairwise CPDs that are not contained in G^ encode uniform distributions. We can choose Z0 as the root node for LSRG to generate the outer regions. The resulting outer region size is three regardless of how large the tree-width of the BN is.
Although we can always use the complete BN representationG^π to obtain a tree-robustness region graph for a sparse BN G^, we are left with a dramatically increased number of outer regions along with the model dimensionality. This is because most of the outer regions are added unnecessarily. Also, all the outer regions contain the root node, resulting in possible numerical issues [15] associated with the region graph.
We will therefore generate the LSRG outer regions directly from the moral graph ofG^ instead of G^π, as shown in Figure 5(b). The resulting region graph will only ensure partial tree-robustness rather than tree-robustness. Following that, we obtain a significant reduction in outer regions while retaining sensible approximations. For example, the number of LSRG outer regions will be reduced from 21 to 6 when using the Figure 5(b) model rather than Figure 5(a) model.
After performing the PWC algorithm, all BN nodes will be connected to no more than two nodes. So, the subgraphs of the BN will be of simple forms, such as planar graphs and/or four-dimensional complete DAGs. It is impossible to contain a complete DAG larger than four dimensions in the subgraph because it violates the PWC algorithm.
If we bound the maximum region size to two, we can use the node pairs connecting an edge in the moral graph as outer regions. In this case, the Kikuchi approximation equals the Bethe approximation. If we bound the maximum region size to three, we will obtain LSRG outer regions with triplet cycles and cycles of length>3. The outer regions found by LSRG are not equal sized. So, the length >3 outer regions will need further decompositions. This is one limitation that the LSRG does not control the outer region sizes.
Our improved LSRG algorithm contains two more steps: (i) add chordal edges to cycles with length>3 to make it chordal. (ii) The chordal edges are added in a sequence following the node ordering π.
For example, in Figure5(b), there is a length four cycle face Z0,B3,Z1,C4 found by LSRG. Following the node ordering π, a chordal edge being added is B3,C4 rather than Z0,Z1 because Z1 is ordered after B3 and C4. The chordal edge connecting B3 and C4 has higher priority than the edge connecting Z0 and Z1. This is one instance of the ordered Markov property that the BN holds for conditional independence Z1⊥Z0pB3,C4. So, the set of outer regions is Z0,B3,C4, and B3,C4,Z1 is more favorable than Z0,B3,Z1 and Z0,C4,Z1. It can be verified that with or without using Z0,B3,Z1 and Z0,C4,Z1 as outer regions, there is no accuracy difference in the approximation. Hence, one significant difference between using a region-based approximation for BN and MN models is that the node ordering information of a BN cannot be ignored.
Given all CPDs are pairwise and singletons for the Figure4(b) BN (with the moral graph shown in Figure 5(b)), the joint distribution of the BN is rewritten as pZ0,…,F7=1/Z∏i,j∈EφijXi,Xj∏i∈VφiXi,i,j=0…7. The improved LSRG region graph for this BN is shown in Table 1.
In Table1, the region counting numbers are listed aside from each region. The inner regions are obtained using the CVM algorithm. There are no third-level regions because their counting numbers are zeros. It can be verified that this region graph satisfies the tree-robustness property because the subgraphs are all planar graphs.
Our improved LSRG algorithm follows directly with the PWC algorithm. So, the proposed PWC-LSRG algorithm is formally presented below in Algorithm1:
We omit the outer region size two scenarios in Algorithm1 because it is straightforward to select all node pairs connecting an edge in the moral graph as the outer regions. The associated energy function is the Bethe function.
The energy function defined by the PWC-LSRG region graph can be either optimized by the usual GBP/CCCP message-passing algorithms or amortized by training a (deep) neural network to receive performance gains. We show the amortized technique in the next section.Figure 5
(a) Complete BNG^π for Figure 4(b) BN G^, with bold dashed edges, a spanning tree of G^π. Each node is indexed with a number. Note that all the CPDs in G^π are pairwise and singleton factors. (b) Moral graphs of Figure 4(b)G^ with the dashed edge added as a chordal edge.
(a)(b)Table 1
Improved LSRG region graph for Figure4(b) BN.Algorithm 1: PWC-LSRG algorithm.
Input: a discrete BNG with CPD for each variable Xi(i=1,…,m)Output: region graphGR(1)
Perform BF (8) to convert G to its binary factorized BN G1, with variables Xj(j=1,…,m1);(2)
Initialize a BNG^=G1;(3)
for each CPDpXjXpaj in G^:(4)
ifpXjXpaj=3:(5)
add an intermediate nodeZ with the CPD defined in (9) and reconnect all nodes;(6)
convertpXpajXpapaj as hidden pairwise factors in G^;(7)
end if;(8)
end for;(9)
define a node orderingπ for all variables Xk(k=1,…,m2) in G^;(10)
convertG^ to its moralised graph MG^;(11)
find tree-robust subgraphsH of MG^; //construct tree-robust cycles(12)
Add simple cyclesB in H as outer regions to GR;(13)
Mark all edges inH as used;(14)
Mark all nodes inH as visited;(15)
for each unused edgee=s,t from a visited node s: //construct partial tree-robust cycles(16)
ift is unvisited:(17)
add an eare1 from s through edge e to some visited node u;(18)
find the shortest pathe2 from s to u through used edges;(19)
add cyclee1∪e2 as outer region to GR;(20)
end if;(21)
end for;(22)
for each size>3 outer region r(r∈R) in GR:(23)
add chordal edges to the cycle corresponding tor in MG^ by following the node ordering π to make the cycle chordal;(24)
replacer with its triplet decompositions r1,r2,…, in GR;(25)
end for;(26)
construct then-level region graph GR=CVMGR and assign all the pairwise/singleton factors into relevant outer regions;(27)
returnGR;


## 3.3. Amortizing

Amortizing [27–29] refers to minimizing the BFE or KFE objectives using (deep) neural networks (DNNs) and their associated loss functions. The parameters to optimize in the BFE/KFE are amortized by optimizing the neural network parameters. The input of the DNN are CPDs of the BN, and the loss function is defined by the BFE/KFE together with a penalty term. The output of the DNN is the predicted marginal beliefs composing the BEF/KFE objective function. The inference is carried out by training the DNN using forward and backward propagations to reconfigure the weights and bias parameters until convergence. The key to the DNN optimization is to define the penalty terms derived from the factor graph (for the BFE) or the region graph (for the KFE).It has been shown [11, 29] that the KFE is superior to the BFE, provided that the region graph is well defined, no matter whether the objective function is minimized by GBP message passing or amortized by neural networks. The BFE is also viewed as a special case of the KFE. So, we omit the BFE and its amortizing formulas [28] for brevity. Amortizing the KFE objective function involves the following steps:(1)
Given ann-level region graph GR, predict the region beliefs for outer regions brr∈Rlevel1 by the following:(10)brxr;ω=softmaxfGR,r;ω,where fGR,r;ω is a vector of scores for all variable configurations in region r, produced by a neural network with parameters ω. br is then predicted by using a SoftMax function of the neural network.(2)
Compute the inner region beliefsbuu∈Rlevel2:n by the following:(11)buxu;ω=1Rpau∑rp∈Rpau∑xp\xubpxp;ω,where Rpau denotes all the parent regions of u and ∙ is the length of the set. Equation (11) means that the inner region belief is the average over the same subset of all the parent region beliefs.(3)
The constrained KFE objective (as shown in (5)) is then minimized as a loss function:(12)minωFKFEbrxr+γ∑r∑rp∈Rpardbr,∑xp\xrbpxp;ω,where γ is a tuning parameter regularizing the belief consistency, and d∙,∙ is a L2 distance measuring the belief mismatch between the current and the estimated belief difference for determining convergence. Suppose the number of penalty terms in d∙,∙ is F, γ is closely related to F.The input of the DNN are pairwise factors converted from the BN CPDs, and the region graph constraints are used as penalties. The amortized technique can be viewed as consuming a region graph and outputting the beliefs as predictions.Specifically, each nodeXi has undergone a transformer [36] layer (the use of a transformer instead of conventional methods, such as LSTM or direct feedforward NN, is that we aim to explicitly model attention mechanism, both for extended use in learning tasks and gaining improved accuracy, i.e., [28] to conventional methods.) to obtain a hidden node representation hi, which is associated with a learnable embedding vector ei and shared with all layers. The outer regions ri,j,k∈R are of size three and their beliefs are obtained by applying an affine layer followed by a SoftMax function, i.e., bi,j,kxi,xj,xk=softmaxWhi;hj;hk+b. The inner region beliefs are calculated by (11). Consistency of the region’s beliefs is carried out by marginalizing irrelevant variables and comparing the difference between the shared variables. The DNN is continuously trained until the region beliefs’ mismatch is reduced to a certain threshold. The gradient of the loss function (12) is computed automatically by the torch [37] build-in functionality.Given all the input factors are pairwise and the proposed PWC-LSRG algorithm generates the region graph, we can directly apply the amortized inference framework (10)–(12) to the BN inference. The resulting amortized PWC-LSRG algorithm is presented in Algorithm 2.Algorithm 2: Amortized PWC-LSRG algorithm for BN inference.
Input: a Discrete BN G with CPD for each variable Xi(i=1,…,m)Output: marginal beliefs of G(1)
perform PWC-LSRG algorithm toG obtaining a pairwise BN G^ with variables Xj(j=1,…,m1), state dimension N, factors ϕk(k=1,…,m2), and its associated region graph GR with outer regions rs(s=1,…,m3);(2)
transform [36] Xj into NN node representation hj and initialize NN unary (m1×RNandouterregionm3×RN∙N∙N) embedding parameters;//attention layer(3)
mlp⟵ linear residual layers;//multilayer perception(4)
assign relevant factorsϕk to each outer region rs;(5)
assign uniform factors to each inner region;(6)
while region consistent error>ϵdo;(7)
obtain region logitslrs by mlp;(8)
obtain region beliefbr⟵softmaxlr;(9)
compute KFE (5) after forward propagation;(10)
compute region consistent error by comparing the parent and child region beliefs using (11);(11)
compute loss function (12) and perform backward propagation to update all parameters;(12)
end while;(13)
compute outer region marginal beliefs using (10);(14)
return the marginals ofXi in G by marginalizing irrelevant variables from the beliefs in #13Algorithm2 presents the algorithm when the outer region size is bounded to three, which is a balance between accuracy and efficiency. If we bound the outer region size to four, the region graph will extend to four levels, with all outer region sizes enforced to four. In our amortized framework, the DNN tensors must be equal sized thus outer regions of size <4 must be merged into size four. This process will introduce extra overheads. We can also lower the outer region size to two and change the #02 step of the algorithm to be m2×RN∙N embedding parameters. The outer regions are then composed of all node pairs connecting an edge in the moral graph of G^. The approximation is equivalent to amortized Bethe approximation, which is more efficient than using size three outer regions.The computational cost of the associated DNN is determined by both the number of nodes,m1 in BN G^, and the number of outer regions m3 generated by the PWC-LSRG. Each node is associated with a fixed state dimension N (i.e., 100). Given the number of transformer layers L is fixed, the node embeddings are of size L∙N∙m1. The region embeddings are then of size 3∙N∙m3, as it is a concatenation (in this case, triplets) of the node embeddings. The embeddings are then undergone the MLP and softmax layers. So, large BN with high tree width is likely more computationally intensive than low tree-width BNs, if there are more outer regions; hence, more penalty terms (m3) are produced. However, the tree-width does not necessarily determine the number of penalty terms because it only measures the length of the largest cluster. The number of penalty terms is determined by the cycles found by the PWC-LSRG algorithm (in Algorithm 1). We show tests in section 4.3 for empirical details.The flowchart of the amortized PWC-LSRG algorithm is shown in Figure6.Figure 6
Flowchart of the amortized PWC-LSRG algorithm.


## 4. Experiments

We validate the amortized PWC-LSRG algorithm and its subalgorithms separately. The PWC and the PWC-LSRG region graph are validated in sections4.1 and 4.2. Then, the amortized PWC-LSRG is tested in section 4.3.


### 4.1. Validation of the PWC Algorithm

The testing environment is Java 1.8, i5-9400H. We use Figure4 BN, the well-known Student [1], and Grid [38] BNs for the testing. All these tests can be accessed online [39]. We index all nodes 1,…,m in the BN using the parent-child ordering and associate observations for a root parent node i and the child node m. The observation setting is the same for all tests unless noted. All the CPDs are randomly generated. The testing purpose is to validate that the PWC algorithm is exact.We use the probability of evidencelogpe [1] for the testing as logpe reflects the overall accuracy of the BN when entering observations to the variables. Table 2 shows that the results obtained by using BN G^ are the same as those obtained by G under the JT inference. These tests also empirically validated Proposition 1.Table 2
logpe comparison of G and its pairwise converted BN G^ using JT inference.
BNslogpe of Glogpe of G^Abs. diff(1) Figure4 BN (m = 6, t.w. 3)−1.4579−1.45790(2) Student (m = 8, t.w. 3)−1.6478−1.64780(3 )3 × 3 grid (m = 9, t.w. 3)−1.8457−1.84570“m” represents the number of variables in the BN; “t.w.” is short for tree-width; “Abs. diff.” is short for the absolute difference.


### 4.2. Validation of the PWC-LSRG Region Graph

The testing environment is Java 1.8, i5-9400H. We use the BNs in section4.1 and two more well-known BNs, Asia [7] and HMM [1], for the testing. The testing purpose is to validate that the regional graph generated by PWC-LSRG yields sensible approximations under CCCP [17] (convergent alternative of GBP) message passing.Table3 shows the exact results obtained by JT and the absolute difference for all approximate algorithms. The PWC-LSRG outperforms other approximate algorithms on accuracy as BP and IJGP degrade significantly when the tree-width increases. For example, they are considerably worse than PWC-LSRG on the 5×5 grid, shown with results in orange. Because PWC-LSRG improves over the LSRG, we omit the results compared between them.Table 3
logpe comparison of PWC-LSRG with competing algorithms.
BNsJT (exact)BP (diff.)IJGP (diff.)PWC-LSRG (diff.)1. Figure4(b) BN (m = 8, t.w. 3)−1.1112.90E − 072.90E − 072.90E − 072. Asia (m = 8, t.w. 3)−0.9592.01E − 1303.69E − 123. Student (m = 8, t.w. 3)−0.5909.91E − 112.19E − 112.64E-104.2 × 3 HMM (m = 6, t.w. 3)−1.4107.58E − 101.39E − 099.72E − 115.5 × 5 grid (m = 25, t.w. 6)−0.3200.0027920.0030361.50E − 04“m” represents the number of variables in the BN; “t.w.” is short for tree-width; “diff.” is short for the absolute difference. Bold values are the smallest and the best results. Italic values are the worst results.


### 4.3. Validation of the Amortized PWC-LSRG Algorithm

Then, the testing environment is Python 3.9, Apple Silicon M1-Max. All test BNs are grouped into low, median, and high by tree-width, as shown in Table4 (model accessible at [39]). Except noted, all test BNs contain binary variables with CPDs generated from uniform0,1. The Barley [5], Diabetes [40], and Pedigree pigs [5] are used in practice for medical and biology analysis. The testing purpose is to validate the efficiency gains obtained by using amortized PWC-LSRG algorithm (code accessible at [39]).Table 4
logpe comparison of amortized PWC-LSRG with competing algorithms.
BNsExactGBP-MP (diff.)Amo. Bethe (diff.)Amo. PWC-LSRG (diff.)Low t.w(1) Figure3(b) BN (m = 5, t.w. 3)−3.060.620.210.16(2) Figure4(b) BN (m = 8, t.w. 3)−2.440.120.710.53(3) Asia (m = 8, t.w. 3)−2.960.090.570.23(4) Student (m = 8, t.w. 3)−3.690.470.280.11Median t.w(5) 3 × 3 DBN (m = 9, t.w. 4)−6.111.472.380.45(6) 5 × 5 grid (m = 25, t.w. 6)−8.733.100.301.18(7) Barley (m = 58, t.w. 7)−19.471.642.273.41High t.w(8) Diabetes (m = 413, t.w. 12)−128.6582.748.169.53(9) Pedigree (m = 441, t.w. 40)−143.0626.6538.7520.17“m” represents the number of variables in the BN; “t.w.” is short for tree-width; “diff.” is short for the absolute difference. Bold values are the smallest and the best results. Italic values are the worst results in their group (due to the large tree-widthof #9 Pedigree BN the exact result cannot be obtained by JT. We instead use the maximum cluster size of 20 to compute the near exact result).The competing algorithms are amortized Bethe approximation [28], GBP with parent-to-child region message passing [11] (it is one of the fastest message-passing themes), and the amortized PWC-LSRG. For fairness, we use the same outer regions generated by PWC-LSRG for the GBP-MP algorithm because GBP does not automatically define outer regions. Also, both amortized Bethe and amortized PWC-LSRG use the same (deep) neural network layers (attention + two residual layers + softmax) for fairness.The different tree-width BNs test results are shown in Figure7. In Figure 8, we perform the worst case-complete BN tests. Finally, in Figure 9, we test the tuning parameter γ variations for different-sized BNs to discuss its effects on accuracy.Figure 7
2 × 2 figure box comparing accuracy on (a) low t.w. BNs and (c) median and high t.w. BNs, using singleton marginal beliefs; comparing efficiency on (b) low t.w. BNs and (d) median and high t.w. BNs in seconds. The competing algorithms are amortized Bethe (short for Amo. Bethe) [28], GBP-MP [11], and amortized PWC-LSRG (short for Amo. PWC-LSRG). (a) Marginal belief comparison. (b) Efficiency comparison with normal factors. (c) Marginal belief comparison. (d) Efficiency comparison with normal factors.
(a)(b)(c)(d)Figure 8
2 × 2 figure box comparing the accuracy (a) and (b) and efficiency (c) and (d) for the competing algorithms: mean-field approximation (short for M-F), amortized Bethe (short for Amo. Bethe) [28], GBP-MP [11], and amortized PWC-LSRG (short for Amo. PWC-LSRG). (a) log PE comparison. (b) Marginal belief comparison. (c) Efficiency comparison with normal factors. (d) Efficiency comparison with extreme factors.
(a)(b)(c)(d)Figure 9
Tuning parameterγ variants for the low, median, and high t.w. BNs on (a) logpe approximation and (b) marginal belief approximation. Optimal γ against the number of penalty terms F on (c) logpe approximation and (d) marginal belief (short for MAR) approximation. (a) log PE comparison. (b) Marginal belief comparison. (c) No of penalties vs tuning parameter on log PE. (d) No of penalties vs tuning parameter on MAR.
(a)(b)(c)(d)In Figure7(a), we set the tuning parameter γ=5 for these low tree-width BNs. Details of the settings for γ will be discussed in Figure 9. Both amortized Bethe and amortized PWC-LSRG performed better than GBP-MP on marginal beliefs accuracy. The GBP-MP has some variables poorly approximated, i.e., KL > 0.1, so the statistics are destroyed. However, GBP-MP achieved the most efficiency for the low tree-width BNs, as shown in Figure 7(b), given the number of variables is small in these tests. Still, all algorithms complete the tests within 1 second. The tuning parameters for the median and high tree-width BNs are set as γ=10 and γ=50, respectively. Although GBP-MP uses the same outer regions as Amo. PWC-LSRG, the inaccuracy presented by GBP-MP is even more evident in Figure 7(c) on median and high tree-width BNs, indicating that the variant of message-passing themes is unstable. The variances of the KL errors on the Barely and Diabetes BNs for GBP-MP are not shown in Figure 7(c) since they are destroyed in the statistics. The Amo. PWC-LSRG is consistently stable in the high tree-width BN tests and compares more favorably than Amo. Bethe on accuracy. The average KL error of the Amo. PWC-LSRG is of 1e − 2 magnitude, indicating its efficacy for practical use. The efficiency gain is significant in using the Amo PWC-LSRG, as shown in Figure 7(d). For example, the Amo PWC-LSRG is almost five times faster compared to the GBP-MP (80 seconds vs. 400 seconds) on the Diabetes BN test, which is a large BN yielding 874 outer regions.Interestingly, the Amo PWC-LSRG computes the Pedigree (m = 441, t.w. 40) BN more efficiently than the Diabetes (m = 413, t.w. 12) BN, although the Pedigree has more nodes and higher t.w. than Diabetes. This is because the penalty terms used in Diabetes BN is twice (874 vs. 384) the number of the Pedigree BN. So, the computational cost is dominated by the number of penalty terms in the amortized Kikuchi function.We also test the Amo. PWC-LSRG on Nvidia Tesla v100 GPU on the Diabetes BN, which completes the task in 120 seconds. M1-Max is the latest high-speed CPU that achieved better efficiency than the GPU for large BNs. But in either platform, the efficiency gain is almost four times of the message-passing algorithm.The correspondinglogpe comparison is shown in Table 4. The discrepancies of the logpe results among all algorithms increase along with the tree-width. In high tree-width BN tests, both GBP-MP and Amo. Bethe present more significant errors than Amo. PWC-LSRG, as shown in orange values.Then, we test on the worst-case complete BNs. We randomly generate pairwise CPDs for complete BNs from 4 to 16 variables without associating any observations. Thus, the tree-width of the BNs ranges from 3 to 15. Although the number of variables is not large, in practice, most BNs will not exceed the tree-width of 15, and the structures are less dense than complete BNs. Thus, the test is still challenging. Except for the BP class of algorithms, we also include the mean-field [9] approximation as one representative of the VI class of algorithms for comparison.Figure8 error bar results of mean and std. are obtained using 20 instances for each subfigure. It is expected that Amo. PWC-LSRG would outperform the competing algorithms in accuracy, as shown in Figures 8(a) and 8(b) for both logpe and marginal beliefs. The Amo. PWC-LSRG is about four times more efficient than GBP-MP using both normal (uniform0,1) and extreme (near zero and one) CPDs, as shown in Figures 8(c) and 8(d). The Amo. Bethe is slightly more efficient than Amo. PWC-LSRG in extreme CPD cases, since it consumes a region graph with outer regions size two rather than three. The mean-field approximation (M-F) is the most efficient algorithm among the competing algorithms and also performed remarkably on the logpe approximation. However, it is significantly worse than Amo. PWC-LSRG on the marginal belief approximation.Finally, we test the tuning parameterγ variants against the number of penalty terms F (shown in (12)) associated with the Kikuchi energy function on both logpe and marginal beliefs. The testing aims to help determine the optimal value of γ for different-sized BNs. The number of penalty terms F can be roughly measured by the overall number of regions in the region graph. In Figures 9(a) and 9(b), we use the Student, Barley, and Diabetes BNs as representatives for the low, median, and high tree-width BNs. The tuning parameter γ is set discrete values in 3,5,10,20,50,100 increasingly. The optimal value of γ for the Student BN is 3 (or 5) for both logpe and marginal belief approximations. But for the Barley and Diabetes BNs, the optimal γ disagrees on logpe and marginal belief approximations, and tends to increase along with the scale of the BN.In Figures9(c) and 9(d), we run all test BNs (in Table 4) multiple times and record the optimal value of γ for both logpe and marginal belief approximations. We then against the optimal γ to the number of penalty terms F and fit the curves. On the logpe approximation, as shown in (c), the optimal γ is proportional to lnF, and the Pearson correlation is 0.78. In contrast, the optimal γ is linear to F on the marginal belief approximation, with a high Pearson correlation of 0.98. Therefore, it is challenging to know the optimal γ for different tasks. However, based on the above analysis, we can set an optimal relative value of γ=5 for low t.w. small-sized BNs and γ=50 for high t.w. large-sized BNs in the marginal belief approximations. Or double these values of γ for their corresponding logpe approximation. Also, our findings of γ against F share some nonlinearity effects with the conclusions from Amo. Bethe algorithm [28], in which γ scales with 1/F2.


## 4.1. Validation of the PWC Algorithm

The testing environment is Java 1.8, i5-9400H. We use Figure4 BN, the well-known Student [1], and Grid [38] BNs for the testing. All these tests can be accessed online [39]. We index all nodes 1,…,m in the BN using the parent-child ordering and associate observations for a root parent node i and the child node m. The observation setting is the same for all tests unless noted. All the CPDs are randomly generated. The testing purpose is to validate that the PWC algorithm is exact.We use the probability of evidencelogpe [1] for the testing as logpe reflects the overall accuracy of the BN when entering observations to the variables. Table 2 shows that the results obtained by using BN G^ are the same as those obtained by G under the JT inference. These tests also empirically validated Proposition 1.Table 2
logpe comparison of G and its pairwise converted BN G^ using JT inference.
BNslogpe of Glogpe of G^Abs. diff(1) Figure4 BN (m = 6, t.w. 3)−1.4579−1.45790(2) Student (m = 8, t.w. 3)−1.6478−1.64780(3 )3 × 3 grid (m = 9, t.w. 3)−1.8457−1.84570“m” represents the number of variables in the BN; “t.w.” is short for tree-width; “Abs. diff.” is short for the absolute difference.


## 4.2. Validation of the PWC-LSRG Region Graph

The testing environment is Java 1.8, i5-9400H. We use the BNs in section4.1 and two more well-known BNs, Asia [7] and HMM [1], for the testing. The testing purpose is to validate that the regional graph generated by PWC-LSRG yields sensible approximations under CCCP [17] (convergent alternative of GBP) message passing.Table3 shows the exact results obtained by JT and the absolute difference for all approximate algorithms. The PWC-LSRG outperforms other approximate algorithms on accuracy as BP and IJGP degrade significantly when the tree-width increases. For example, they are considerably worse than PWC-LSRG on the 5×5 grid, shown with results in orange. Because PWC-LSRG improves over the LSRG, we omit the results compared between them.Table 3
logpe comparison of PWC-LSRG with competing algorithms.
BNsJT (exact)BP (diff.)IJGP (diff.)PWC-LSRG (diff.)1. Figure4(b) BN (m = 8, t.w. 3)−1.1112.90E − 072.90E − 072.90E − 072. Asia (m = 8, t.w. 3)−0.9592.01E − 1303.69E − 123. Student (m = 8, t.w. 3)−0.5909.91E − 112.19E − 112.64E-104.2 × 3 HMM (m = 6, t.w. 3)−1.4107.58E − 101.39E − 099.72E − 115.5 × 5 grid (m = 25, t.w. 6)−0.3200.0027920.0030361.50E − 04“m” represents the number of variables in the BN; “t.w.” is short for tree-width; “diff.” is short for the absolute difference. Bold values are the smallest and the best results. Italic values are the worst results.


## 4.3. Validation of the Amortized PWC-LSRG Algorithm

Then, the testing environment is Python 3.9, Apple Silicon M1-Max. All test BNs are grouped into low, median, and high by tree-width, as shown in Table4 (model accessible at [39]). Except noted, all test BNs contain binary variables with CPDs generated from uniform0,1. The Barley [5], Diabetes [40], and Pedigree pigs [5] are used in practice for medical and biology analysis. The testing purpose is to validate the efficiency gains obtained by using amortized PWC-LSRG algorithm (code accessible at [39]).Table 4
logpe comparison of amortized PWC-LSRG with competing algorithms.
BNsExactGBP-MP (diff.)Amo. Bethe (diff.)Amo. PWC-LSRG (diff.)Low t.w(1) Figure3(b) BN (m = 5, t.w. 3)−3.060.620.210.16(2) Figure4(b) BN (m = 8, t.w. 3)−2.440.120.710.53(3) Asia (m = 8, t.w. 3)−2.960.090.570.23(4) Student (m = 8, t.w. 3)−3.690.470.280.11Median t.w(5) 3 × 3 DBN (m = 9, t.w. 4)−6.111.472.380.45(6) 5 × 5 grid (m = 25, t.w. 6)−8.733.100.301.18(7) Barley (m = 58, t.w. 7)−19.471.642.273.41High t.w(8) Diabetes (m = 413, t.w. 12)−128.6582.748.169.53(9) Pedigree (m = 441, t.w. 40)−143.0626.6538.7520.17“m” represents the number of variables in the BN; “t.w.” is short for tree-width; “diff.” is short for the absolute difference. Bold values are the smallest and the best results. Italic values are the worst results in their group (due to the large tree-widthof #9 Pedigree BN the exact result cannot be obtained by JT. We instead use the maximum cluster size of 20 to compute the near exact result).The competing algorithms are amortized Bethe approximation [28], GBP with parent-to-child region message passing [11] (it is one of the fastest message-passing themes), and the amortized PWC-LSRG. For fairness, we use the same outer regions generated by PWC-LSRG for the GBP-MP algorithm because GBP does not automatically define outer regions. Also, both amortized Bethe and amortized PWC-LSRG use the same (deep) neural network layers (attention + two residual layers + softmax) for fairness.The different tree-width BNs test results are shown in Figure7. In Figure 8, we perform the worst case-complete BN tests. Finally, in Figure 9, we test the tuning parameter γ variations for different-sized BNs to discuss its effects on accuracy.Figure 7
2 × 2 figure box comparing accuracy on (a) low t.w. BNs and (c) median and high t.w. BNs, using singleton marginal beliefs; comparing efficiency on (b) low t.w. BNs and (d) median and high t.w. BNs in seconds. The competing algorithms are amortized Bethe (short for Amo. Bethe) [28], GBP-MP [11], and amortized PWC-LSRG (short for Amo. PWC-LSRG). (a) Marginal belief comparison. (b) Efficiency comparison with normal factors. (c) Marginal belief comparison. (d) Efficiency comparison with normal factors.
(a)(b)(c)(d)Figure 8
2 × 2 figure box comparing the accuracy (a) and (b) and efficiency (c) and (d) for the competing algorithms: mean-field approximation (short for M-F), amortized Bethe (short for Amo. Bethe) [28], GBP-MP [11], and amortized PWC-LSRG (short for Amo. PWC-LSRG). (a) log PE comparison. (b) Marginal belief comparison. (c) Efficiency comparison with normal factors. (d) Efficiency comparison with extreme factors.
(a)(b)(c)(d)Figure 9
Tuning parameterγ variants for the low, median, and high t.w. BNs on (a) logpe approximation and (b) marginal belief approximation. Optimal γ against the number of penalty terms F on (c) logpe approximation and (d) marginal belief (short for MAR) approximation. (a) log PE comparison. (b) Marginal belief comparison. (c) No of penalties vs tuning parameter on log PE. (d) No of penalties vs tuning parameter on MAR.
(a)(b)(c)(d)In Figure7(a), we set the tuning parameter γ=5 for these low tree-width BNs. Details of the settings for γ will be discussed in Figure 9. Both amortized Bethe and amortized PWC-LSRG performed better than GBP-MP on marginal beliefs accuracy. The GBP-MP has some variables poorly approximated, i.e., KL > 0.1, so the statistics are destroyed. However, GBP-MP achieved the most efficiency for the low tree-width BNs, as shown in Figure 7(b), given the number of variables is small in these tests. Still, all algorithms complete the tests within 1 second. The tuning parameters for the median and high tree-width BNs are set as γ=10 and γ=50, respectively. Although GBP-MP uses the same outer regions as Amo. PWC-LSRG, the inaccuracy presented by GBP-MP is even more evident in Figure 7(c) on median and high tree-width BNs, indicating that the variant of message-passing themes is unstable. The variances of the KL errors on the Barely and Diabetes BNs for GBP-MP are not shown in Figure 7(c) since they are destroyed in the statistics. The Amo. PWC-LSRG is consistently stable in the high tree-width BN tests and compares more favorably than Amo. Bethe on accuracy. The average KL error of the Amo. PWC-LSRG is of 1e − 2 magnitude, indicating its efficacy for practical use. The efficiency gain is significant in using the Amo PWC-LSRG, as shown in Figure 7(d). For example, the Amo PWC-LSRG is almost five times faster compared to the GBP-MP (80 seconds vs. 400 seconds) on the Diabetes BN test, which is a large BN yielding 874 outer regions.Interestingly, the Amo PWC-LSRG computes the Pedigree (m = 441, t.w. 40) BN more efficiently than the Diabetes (m = 413, t.w. 12) BN, although the Pedigree has more nodes and higher t.w. than Diabetes. This is because the penalty terms used in Diabetes BN is twice (874 vs. 384) the number of the Pedigree BN. So, the computational cost is dominated by the number of penalty terms in the amortized Kikuchi function.We also test the Amo. PWC-LSRG on Nvidia Tesla v100 GPU on the Diabetes BN, which completes the task in 120 seconds. M1-Max is the latest high-speed CPU that achieved better efficiency than the GPU for large BNs. But in either platform, the efficiency gain is almost four times of the message-passing algorithm.The correspondinglogpe comparison is shown in Table 4. The discrepancies of the logpe results among all algorithms increase along with the tree-width. In high tree-width BN tests, both GBP-MP and Amo. Bethe present more significant errors than Amo. PWC-LSRG, as shown in orange values.Then, we test on the worst-case complete BNs. We randomly generate pairwise CPDs for complete BNs from 4 to 16 variables without associating any observations. Thus, the tree-width of the BNs ranges from 3 to 15. Although the number of variables is not large, in practice, most BNs will not exceed the tree-width of 15, and the structures are less dense than complete BNs. Thus, the test is still challenging. Except for the BP class of algorithms, we also include the mean-field [9] approximation as one representative of the VI class of algorithms for comparison.Figure8 error bar results of mean and std. are obtained using 20 instances for each subfigure. It is expected that Amo. PWC-LSRG would outperform the competing algorithms in accuracy, as shown in Figures 8(a) and 8(b) for both logpe and marginal beliefs. The Amo. PWC-LSRG is about four times more efficient than GBP-MP using both normal (uniform0,1) and extreme (near zero and one) CPDs, as shown in Figures 8(c) and 8(d). The Amo. Bethe is slightly more efficient than Amo. PWC-LSRG in extreme CPD cases, since it consumes a region graph with outer regions size two rather than three. The mean-field approximation (M-F) is the most efficient algorithm among the competing algorithms and also performed remarkably on the logpe approximation. However, it is significantly worse than Amo. PWC-LSRG on the marginal belief approximation.Finally, we test the tuning parameterγ variants against the number of penalty terms F (shown in (12)) associated with the Kikuchi energy function on both logpe and marginal beliefs. The testing aims to help determine the optimal value of γ for different-sized BNs. The number of penalty terms F can be roughly measured by the overall number of regions in the region graph. In Figures 9(a) and 9(b), we use the Student, Barley, and Diabetes BNs as representatives for the low, median, and high tree-width BNs. The tuning parameter γ is set discrete values in 3,5,10,20,50,100 increasingly. The optimal value of γ for the Student BN is 3 (or 5) for both logpe and marginal belief approximations. But for the Barley and Diabetes BNs, the optimal γ disagrees on logpe and marginal belief approximations, and tends to increase along with the scale of the BN.In Figures9(c) and 9(d), we run all test BNs (in Table 4) multiple times and record the optimal value of γ for both logpe and marginal belief approximations. We then against the optimal γ to the number of penalty terms F and fit the curves. On the logpe approximation, as shown in (c), the optimal γ is proportional to lnF, and the Pearson correlation is 0.78. In contrast, the optimal γ is linear to F on the marginal belief approximation, with a high Pearson correlation of 0.98. Therefore, it is challenging to know the optimal γ for different tasks. However, based on the above analysis, we can set an optimal relative value of γ=5 for low t.w. small-sized BNs and γ=50 for high t.w. large-sized BNs in the marginal belief approximations. Or double these values of γ for their corresponding logpe approximation. Also, our findings of γ against F share some nonlinearity effects with the conclusions from Amo. Bethe algorithm [28], in which γ scales with 1/F2.


## 5. Conclusion and Future Works

We empirically validated that our proposed amortized PWC-LSRG algorithm is significantly more efficient than the conventional algorithms. Its accuracy is also of practical use and compares more favorably to the related amortized Bethe algorithm. More importantly, it avoids the iterative message-passing and sampling mechanisms, which are typical procedures used in other methods for carrying out BN inferences. The subalgorithm PWC-LSRG is also more advanced than other structural information-based region methods, which contributes to improving the accuracy. Future work will involve using the amortized PWC-LSRG framework for learning from data and models presenting latent variables. We will further study the optimal tuning parameter with different DNN layers and settings.



---

*Source: 2131915-2023-03-01.xml*