# SOLBP: Second-Order Loopy Belief Propagation for Inference in Uncertain
Bayesian Networks ††thanks: This research was partially supported by the Army
Research Office under contract W911NF-15-1-0479 and partially supported by the
U.S. Army Research Laboratory and the U.K. Ministry of Defence under agreement
W911NF-16-3-0001. The views and conclusions contained in this document are
those of the authors and should not be interpreted as representing the
official policies, either expressed or implied, of the U.S. Army Research
Laboratory, the U.S. Government, the U.K. Ministry of Defence or the U.K.
Government. The U.S. and U.K. Governments are authorized to reproduce and
distribute reprints for Government purposes notwithstanding any copyright
notation hereon.

Conrad D. Hougen1, Lance M. Kaplan2, Magdalena Ivanovska3, Federico Cerutti4,
Kumar Vijay Mishra2 and Alfred O. Hero III1
1University of Michigan, Ann Arbor, MI, USA
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>2 DEVCOM Army Research Laboratory,
Adelphi, MD, USA
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>3 BI Norwegian Business
School, Oslo, Norway
Email<EMAIL_ADDRESS>4University of Brescia, Brescia, Italy
Email<EMAIL_ADDRESS>
###### Abstract

In second-order uncertain Bayesian networks, the conditional probabilities are
only known within distributions, i.e., probabilities over probabilities. The
delta-method has been applied to extend exact first-order inference methods to
propagate both means and variances through sum-product networks derived from
Bayesian networks, thereby characterizing epistemic uncertainty, or the
uncertainty in the model itself. Alternatively, second-order belief
propagation has been demonstrated for polytrees but not for general directed
acyclic graph structures. In this work, we extend Loopy Belief Propagation to
the setting of second-order Bayesian networks, giving rise to Second-Order
Loopy Belief Propagation (SOLBP). For second-order Bayesian networks, SOLBP
generates inferences consistent with those generated by sum-product networks,
while being more computationally efficient and scalable.

##  I Introduction

Probabilistic graphical models (PGMs), such as Bayesian networks (BNs) and
Markov random fields, encode conditional dependencies between variables of
interest. In particular, BNs [1] (Section II) encode dependencies between
variables using directed acyclic graphs (DAGs), parameterized by conditional
probabilities. BNs are often used to infer the probability of values of
unobserved variables, given values from a set of evidence variables. Standard
BNs only consider first-order uncertainty, where the network conditional
probabilities are assumed to be known precisely. In these BNs, first-order
inference methods are employed to predict the conditional probabilities of the
network as point values.

Conversely, uncertain BNs (Section III) are BNs where the conditional
probabilities are not known, or are known only up to intervals or
distributions. The type of uncertainty in BNs can take several forms, leading
to different categories of uncertain BNs, such as credal networks [2],
valuation-based systems [3], subjective BNs [4], and second-order BNs [5]. A
comprehensive survey of learning and reasoning for uncertain BNs is provided
in [6].

Second-order BNs [5], in particular, consider second-order probabilities where
knowledge of the conditional probabilities is encoded by posterior
distributions, which here are distributions over probabilities. The goal of
second-order inference is to capture epistemic uncertainty in the model. Not
only do we wish to make accurate inferences: we would also like to quantify
our confidence level in the inferred values. Exact second-order inference for
second-order BNs is intractable. Recent works [7, 8] have promoted methods for
second-order inference which involve translating them into sum-product
networks (SPNs) (Section II-B) and then using the delta-method (Section III-A)
to approximate the variance of the queried probabilities by performing a
first-order Taylor expansion [9, 10]. The propagation of the means and
variances by extending exact first-order inference methods via the delta-
method leads to good approximations of the distributions of the queried
probabilities [5, 7]. However, sum-product networks do not scale well as the
size of the BN grows. As an alternative to the delta-method approximation,
Monte Carlo sampling of first-order SPN inference could more accurately
estimate the queried distribution but at even more computational expense.

One of the most efficient inference mechanisms for operating on arbitrary BNs
is a generalization of belief propagation (Section II-A) to networks
containing loops, the so-called loopy belief propagation (LBP) [11]. Here, a
loop is any multi-path between a given pair of nodes. LBP usually, but not
necessarily, converges to the correct inference [12].

Our primary contribution, Section IV, is an extension of LBP to the setting of
second-order BNs with discrete-valued variables, where the conditional
probabilities are Dirichlet-distributed. We call the resulting method Second-
Order Loopy Belief Propagation (SOLBP). Indeed, using SPNs with delta-method
approximation will constrain the size of the BN that we wish to perform
inference on. Techniques that are commonly used for second-order inference,
such as belief propagation, have so far been limited to polytree networks,
where the graph structures are disallowed from containing loops. The
experimental results in Section V show that the inferences generated by SOLBP
match those generated by the delta-method in SPNs when LPB correctly
converges, confirming correctness of SOLBP, at a lesser computational cost.
Finally, we draw our conclusions in Section VI.

##  II Background in Bayesian Networks

Associated with any BN is a directed acyclic graph (DAG), where each node is a
random variable, and links represent conditional dependencies [13]. Each child
variable, $Y$, with parents $\mathbf{X}=[X_{1},...,X_{N_{p}}]$, has an
associated conditional probability table
$p(Y=y|X_{1}=x_{1},...,X_{N_{p}}=x_{N_{p}})$. Here, we have assumed that the
BN is composed of discrete-valued variables, which is the common case. A
variable, $Y$, with no ancestors, is then associated with an unconditional
probability table $p(Y)$. Each variable $Y$ has a finite set of mutually
exclusive states, $\mathbb{Y}$.

A given graph with fixed structure—i.e., the set of nodes and edges—represents
a family of BNs parameterized by the probabilities $\bm{\theta}$, which we
interchangeably refer to as the parameters or conditionals. Each parameter
corresponds to a particular entry in the conditional or unconditional
probability tables of the BN.

The joint probability of variables $Y_{1},...,Y_{N}$ is

$p(y_{1},...,y_{N};\bm{\theta})=\prod_{i=1}^{N}\theta_{y_{i}|\mathbf{x}_{i}},$
(1)

where $y_{i}\in\mathbb{Y}_{i}$ is the value that $Y_{i}$ is assigned, and
$\mathbf{x}_{i}\in\prod_{j=1}^{N_{p,i}}\mathbb{X}_{i,j}$ represents a specific
assignment of the parents of variable $Y$. For any assignment of the parents
of $Y$, the conditionals, $\theta_{y_{i}|\mathbf{x}_{i}}$ form a distribution,
i.e.

$\sum_{y_{i}\in\mathbb{Y}_{i}}\theta_{y_{i}|\mathbf{x}_{i}}=1,\forall\mathbf{x}_{i}\in\prod_{j}\mathbb{X}_{i,j},\forall
i.$ (2)

One of the main inference tasks in BNs is to determine the conditional
probabilities of unobserved variables, given observations of evidence
variables. We will focus here on the sub-task of computing the marginals of
the nodes, i.e. determining the probability distribution of each variable in
the network $Y_{i}$ upon observing the values of a set of variables
$\mathbf{E}$, i.e. computing $p(y_{i}|\textbf{e})$, for each
$y_{i}\in\mathbb{Y}_{i}$. We denote this distribution by
$p(Y_{i}|\mathbf{e})$. If $\mathbf{E}=\emptyset$, then
$p(y_{i}|\textbf{e})=p(y_{i})$. In general the computation of a queried
probability for node $y_{i}$ requires calculations of the marginal

$p(y_{i},\mathbf{e})=\sum_{\mathbf{y}\sim\mathbf{e}\cup\\{y_{i}\\}}\left(\prod^{N}_{j=1}\theta_{y_{j}|\mathbf{x}_{j}}\right)\;,$
(3)

where “$\sim$” is used to denote compatibility of variable assignments. By the
rule of conditional probability,

$p(y_{i}|\mathbf{e})=\frac{p(y_{i},\mathbf{e})}{\sum_{y_{i}\in\mathbb{Y}_{i}}p(y_{i},\mathbf{e})}.$
(4)

Note that the computation of the probability is a nonlinear function of the
set of conditional probabilities $\Theta$ for the BN. The following
subsections review methods to compute (4).

###  II-A Belief Propagation

Belief propagation, introduced in [1], uses message passing to efficiently
compute exact inferences over polytrees. Even when the BN is not a tree, it
may be possible to segment the BN into multiple tree structures, giving rise
to the junction-tree algorithm [14], which reduces to belief propagation for
networks without loops. In belief propagation, each node in the BN both passes
and receives messages to or from neighboring nodes. Specifically, each node
receives $\lambda$-messages from its children and $\pi$-messages from its
parents, which are then used to calculate outgoing $\lambda$ and $\pi$
messages.

The $\pi$ message from a parent node $X_{1}$ to a child node $Y$ encodes the
probability of the parent variable, conditioned on the evidence coming from
$X_{1}$. If the evidence includes a particular value for $Y$, then the
internal $\bm{\pi}_{Y}$ value is a one-hot vector, i.e.,
$\pi_{Y}(y)=\delta_{y,e_{Y}}$ where $e_{Y}$ is the value of $Y$ in the
evidence $\mathbf{e}$. Otherwise, node $Y$ collects its $\pi$-messages from
its parents $\mathbf{X}=[X_{1},\ldots,X_{N_{p}}]$. Once node $Y$ receives
these messages, the node $Y$ can compute its internal $\pi$-values:

$\pi_{Y}(y)=\sum_{\bm{X}}\theta_{y|x_{1},\ldots,x_{N_{p}}}\prod_{i=1}^{N_{p}}\pi_{Y,X_{i}}(x_{i}).$
(5)

Similarly, the internal $\bm{\lambda}_{Y}$ value for node $Y$ is a one-hot
vector if the evidence includes the value of $Y$, i.e.,
$\lambda_{Y}(y)=\delta_{y,e_{y}}$. Otherwise once the $\lambda$-messages are
received by node $Y$ from all of its children,
$\mathbf{Z}=[Z_{1},\ldots,Z_{N_{c}}]$, the internal value is accumulated as a
product:

$\lambda_{Y}(y)=\prod_{j=1}^{N_{c}}\lambda_{Y,Z_{j}}(y).$ (6)

Once a node can compute its internal $\bm{\pi}_{Y}$ and $\bm{\lambda}_{Y}$
values, it can compute its value probabilities as

$p(y|\mathbf{e})=\dfrac{\lambda_{Y}(y)\pi_{Y}(y)}{\sum_{y^{\prime}\in\mathbb{Y}}\lambda_{Y}(y^{\prime})\pi_{Y}(y^{\prime})}.$
(7)

A node can only send out a message to a neighbor once it has received its
messages from all other neighbors. The node sends out a $\pi$-message to a
child node $Z_{j}$ once it has computed its internal $\bm{\pi}_{Y}$ value and
received $\lambda$-messages from all the other children in $\mathbf{Z}$. This
message is computed as

$\pi_{Z_{i},Y}(y)=\pi_{Y}(y)\prod_{j\neq i}\lambda_{Y,Z_{j}}(y),$ (8)

where $\bm{\lambda}_{Y,Z_{j}}$ is the $\lambda$-message that child $Z_{j}$
sent to node $Y$. Likewise, node $Y$ can send out a $\lambda$ message to a
parent node $X_{i}$ once it has computed its internal $\bm{\lambda}_{Y}$ value
and received $\pi$-messages from all the other parents in $\mathbf{X}$. The
message is computed as

$\lambda_{X_{i},Y}(x_{i})=\sum_{y\in\mathbb{Y}}\lambda_{Y}(y)\sum_{\bm{X}\sim
x_{i}}\theta_{y|x_{1},\ldots,x_{N_{p}}}\prod_{j\neq i}\pi_{Y,X_{i}}(x_{j}),$
(9)

where $\bm{\pi}_{Y,X_{i}}$ is the $\pi$-message that parent $X_{i}$ sent to
node $Y$.

Initially, only nodes that do not have parents can compute their internal
$\bm{\pi}_{Y}$ values where $(\ref{fprop})$ simplifies to
$\pi_{Y}(y)=\theta_{y}$. Likewise, nodes that have no children can assign
their internal $\bm{\lambda}_{Y}$ values to the non-informative identity value
for product fusion, i.e., $\lambda_{Y}(y)=1$.

Nodes which do not have parents do not send $\lambda$-messages, and similarly,
nodes without children do not send $\pi$-messages. Once all of the nodes in
the network have sent their respective $\lambda$ and $\pi$ messages, each node
$Y$ will have an associated belief, $p(y|\mathbf{e})$, given the instantiated
variables. The instantiated variables are also called evidence, which we have
denoted by the boldface $\mathbf{e}$, where we use the boldface notation to
indicate that the evidence is a vector of multiple variables.

Belief propagation is designed to provide exact inference for polytree BNs. It
can be shown that for polytree networks, all nodes will be able to send their
$\pi$\- and $\lambda-$message to each other, and the probabilities computed by
(7) are equivalent to (4). An extension – Loopy Belief Propagation [11] – is
based on the idea to run Belief Propagation on a graph containing loops using
a fixed-point iteration procedure, despite the fact that the presence of loops
does not guarantee convergence. As our contribution directly expands on LBP,
we will illustrate its details in Section IV in the context of our proposal.

###  II-B Sum-Product Networks for Inference

Another approach to inference in BNs takes advantage of a known technique for
translating BNs into Sum-Product Networks (SPNs) as pioneered by Darwiche
[15]. An SPN consists of a rooted DAG of interior operator nodes (e.g. sum or
product operations) and leaf nodes which are associated with indicator
variables. In order to convert a BN into an SPN, the basic idea is to
introduce indicator variables $\lambda_{y_{i}};\ i=1,\ldots,n$ such that

$\lambda_{y_{i}}=\begin{cases}1,&\text{if}\ Y_{i}=y_{i}\\\
0,&\text{otherwise}\end{cases}.$ (10)

When combined with the BN parameter variables $\bm{\theta}$, we may write down
a canonical polynomial for any particular BN. Given the canonical polynomial,
it is relatively straightforward to form an associated SPN graphical model, by
means of variable elimination. The specific SPN structure is not unique and
will depend on the order in which variables are eliminated.

The indicator variables of the leaf nodes in an SPN represent the assignment
states of variables from the original BN from which the SPN was derived. By
setting the indicator nodes of the SPN to either 1 or 0, we can efficiently
compute $p(\mathbf{e})$, where $\mathbf{e}$ represents the evidence, i.e. the
values of the observed variables. In order to accomplish this, we set the
values of the leaf indicators to 0 if the conditional probabilities associated
with each leaf are inconsistent with the evidence. A forward pass from the
leaf nodes to the root of the SPN performs marginalization of (1) over all
variables which are unobserved:

$p(\mathbf{e})=\sum_{\mathbf{y}\sim\mathbf{e}}\prod_{i=1}^{N}\theta_{y_{i}|\mathbf{x}_{i}}.$
(11)

An SPN readily performs this operation during a single forward pass.

Additionally, a backward pass through an SPN can be used to efficiently
compute partial derivatives of the evidence likelihood with respect to each of
the BN parameters, i.e. partial derivatives of the form

$\frac{\partial
p(\mathbf{e};\bm{\theta})}{\partial\theta_{y_{k}|{\mathbf{x}_{k}}}}.$

These partial derivatives are beneficial for computing joint probabilities
$p(y_{k},\mathbf{e})$, for each of the variables $y_{k}$. These values are the
outputs of the leaf nodes of the SPN:

$p(y_{k},\mathbf{e})=\sum_{\mathbf{x}_{k}\in\prod_{j}\mathbb{X}_{k,j}}\frac{\partial
p(\mathbf{e};\bm{\theta})}{\partial\theta_{y_{k}|{\mathbf{x}_{k}}}}\theta_{y_{k}|\mathbf{x}_{k}}.$
(12)

##  III Background in Uncertain Bayesian Networks

The conditional probabilities for BNs are either subjective in that they are
provided by domain experts or they are learned from historical data. Here,
historical data represents instantiations of the variable values sampled via
the distribution described by the BN. As a result of learning with finite
datasets, the conditional probabilities are not known precisely. Nevertheless,
the conditional probabilities are estimated from the data, usually via maximum
likelihood estimation, and are treated as exact values during inference.

Second-order inference methods instead consider that knowledge of the
conditional probabilities exhibit epistemic uncertainty encoded as the
posterior distributions of the conditionals, given training data. As described
in [13], when the dataset is complete in the sense that the values of all the
variables are visible for each instantiation, the posterior distributions for
the conditional probabilities factor as a product of Dirichlet distributions
such that the posterior for $\bm{\theta}_{Y_{i}|\mathbf{x}_{i}}$ is

$f(\bm{\theta}_{Y_{i}|\mathbf{x}_{i}})=\frac{1}{B(\bm{\alpha}_{\theta_{Y_{i}|\mathbf{x}_{i}}})}\prod_{k=1}^{|\mathbb{Y}_{i}|}\left(\theta_{y_{i,k}|\mathbf{x}_{i}}\right)^{\alpha_{\theta_{Y_{i}|\mathbf{x}_{i},k}}-1},$
(13)

where $k$ indexes the possible values of $Y_{i}$, $B(\cdot)$ is the
multidimensional beta function

$B(\bm{\alpha})=\frac{\prod_{k}\Gamma(\alpha_{k})}{\Gamma(\sum_{k}\alpha_{k})},$

and $\Gamma(\cdot)$ is the gamma function. The Dirichlet parameter
$\alpha_{k}$ is one plus the number of times in the training data that node
$Y_{i}$ took on a value of $y_{i,k}$ when the values of the nodes
$\mathbf{X}_{i}$ took on values given by $\mathbf{x}_{i}$. Remarkably,
$\bm{\theta}_{Y_{i}|\mathbf{x}_{i}}$ is independent of all other conditionals
$\bm{\theta}_{Y_{i}|\mathbf{x}^{\prime}_{i}}$ and
$\bm{\theta}_{Y_{j}|\mathbf{x}_{i}}$ for
$\mathbf{x}_{i}\neq\mathbf{x}_{j}^{\prime}$ and $i\neq j$.

In light of the posterior distributions for the conditional probabilities,
second-order BN inference computes the distribution of the queried
probabilities in (4). Exact computation is intractable, but it can be
approximated by leveraging the delta-method (see next subsection) that
approximates the mean and covariance matrix for the queried values. To that
end, we approximate the posteriors as Gaussian. In light of the Dirichlet
distribution given by (13), the means for the values of
$\bm{\theta}_{Y_{i}|\mathbf{x}_{i}}$ are given by

$\mu[\theta_{y_{i,k}|\mathbf{x}_{i}}]=\frac{\alpha_{\theta_{Y_{i}|\mathbf{x}_{i},k}}}{S_{Y_{i}|\mathbf{x}_{i}}}$
(14)

where the denominator

$S_{Y_{i}|\mathbf{x}_{i}}=\sum_{k=1}^{|\mathbb{Y}_{i}|}\alpha_{\theta_{Y_{i}|\mathbf{x}_{i},k}}$
(15)

is known as the Dirichlet strength. The covariance for the terms in
$\bm{\theta}_{Y_{i}|\mathbf{x}_{i}}$ is

$\sigma[\theta_{y_{i}|\mathbf{x}_{i}},\theta_{y^{\prime}_{i}|\mathbf{x}_{i}}]=\frac{\mu[\theta_{y_{i}|\mathbf{x}_{i}}](\delta_{y_{i},y^{\prime}_{k}}-\mu[\theta_{y^{\prime}_{i}|\mathbf{x}_{i}}])}{S_{Y_{i}|\mathbf{x}_{i}}+1}\;.$
(16)

###  III-A The Delta-Method

In general, reasoning in a BN can be viewed as passing the conditional
probabilities through a nonlinear operation, e.g. (4). The delta-method
introduced in [9] approximates the mean and variance of the output of the
nonlinear operation via a first-order Taylor series.

The MeanVAR introduced in [10] uses the delta-method to compute (4) directly
from the mean and variances of the conditional probabilities. Alternatively,
the second-order belief propagation (SOBP) in [5] provides a more scalable
solution by applying the delta-method to each step to compute the mean and
covariances of the $\pi-$ and $\lambda-$messages (see (5)-(9)). Finally, [7,
8] apply the delta-method for each sum and product operation in the
probabilistic circuit.

In general, consider a set of $N$ input messages
$\bm{\chi}_{i}=[\chi_{i}(x):x\in\mathbb{X}_{i}]$ for $i=1,\ldots,N$ with means
$\mu[\bm{\chi}_{i}]$ and covariance $\sigma[\bm{\chi}_{i}]$ that are
statistically independent of each other. Let the output message be related to
the inputs via

$\bm{\upsilon}=\mathbf{f}(\bm{\chi}_{1},\ldots,\bm{\chi}_{N}),$ (17)

where $\mathbf{f}=[f_{\upsilon(u)}(\cdot):u\in\mathbb{U}].$ The delta-method
linearizes the function $\mathbf{f}$ about the expected values of the input
messages to approximate the output message as

$\bm{\upsilon}\approx\mathbf{f}(\mu[\bm{\chi}_{1}],\ldots,\mu[\bm{\chi}_{N}])+\sum_{i=1}^{N}\mathbf{J}_{\chi_{i}}(\bm{\chi}_{i}-\mu[\bm{\chi}_{i}]),$
(18)

where $\mathbf{J}_{\chi_{i}}$ is the Jacobian of the function $\mathbf{f}$
with elements representing the partial derivatives $\frac{\partial
f_{\upsilon(u)}}{\partial\chi(x)}$ evaluated at the mean values of the
$\bm{\chi}_{i}$ messages.

Because of the statistical independence between the messages, it follows that
the mean and covariance matrices for the output message are

$\mu[\bm{\upsilon}]=\mathbf{f}(\mu[\bm{\chi}_{1}],\ldots,\mu[\bm{\chi}_{N}]),~{}\sigma[\bm{\upsilon}]=\sum_{i=1}^{N}\mathbf{J}_{\chi_{i}}\sigma[\bm{\chi}_{i}]\mathbf{J}_{\chi_{i}}^{T}.$
(19)

The delta-method can extend an exact first-order BN inference engine into a
second-order inference method that computes the means and covariances of the
queried variables. In light of (19), the mean values are treated as the point
probabilities in first-order inference.

Any exact inference method can be extended to a second-order method via the
delta-method. For instance, second-order belief propagation was first
developed in [5] for binary-valued polytree BNs. Similarly, [7] provides the
second-order extension of sum-product networks for general discrete-valued
BNs. Note that any exact first-order BN inference method is a composition of
operations to transform the conditional probabilities into the inference. As
argued in [8], the variances computed by the second-order extensions of these
inference methods through the delta-method are equivalent due to the chain
rule of calculus.

##  IV SOLBP Method

The main contribution of this paper is the extension of loopy belief
propagation (LBP) for general BNs via the delta-method, which we refer to as
SOLBP. This is motivated by the fact that belief propagation is more scalable
than inference over sum-product networks. When LBP does converge to the exact
inference, it represents a composition of operations leading to the exact
inference function. Then, by chain-rule arguments in [8], SOLBP should lead to
the same variances as other exact second-order extensions.

The basic concept is that all parent-to-child and child-to-parent messages are
indexed. First, the messages are initialized to neutral values that are
uninformative. A message is selected at random without replacement. The mean
and covariance of that message is updated based on the current values of its
dependent messages. Once the selection of messages have been exhausted, the
round is completed. If the change in the mean values of the messages is small,
the computation terminates. Otherwise, another round is initiated. During a
round, once all the parent or children input message have been updated, the
node updates its internal $\pi-$ or $\lambda-$ messages, respectively.
Algorithm 1 summarizes the SOLBP steps.

Algorithm 1 Second-order loopy belief propagation.

1: Input: Evidence $\mathbf{e}$, $\mu[\theta_{Y|\mathbf{x}}]$ and
$\sigma[\theta_{Y|\mathbf{x}}]$ for all $\theta_{Y|\mathbf{x}}\in\Theta$

2: Output: $\mu[\mathbf{p}(Y|\mathbf{e})]$ and
$\sigma[\mathbf{p}(Y|\mathbf{e})]$ for all nodes $Y$ in the BN.

3:

4: Set $\mathcal{M}^{*}$ to be the set of all possible inter-node messages

5: Set $\mathcal{M}\leftarrow\mathcal{M}^{*}$

6: For all $\bm{m}\in\mathcal{M}$, set
$\mu[\mathbf{m}]=\left\\{\begin{array}[]{cc}\frac{1}{|\mathbb{M}|}\mathbf{1}&\mbox{for
$\pi$-message},\\\ \mathbf{1}&\mbox{for
$\lambda$-message},\end{array}\right.\sigma[\mathbf{m}]=\mathbf{0},\mu^{*}[\mathbf{m}]=\mathbf{0},$

7: for each node $Y$ represented in $\mathbf{e}$ set
$\mu[\pi_{Y}(y)]=\mu[\lambda_{Y}(y)]=\delta_{y,e_{y}}$ and
$\sigma[\pi_{Y}(y)]=\sigma[\lambda_{Y}(y)]=0$ for all $y\in\mathbb{Y}$

8: while
$\max\\{|\mu[\bm{m}]-\mu^{*}[\bm{m}]|:\mathbf{m}\in\mathcal{M}^{*}\\}>\epsilon$
do

9: $\mu^{*}[\bm{m}]\leftarrow\mu[\bm{m}]$ for all $m\in\mathcal{M}^{*}$

10: while $\mathcal{M}\neq\emptyset$ do

11: Randomly select $\bm{m}\in\mathcal{M}$,
$\mathcal{M}\leftarrow\mathcal{M}\setminus\\{\bm{m}\\}$

12: if $\bm{m}$ is a $\pi$-message then

13: Update $\mu[\bm{m}]$ and $\sigma[\bm{m}]$ via (8) and (24), respectively

14: if all $\pi$-messages to node $Z$ have been updated in current round and
$Z$ is not in $\mathbf{e}$ then

15: update $\mu[\bm{\pi}_{Z}]$ and $\sigma[\bm{\pi}_{Z}]$ via (5) and (26),
respectively

16: end if

17: else

18: Update $\mu[\bm{m}]$ and $\sigma[\bm{m}]$ via (9) and (20), respectively

19: if all $\lambda$-messages to node $X$ have been updated in current round
and $X$ is not in $\mathbf{e}$ then

20: update $\mu[\bm{\lambda}_{X}]$ and $\sigma[\bm{\lambda}_{X}]$ via (6) and
(29), respectively

21: end if

22: end if

23: end while

24: $\mathcal{M}\leftarrow\mathcal{M}^{*}$

25: end while

26: Compute $\mu[\mathbf{p}(Y|\mathbf{e})]$ and
$\sigma[\mathbf{p}(Y|\mathbf{e})]$ via (7) and (31), respectively, for all
nodes in the BN

The initial $\pi$-messages that each node $Y$ sends to each child $Z$ are set
so that the mean and covariance matrix of an individual message is

$\mu[\bm{\pi}_{Z,Y}]=\frac{1}{|\mathbb{Y}|}\mathbf{1}_{\mathbb{|Y|}},\hskip
3.61371pt\mbox{and}\hskip 3.61371pt\sigma[\bm{\pi}_{Z,Y}]=\mathbf{0},$

respectively. Note that the mean values of the message must sum to one. Also,
$\mathbf{1}_{N}$ is an $N$-length vector of all ones. Similarly, the initial
$\lambda$-messages that each node sends to each parent $X$ are set so that the
mean and covariance matrix of a single message is

$\mu[\bm{\lambda}_{X,Y}]=\mathbf{1}_{\mathbb{|X|}},\hskip
3.61371pt\mbox{and}\hskip 3.61371pt\sigma[\bm{\lambda}_{X,Y}]=\mathbf{0},$

respectively. Unlike the $\pi$-messages, the $\lambda$-messages do not need to
be normalized to have their mean values to sum to one. These message values do
not have any influence on biasing beliefs to any particular variable value in
the execution on the belief propagation operations. The zero covariances also
do not influence the variances of the belief. Initially, all the influence of
beliefs lies in the mean and covariances of the conditional probabilities and
the values of the evidence $\mathbf{e}$.

The update of the message to a parent from a node is accomplished by the back
propagation operation in (9). The means are propagated in the obvious way
(i.e., the input messages in (9) are now the respective means). The covariance
for message from node $Y$ to parent $X_{i}$ is

$\displaystyle\sigma[$
$\displaystyle\lambda_{X_{i},Y}(x_{i})]=\sum_{\mathbf{X}\sim
x_{i}}\mathbf{J}_{\theta_{Y|\mathbf{x}}}\sigma[\bm{\theta}_{Y|\mathbf{x}}]\mathbf{J}_{\theta_{Y|\mathbf{x}}}^{T}$
(20) $\displaystyle+\sum_{j\neq
i}\mathbf{J}_{\pi_{Y,X_{j}}}\sigma[\bm{\pi}_{Y,X_{i}}]\mathbf{J}_{\pi_{Y,X_{j}}}^{T}+\mathbf{J}_{\lambda_{Y}}\sigma[\bm{\lambda}_{Y}]\mathbf{J}_{\lambda_{Y}}^{T}\;,$

where the elements of the Jacobians $\mathbf{J}_{\theta_{Y|\mathbf{x}}}$,
$\mathbf{J}_{\pi_{Y,X_{j}}}$, and $\mathbf{J}_{\lambda_{Y}}$ are given by:

$\frac{\partial\lambda_{X_{i},Y}(x_{i})}{\partial\theta_{y|\mathbf{x}}}=\mu[\lambda_{Y}(y)]\prod_{j\neq
i}\mu[\pi_{Y,X_{j}}(x_{j})]\;,$ (21)
$\displaystyle\frac{\partial\lambda_{X_{i},Y}(x_{i})}{\partial\pi_{Y,X_{j}}(x_{j})}=$
(22) $\displaystyle\sum_{Y}\mu[\lambda_{Y}(y)]\left(\sum_{\mathbf{x}\sim
x_{i},x_{j}}\mu[\theta_{y|\mathbf{x}}]\prod_{\ell\neq
i,j}\mu[\pi_{Y,X_{\ell}}(x_{\ell})]\right)\;,$

and

$\frac{\partial\lambda_{Y,X_{i}}(x_{i})}{\partial\lambda_{Y}(y)}=\sum_{\mathbf{X}\sim
x_{i}}\mu[\theta_{y|\mathbf{x}}]\prod_{j\neq i}\mu[\pi_{Y,X_{j}}(x_{j})]\;,$
(23)

respectively.

Similarly, the computation for the mean of the message from node $Y$ to its
child $Z_{j}$ is computed by (8) in an obvious manner. The covariance is
computed as

$\sigma[\bm{\pi}_{Z_{i},Y}]=\mathbf{J}_{\pi_{Y}}\sigma[\bm{\pi}_{Y}]\mathbf{J}_{\pi_{Y}}+\sum_{j\neq
i}\mathbf{J}_{\lambda_{Y,Z_{j}}}\sigma[\bm{\lambda}_{Y,Z_{j}}]\mathbf{J}_{\lambda_{Y,Z_{j}}}^{T},$
(24)

where again the Jacobians $\mathbf{J}_{\pi_{Y}}$ and
$\mathbf{J}_{\lambda_{Y,Z_{j}}}$ are diagonal with non-zero values

$\frac{\partial\pi_{Z_{i},Y}(y)}{\partial\pi_{Y}(y)}=\frac{\mu[\pi_{Z_{i},Y}(y)]}{\mu[\pi_{Y}(y)]}\hskip
7.22743pt\mbox{and}\hskip
7.22743pt\frac{\partial\pi_{Z_{i},Y}(y)}{\partial\lambda_{Y,Z_{j}}(y)}=\frac{\mu[\pi_{Z_{i},Y}(y)]}{\mu[\lambda_{Y,Z_{j}}(y)]},$
(25)

respectively.

The update of an internal $\pi$ value is accomplished by the forward
propagation in (5). Again, the means are propagated in the obvious way, and
the covariance is updated as

$\displaystyle\sigma[\bm{\pi}_{Y}]=$
$\displaystyle\sum_{\bm{X}}\mathbf{J}_{\theta_{Y|\mathbf{x}}}\sigma[\bm{\theta}_{Y|x_{1},\ldots,x_{N_{p}}}]\mathbf{J}_{\theta_{Y|\mathbf{x}}}^{T}+$
(26)
$\displaystyle\sum_{i=1}^{N_{p}}\mathbf{J}_{\pi_{Y,X_{i}}}\sigma[\bm{\pi}_{X_{i}}]\mathbf{J}_{\pi_{Y,X_{i}}}^{T}\;,$

where the Jacobian

$\mathbf{J}_{\theta_{Y|\mathbf{x}}}=\left(\prod_{i=1}^{N_{p}}\mu[\pi_{Y,X_{i}}(x_{i})]\right)\mathbf{I},$
(27)

$\mathbf{I}$ is the identity matrix, and the elements of the Jacobian
$\mathbf{J}_{\pi_{Y,X_{i}}}$ are

$\frac{\partial\pi_{Y}(y)}{\partial\pi_{Y,X_{i}}(x_{i})}=\sum_{\mathbf{X}\sim
x_{i}}\mu[\theta_{y|\mathbf{x}}]\prod_{j\neq i}\mu[\pi_{Y,X_{j}}(x_{j})].$
(28)

Similarly, the internal $\lambda-$value is updated via (6) for the mean, and
the covariance update is given by

$\sigma[\bm{\lambda_{Y}}]=\sum_{j=1}^{N_{c}}\mathbf{J}_{\lambda_{Y,Z_{j}}}\sigma[\bm{\lambda}_{Y,Z_{j}}]\mathbf{J}_{\lambda_{Y,Z_{j}}}^{T}\;,$
(29)

where the Jacobians $\mathbf{J}_{\lambda_{V}(X_{i})}$ are diagonal with
diagonal terms

$\frac{\partial\lambda_{Y}(y)}{\partial\lambda_{Y,Z_{j}}(y)}=\frac{\mu[\lambda_{Y}(y)]}{\mu[\lambda_{Y,Z_{i}}(y)]}.$
(30)

Once all second-order statistics between the $\pi-$ and $\lambda$ messages
between the nodes have been updated along with the internal values, the round
is completed. Then, the mean of the queried probabilities for each node is
computed via (7), and the covariances are given by

$\sigma[p(Y|\mathbf{e})]=\mathbf{J}_{\lambda_{Y}}\sigma[\bm{\lambda}_{Y}]\mathbf{J}_{\lambda_{Y}}+\mathbf{J}_{\pi_{Y}}\sigma[\bm{\pi}_{Y}]\mathbf{J}_{\pi_{Y}}^{T}$
(31)

where the elements of the Jacobians $\mathbf{J}_{\lambda_{Y}}$ and
$\mathbf{J}_{\pi_{Y}}$ are given by

$\frac{\partial{p(y|\mathbf{e})}}{\partial\lambda_{Y}(y^{\prime})}=\frac{\mu[p(y|\mathbf{e})](\delta_{y,y^{\prime}}-\mu[p(y^{\prime}|\mathbf{e})])}{\mu[\lambda_{Y}(y^{\prime})]}\hskip
14.45377pt\mbox{and}$ (32)
$\frac{\partial{p(y|\mathbf{e})}}{\partial\pi_{Y}(y)}=\frac{\mu[p(y|\mathbf{e})](\delta_{y,y^{\prime}}-\mu[p(y^{\prime}|\mathbf{e})])}{\mu[\pi_{Y}(y^{\prime})]},$
(33)

respectively.

##  V Numerical Evaluation

This section provides empirical evidence that SOLBP determines the mean and
variance of queries in the form of $p(x_{i}|\mathbf{e})$ for each unobserved
variable that are consistent with those computed via second-order SPN (SOSPN)
as described in [7, 8]. Furthermore, empirical results demonstrate the
computational efficiency of SOLBP relative to SOSPN, thus showing the
advantages of SOLBP for larger DAG BNs.

In order to make the appropriate comparisons, we set up a test framework in
MATLAB which takes BNs as input. During each experimental trial, conditional
probabilities, $\bm{\theta}$, of the ground truth BN were individually sampled
from a uniform Dirichlet distribution. Then, a sparse but complete set of
observations of the network variables were sampled from the appropriate
discrete distributions associated with the ground truth BN. These
$N_{\textrm{train}}$ examples form the training data for an exact Bayesian
learning process. During each trial, the first step utilizes the
$N_{\textrm{train}}$ number of complete observations of the network to compute
the parameters of the Dirichlet distribution by counting the number of
examples consistent with variable values (see Section III).

For the inference step, we instantiate a specific set of BN variables as the
observed evidence, $\mathbf{e}$, leaving the remaining variables unobserved.
Then, we infer $p(Y_{i}|\mathbf{e})$ via SPN over the BN with ground truth
conditional probabilities to establish the ground truth inferences. Finally,
we infer the mean and variances of the query via SOLBP or SOSPN over the
uncertain BN learned via the $N_{\textrm{train}}$ samples. These inferred
means and variances are recorded over $N_{\textrm{runs}}$ number of trials,
where we used $N_{\textrm{runs}}=1000$, and a new ground truth BN was
generated after every $100$ trials. Note that, while the ground truth BN
remained static for $100$ trials at a time, a new uncertain BN is learned on
each trial since new Monte Carlo observations of the variables are generated
on a per-trial basis.

In order to evaluate accuracy, the recorded means and variances from both the
SOSPN and SOLBP methods are plotted against each other. Ideally, we would like
to see that the inferences match in all cases, resulting in a diagonal line of
slope one. Additionally, we follow the method from [5] for generating plots of
the desired confidence bound divergence (DeCBoD) for each second-order
inference method by establishing confidence bounds from the means and
variances. The DeCBoD determines the fractional rate at which the ground truth
probabilities fall within confidence intervals of significance level $\gamma$
around the inferred probabilities. For $\gamma\in[0,0.99]$, we compute the
fractional in-bounds rate and plot it against $\gamma$.

We tested SOLBP on three sets of BNs, varying types and sizes: (1) three-node
networks; (2) small loopy networks; and (3) a large loopy network. Recall, one
of the primary motivations for SOLBP is that, similar to the SOSPN method from
[7], SOLBP can be applied to DAGs with loops, in addition to polytrees.
However, we first want to verify that SOLBP can perform equivalently to the
SOSPN method on polytrees, so we show the results for both polytrees and DAGs
with loops in the following subsections. We provide a discussion on the
relative computational efforts in Section V-D.

###  V-A Three-Node Networks

We first tested three simple three-node networks in Figure 1, with binary-
valued variables. Values of the means and variances of the inferred
probabilities from SOLBP were plotted against the means and variances inferred
by the SOSPN method. Evidence variables were chosen as the nodes incident with
exactly one edge, regardless of edge direction, which are colored grey in the
BN figures for clarity. Figure 2, compare the inferred means and variances for
the three-node chain network over the two methods. The other two networks lead
to similar plots.

Figure 1: (a) 3-node chain, (b) 3-node tent, (c) 3-node V.

(a)

(b)

Figure 2: 3-node chain inferred means (left) and variances (right).

The results demonstrate that SOLBP generates inferences which exactly match
those generated by the SPN method, up to numerical precision, for these three
polytree BNs. Given that the inferred means and variances match, we also
expect matching DeCBoD curves. Figure 3 shows a sample DeCBoD plot for the
three-node chain network from Figure 1(a), with the results from SOLBP
identified by red crosses and the SPN results identified by the blue lines.
DeCBoD plots for the remaining polytrees are similar. The DeCBOD curves show
that actual confidence matches well with the desired confidence, meaning that
the inferred distributions (mean and variance) is well calibrated to the true
uncertainty.

Figure 3: DeCBoD plot for the 3-node chain network.

###  V-B Small Loopy Networks

Next, we verified inference accuracy on two basic networks with loops, shown
in Figure 4. In the case of loopy networks, the number of iterations of SOLBP
may affect the resulting inferences, depending on the stopping condition used.
Here, we ran SOLBP until the absolute maximum difference between any inferred
means from the current and prior iteration were bounded by $\epsilon=1\times
10^{-8}$. The observed nodes in the evidence set are shaded grey in Figure 4.

Figure 5 compares the resulting inferred means and variances for SOLBP and
SOSPN over the diamond network. The triangle network produces similar plots.
Figure 6 shows the DeCBoD plots for the triangle and diamond BNs.

Figure 4: (a) Triangle loopy BN, (b) Diamond loopy BN.

(a)

(b)

Figure 5: Diamond BN inferred means (left) and variances (right).

(a)

(b)

Figure 6: DeCBoD plots for triangle loopy BN (left) and diamond loopy BN
(right)

Figure 5 show that the SOLBP and SPN methods produce inferences of means and
variances which are similar. However, unlike the polytree networks, the values
do not all match completely, meaning that some of the scatter points lie off
the diagonal. In this particular case with the chosen evidence variables,
simply using a more strict stopping condition would produce a complete match,
in a finite number of iterations, which we verified. However, for general,
multi-loop loopy networks, SOLBP may never converge to the correct values, as
is the case with first-order LBP [12]. On the other hand, for polytrees,
within a few runs of message passing, all nodes will have sent and received
$\lambda$ and $\pi$-messages with means and covariances consistent with exact
inference methods. For single-loop BNs, we would expect eventual convergence
of SOLBP over an infinite number of iterations, similar to first-order LBP.
Despite the theoretical limitations, empirical evidence demonstrates that the
stopping condition we chose for SOLBP produces approximate means and variances
that are close to the values generated by the SPN method for small, single-
loop networks.

###  V-C A Large Loopy Network

On a more practical front, we wanted to verify inference approximation
accuracy of SOLBP for a larger network with more parameters and multiple
loops. Therefore, we created a BN with 21 nodes, as shown in Figure 7.
Approximately half of the BN nodes were chosen randomly as evidence variables,
which are shaded grey in the figure. While the previous experiments were
limited to binary-valued BN variables, we allowed the BN variables to be
either binary or tri-valued, to increase the number of total parameters in the
BN. We are interested in testing more parameters, both due to potential noise
introduced by the increased degrees of freedom as well as to demonstrate
potential limitations of either SOLBP or the SOSPN method.

Figure 7: A larger 21-node BN with multiple loops. Grey nodes were used as
evidence for the experimental results shown.

Running SOLBP using the same experimental setup on the 21-node network yielded
inferred means and variances as shown in Figures 8 and 9. We see that the
SOLBP inferences continue to approximate the inferred values generated by the
SOSPN method with relatively high fidelity, though there is a noticeable
amount of noise introduced in the larger network setting. As discussed in [11]
and [12] for first-order LBP, convergence is not guaranteed in BNs with
multiple loops, and this may even depend on the values of the conditional
probabilities, in addition to the structure. This applies to SOLBP as well.
Furthermore, the selection of evidence nodes can impact the rate at which
SOLBP converges. It turns out that one can simply decrease the threshold
$\epsilon$ used as the stopping condition, or force SOLBP to run for
additional iterations, and this will improve the approximation accuracy of
SOLBP.

(a)

(b)

Figure 8: 21-node BN inferred means (left) and variances (right). Figure 9:
DeCBoD plot for the 21-node network.

###  V-D Computational Effort

Importantly, SOLBP produces inferences that are nearly identical to those
produced by the SOSPN method, and SOLBP also scales much better in terms of
computational effort. The size of a SPN grows exponentially with respect to
the size of the BN as described in [15]. Table I enumerates run-time metrics
indicating that SOLBP is able to execute the second-order inference process
multiple times faster than the SOSPN method for moderately sized BNs. This
confirms our intuition that SOLBP should be much more scalable than the SPN
method.

TABLE I: Comparison of compute time between SOSPN and SOLBP for multiple BN
structures.

##  VI Conclusion

We demonstrated second-order inference using loopy belief propagation. This
enables the scalability of second-order belief propagation while in turn
allowing for inference over all BNs. Empirical results demonstrate that SOLBP
accurately captures the epistemic uncertainty due to limited training data.
The calibration accuracy is comparable to that of SOSPN. Moreover, SOLBP is
signicantly more scalable than SOSPN.

The consistency between SOLBP and SOSPN does increase with more iterations of
SOLPB, with a tradeoff of computational efficiency. Fortunately, the DeCBoD
plots demonstrate sufficient accuracy of the epistemic uncertainty even for a
small number of iterations. Future work will systematically study the accuracy
and efficiency of SOLBP over larger and more complex BNs. Futhermore, we plan
to develop theoretical analysis to show why variances from SOLBP should be
comparable to that of SOSPN.

## References

  * [1] J. Pearl, “Fusion, propagation, and structuring in belief networks,” _Artificial Intelligence_ , vol. 29, no. 3, pp. 241–288, Sep. 1986.
  * [2] M. Zaffalon and E. Fagiuoli, “2U: An exact interval propagation algorithm for polytrees with binary variables,” _Artificial Intelligence_ , vol. 106, no. 1, pp. 77–107, 1998.
  * [3] P. P. Shenoy, “A valuation-based language for expert systems,” _Int. Journal of Approximate Reasoning_ , vol. 3, no. 2, pp. 383–411, 1989.
  * [4] M. Ivanovska, A. Jøsang, L. Kaplan, and F. Sambo, “Subjective networks: Perspectives and challenges,” in _Proc. of the 4th International Workshop on Graph Structures for Knowledge Representation and Reasoning_ , Buenos Aires, Argentina, Jul. 2015, pp. 107–124.
  * [5] L. Kaplan and M. Ivanovska, “Efficient belief propagation in second-order Bayesian networks for singly-connected graphs,” _International Journal of Approximate Reasoning_ , vol. 93, pp. 132–152, 2018.
  * [6] J. Rohmer, “Uncertainties in conditional probability tables of discrete Bayesian belief networks: A comprehensive review,” _Engineering Applications of Artificial Intelligence_ , vol. 88, p. 103384, 2020.
  * [7] F. Cerutti, L. M. Kaplan, A. Kimmig, and M. Şensoy, “Handling epistemic and aleatory uncertainties in probabilistic circuits,” _Machine Learning_ , pp. 1–43, 2022.
  * [8] L. Kaplan, M. Ivanovska, K. V. Mishra, and F. Cerutti, “Inference in second-order Bayesian networks,” in preparation.
  * [9] G. D. Kleiter, “Propagating imprecise probabilities in Bayesian networks,” _Artificial Intelligence_ , vol. 88, no. 1-2, pp. 143–161, 1996.
  * [10] T. Van Allen, A. Singh, R. Greiner, and P. Hooper, “Quantifying the uncertainty of a belief net response: Bayesian error-bars for belief net inference,” _Artificial Intelligence_ , vol. 172, no. 4-5, pp. 483–513, 2008\.
  * [11] K. P. Murphy, Y. Weiss, and M. I. Jordan, “Loopy belief propagation for approximate inference: An emprical study,” in _Proc. of the 15th Conf. on Uncertainty in Artificial Intelligence_ , 1999, pp. 467–475.
  * [12] A. T. Ihler, J. W. Fisher III, and A. S. Willsky, “Loopy belief propagation: Convergence and effects of message errors,” _Machine Learning Research_ , vol. 6, no. 31, pp. 905–936, 2005.
  * [13] F. V. Jensen and T. D. Nielsen, _Bayesian Networks and Decision Graphs_. New York: Springer, 2007.
  * [14] A. L. Madsen and F. V. Jensen, “Lazy propagation: A junction tree inference algorithm based on lazy evaluation,” _Artificial Intelligence_ , vol. 113, no. 1-2, pp. 203–245, 1999.
  * [15] A. Darwiche, “A differential approach to inference in Bayesian networks,” _Journal of the ACM_ , vol. 50, no. 3, pp. 280–305, 2003.