Average-Case Analysis of a Search Algorithm for Estimating Prior and Posterior Probabilities in Bayesian Networks with Extreme Probabilities

1 Introduction

This paper provides a general purpose search-based tech­ nique for computing posterior probabilities in arbitrarily structured discrete 1 Bayesian networks.

Implementations of Bayesian networks have been placed into three classes [Pearl, 1988; Henrion, 1990]:

1. Exact methods that exploit the structure of the network to allow efficient propagation of evidence [Pearl, 1988; Lauritzen and Spiegelhalter, 1988; Jensen et al., 1990].
2. Stochastic simulation methods that give estimates of probabilities by generating samples of instantia­ tions of the network, using for example Monte Carlo techniques (see [Henrion, 1990]).
3. Search-based approximation techniques that search through a space of possible values to estimate prob­ abilities.

At one level, the method in this paper falls into the ex­ act class; if it is allowed to run to completion, it will have computed the exact conditional probability in a Bayesian network. It, however has the extra feature that it can be stopped before completion to give an answer, with a known error. Under certain distribution assumptions We write w I= X= c (read "X has value c in world w") if variable X is assigned value c in world w. Each variable has exactly one value in a possible world; thus we can see a random variable as a set of mutually exclusive and covering propositions (i.e., the propositions of the form X= v for each v E vals(X)). A proposition is a logical formula made up of atomic propositions and the usual logical connectives. We define the truth of propositions in a world using the normal truth tables (e.g., w I= a /1. {3 if w I= a and w I= {3).

*Scholar, Canadian Institute for Advanced Research.

1 All of the variables have a discrete, and here even finite, set of possible values.

2 This could have also been presented as joint distribu­ tions, with probabilistic assignments to the possible worlds corresponding to joint distributions. If that view suits you, then please read 'possible worlds' as 'elementary events in a joint distribution' [Pearl, 1988, p. 33).

Probability is a function from propositions to [0, 1], defined by

Conditional probability is defined by:

If {3 is a conjunction of observations, then P( ai,B) is the posterior probability of a. P( a) is the prior prob­ ability of formula a.

2.2 Searching possible worlds

For a finite number of variables with a finite number of values, we can compute the probabilities directly, by enumerating the possible worlds. This is, however, com­ putationally expensive as there are exponentially many of these (the product of the sizes of the domains of the variables).

The idea behind the search method presented in this paper is motivated by considering the questions:

• Can we estimate the probabilities by only enumer­ ating a few of the possible worlds?
• How can we enumerate just a few of the most probable possible worlds?
• Can we estimate the error in our probabilities?
• For what cases does the error get small quickly?
• How fast does it converge to a small error?

This paper sets out to answer these questions, for the case where the distribution is given in terms of Bayesian networks.

2.3 Bayesian Networks

A Bayesian network [Pearl, 1988] is a graphical repre­ sentation of (in)dependence amongst random variables. A Bayesian network is a directed acyclic graph where the nodes represent random variables. If there is an arc from variable B to variable A, B is said to be a par­ ent of A. The independence assumption of a Bayesian network says that each variable is independent of its non­ descendents given its parents.

Suppose we have a Bayesian network with random variables X1 , ... ,Xn. The parents of X; are written as IIx, = (x;1 ,···,X;~,)·

Associated with the Bayesian network are conditional probability tables which gives the marginal probabilities of the values of X; depending on the values of its parents Ilx,. This consists of, for each Vj E vals(Xj), probabili­ ties of the form

For any probability distribution, we can compute a joint distribution by

This is often given as the formal definition of a Bayesian network and is shorthand for the propositions that the variables have particular values. That is

An assignment of values to all the variables corre­ sponds to a possible world. For the rest of this paper we will equate a possible world with an assignment of values to all of the variables, and use the two interchangeably (e.g., in using the complete description on the left hand side of the 'I=').

3 Enumerating possible worlds in order

3.1 Ordering the variables

The first thing we do is to impose a total ordering on the variables that is consistent with the ordering of the Bayesian network. We index the random variables X1 , ... ,Xn such that the parents of a node have a lower index than the node. This can always be done as the nodes in a Bayesian network form a partial ordering. If the parents of X; are Ilx, = ( X; 1 , · · ·, X; 1 , ). the total ordering preserves ii < i.

3.2 Search Tree

We are now in a position to determine a search tree for Bayesian networks 3 .

Definition 3.1 A partial description is a tuples of values (v1 , ···,vi} where v; is an element of the domain of variable X;.

The search tree has nodes labelled with partial de­ scriptions, and is defined as follows:

• The root is labelled with the empty tuple (}.
• The children of node labelled with (v1 , · · · , vi} are the nodes labelled with (vi.···, Vj, v} for each v E vals(Xi+d· In other words, the children of a node correspond to the possible values of the next vari­ able in the total ordering.
• The leaves of the tree are tuples of the form (vi. .. ·, Vn}·

3 This search tree is the same as the probability tree of [Howard and Matheson, 1981) and corresponds to the se­ mantic trees used in theorem proving [Chang and Lee, 1973, Section 4.4), but with random variables instead of comple­ mentary literals.

Figure 1: Basic search algorithm

Partial description (vi,···, v;) corresponding to the variable assignment

There is a one to one correspondence between leaves of the tree and possible worlds (or complete assignments to the variables).

We associate a probability with each node in the tree. The probability of the partial description (vi,···, v;) is the probability of the corresponding proposition:

This is well defined as, due to our variable ordering, all of the parents of each variable X;, 0 < i ~ j has a value in the partial description.

The following lemma can be trivially proved, and is the basis for the search algorithm.

Lemma 3.2 The probability of a node is equal to the sum of the probabilities of the leaves that are descendents of the node.

This lemma lets us bound the probabilities of possible worlds by only generating a few of the possible worlds and placing bounds on the sizes of the possible worlds we have not generated.

3.3 Searching the Search Tree

To implement the computation of probabilities, we carry out a search on the search tree, and generate some of the most likely possible worlds. There are man:y different search methods that can be used [Pearl, 1984].

Figure 1 gives a generic search algorithm that can be varied by changing which element is chosen from the queue. There is a priority queue Q of nodes, and a set W of generated worlds. We remove a node (e.g., the most likely); either it is a leaf (if j = n) in which case it is added to W, or else its children are added to Q.

Note that each partial description can only be gener­ ated once. There is no need to check for multiple paths or loops in the search. This simplifies the search, in that we do not need to keep track of a CLOSED list or check whether nodes are already on the OPEN list ( Q in Figure 1) [Pearl, 1984].

No matter which element is chosen from the queue at each time, this algorithm halts and when it halts W is the set of all tuples corresponding to possible worlds.

4 Estimating the Probabilities

If we let the above algorithm run to completion we have an exponential algorithm for enumerating the possible worlds that can be used for computing the prior proba­ bility of any proposition or conjunction of propositions. This is not, however, the point of this algorithm. The idea is that we want to stop the algorithm part way through, and determine any probability we want to com­ pute.

We use W, at the start of an iteration of the while loop, as an approximation to the set of all possible worlds. This can be done irrespective of the search strategy used.

4.1 Prior Probabilities

Suppose we want to compute P(g). At any stage (at the start of the while loop), the possible worlds can be divided into those that are in W and those that will be generated from Q.

We can easily compute the first of these sums, and can bound the second. The second sum is greater than zero and is less than the sum of the probabilities of the par­ tial descriptions on the queue (using Lemma 3.2). This means that we can bound the probabilities of a propo­ sition based on enumerating just some of the possible worlds. Let

Lemma 4.1 Pfv ~ P(g) ~ Pfv + PQ.

As the computation progresses, the probability mass in the queue PQ approaches zero and we get a better refinements on the value of P(g). Note that PQ is mono­ tonically non-increasing through the loop (i.e PQ stays the same or gets smaller through the loop). This thus forms the basis of an "anytime" algorithm for Bayesian networks.

4.2 Posterior Probabilities

The above analysis was for finding the prior probability of any proposition. If we want to compute the posterior probability of some g given some observations obs, we can use the definition of conditional probability, and use

We can estimate the conditional probability from our our estimates of P(g A obs) and P(obs), (namely Pf:ob· and Pft') by noticing that each element of the queue can go towards implying obs A -.g, obs A g or -.obs. We can easily prove the inequality:

Lemma 4.2

It can be proved that P(globs) has the following bound:

Theorem 4.3

For a proof see Appendix A.

If we choose the midpoint as an estimate, the maxi­ mum error is

What is interesting about this is that the error is in­ dependent of g. Thus when we are generating possible worlds for some observation, and want to have poste­ rior estimates within some error, we can generate the required possible worlds independently of the proposi­ tion that we want to compute the probability of.

5 Search Strategies

The above analysis was independent of the search strat­ egy (i.e., independent of which element we remove from the queue).

We can carry out various search strategies, to enu­ merate the most likely possible worlds. For example, [Poole, 1993] discusses a multiplicative version of A · (Pearl, 1984], with conflicts used to refine a heuristic function.

In this paper the search strategy we use is where the element of the queue with highest prior probability is chosen at each stage. This paper does not study various search strategies, but analyses one. I make no claim that this is the best search strategy (see Section 7), but it forms a base case with which to compare other strategies.

6 Complexity

The problem of finding the posterior probability of a proposition in a Bayesian network is NP hard [Cooper, 1990]. Thus we should not expect that our algorithms will be good in the worst case. Our algorithm, when run to completion, is exponential in computing the exact prior and posterior probability of a hypothesis.

Because of the the "anytime" nature of our algorithm, which trades search time for accuracy, we should not consider run time independently of error. It is interesting to estimate how long it takes on average to get within some error, or how accurate we can expect (or guarantee as as asymptotic behaviour) to be within a certain run time.

As we have probabilities it is possible to carry out an average case complexity analysis of our algorithm.

If we make no assumptions about the probability dis­ tributions, the average case of finding the most likely explanation or prior probability within error f < 1 is ex­ ponential in the size n of the Bayesian network [Provan, 1988]. This can be seen by noticing that the size of com­ plete descriptions are linear in n and so the probability of explanations is exponentially small in n. This means that we need to consider exponentially many explana­ tions to cover any fixed proportion of the probability mass.

This is not always a reasonable distribution assump­ tion, for example, when using this for diagnosis of a sys­ tem that basically works we would like to assume that the underlying distribution is such that there is one as­ signment of values (the "normal values" ) that dominates the probability mass. This also may be appropriate for commonsense reasoning tasks where normality assump­ tions (allegedly) dominate (i.e., we assume that abnor­ mality is rare [McCarthy, 1986]).

For our analysis we assume that we have extreme prob­ abilities for each conditional probability given. For each value of the parents of variable X;, we assume that one of the values for X; is close to one, and the other values are thus close to zero. Those that are close to one we call normality values; those that are close to zero we call faults:

Definition 6.1 If we have extreme probabilities (i.e., all of the conditional probabilities are close to one or zero), then the faults of partial description (v1, · · ·, vk) are those v; such that P(X; = v;jX; 1 = v; 1 1\ · ··I\ X;~ . = v;~,) is close to zero. '

Which values are faults is thus context dependent and depends on the values of the predecessor variables. The term "fault" is derived from its use in diagnosis, where deviation from normality is a fault. The non-monotonic reasoning community has called such things abnormali­ ties [McCarthy, 1986].

Each element of the queue consists of a set of normality conditions and a set of faults. For each set of faults there can be at most one element of Q or W whose faults are exactly that set. A k-fault partial description is one with exactly k faults.

For our complexity analysis, let p (p ~ 1) be the prob­ ability of the least likely normality condition. We can interpret our analysis as being for the case where all of the normality conditions have value p or as a bound for the cases where all normality conditions are greater than p. Let f = 1 - p. f is a bound on the probability of the faults.

Let b + 1 be a bound on the number of values of a variable. Thus b is the maximum number of fault values of any variable given its predecessors. We assume that b is fixed (it is not a function of the number of variables) .

6.1 Error convergence

Suppose we have n variables in our Bayesian network. Consider the case where we are about to choose the first k-fault partial description from the queue. We assume for this analysis that k « n.

At this time there can be no (k + 1)-fault node on the queue (as we can only generate (k + 1)-fault nodes from k-fault nodes and we have not yet expanded out first k-fault node), so there are at most ( b; ) possible

combinations of faults on the priority queue (choosing k faults from the set of bn fault assumptions). Each partial description on the queue has a probability less than Jk

(not because they all have k faults, but because some­ thing with probability less than fk has been chosen from the priority queue and it was the element with highest probability). Thus the probability mass in the queue can

be bounded by ( b: ) fk. So

Thus, we have convergence (as the computation pro­ ceeds and k increases) when bnf < 1.

In order to interpret the value of nf, we use the fol­ lowing lemma:

Lemma 6.2 If nf is small then nf ~ 1 - pn.

pn is the probability that there are no faults. Thus 1-pn is the prior probability of some fault in the system. This is low when we are diagnosing a system which works most of the time, and has rare errors. We will carry out our analysis for the class of systems where the error rate of the total system is low. That is where 6 = nf is low (at least less than t). The analysis is parametrized by 6. This becomes the class of Bayesian networks that we analyse. N.B., with 6 fixed, f --+ 0 as n --+ oo.

6.2 How fast is convergence?

The time to compute the posterior probability of any formula (of bounded size) given some observations using this search algorithm is dominated by the search. This is because the search already includes checking every gen­ erated possible world. Thus we just consider generating the possible worlds.

To see how fast we reach convergence, consider how long it takes to reach the stage where we are about to choose the first k-fault hypothesis from the queue. We assume for this analysis that k « n.

There are at most (bn+1)k ways to generate k or fewer faults. For each of these combinations of faults, we go through the while-loop of Figure 1, at most n times (for the other normality assumptions). The only thing in the while loop that depends on the size of n, is adding and removing elements from the priority queue, which can be done in log IQI time.

(for fixed band k). Thus we can reach the stage where we are taking the first k-fault hypothesis off the queue in O(nk+ 1 Iogn) time.

We can combine the above results on the queue size and time to reach the stage where where we are taking the first k-fault hypothesis off the queue to give us:

Lemma 6.3 For fixed b and k, if 6 < t we can attain an accuracy of (bZ( in the estimate of prior probabilities in O(nk+ 1 logn) time.

For example, if b = 1 (there is only one fault mode) we can attain an accuracy of 6 in 0( n 2 log n) time. We can attain an accuracy of 6 ; in 0( n 3 log n) time. We can attain an accuracy of 6 : in 0( n 4 log n) time.

Theorem 6.4 We can obtain an accuracy of f in time

See Appendix A for a proof of this theorem.

6.3 Posterior Probabilities

As discussed in Section 4.2, at any stage through the loop in the algorithm of Figure 1, we can estimate the posterior probability of g given obs with an error that is independent of g.

To make the posterior error less than f, we require

which can be ensured if we make sure that

We can use the analysis for the prior probability, but multiplying the error bound by a factor that is an esti­ mate of P(obs). As it is unlikely that the observations have a low probability, it is unlikely to have a situation where the error term required is dominated by the prob­ ability of the observation. This observation is reflected in Theorem 6.5 below.

The following theorem gives a PAC (probably approx­ imately correct) characterization of the complexity 4 ·

Theorem 6.5 In the space of all systems, to compute the posterior probability of any proposition (of bounded size) given observation obs, we can guarantee an error of less than f ( f < t), for at least 1 - 1/J of the c~es in time

See Appendix A for a proof of this theorem.

Note that in this theorem we are considering "the space of all systems". For diagnosis, this means that we consider a random artifact. Most of these have no faults; and presumably would not be the subject of di­ agnosis. Thus the space of all systems is probably not the space of all systems that we are likely to encounter is a diagnostic situation. A more realistic space of sys­ tems by which to judge our average-time behaviour is the space of all broken systems, that is those that have at least one fault 5 · We are thus excluding all but 6 of

4 This has the extra property that we know when we are in a case for which we cannot guarantee the error; when we have run our algorithm we know our error.

5 1t could also be argued that this is also inappropriate; we would rather consider the space of systems that exhibit faulty behaviour. This would be much harder to analyse here, as we have no notion of the observable variables developed in this paper. The space of broken devices seems like a reasonable approximation.

this occurs when the systems. To exclude all but 1/J of the broken systems we have to exclude all but 61/; of the total number of systems. We thus have the following corollary.

Corollary 6.6 In the space of all broken systems (with at least one fault), to compute the posterior probability of any proposition (of bounded size) given observation obs, we can guarantee an error of less than f ( f < ~ ), for at least 1 - 1/J of the cases in time

7 Refinements

There are a number of refinements that can be carried out to the algorithm of Figure 1. Some of these are straightforward, and work well. The most straightfor­ ward refinements are:

If we are trying to determine the value of P( a), we can stop enumerating the partial descriptions once it can be determined whether a is true in that partial description. When conditioning on an observation we can prune any partial description that is inconsistent with the observa­ tion.

We do not really require that we find the most likely possible worlds in order, as we are just summing over them anyway. One way to improve the algorithm is to carry out .a depth-first depth-bounded search. We can guess the probability of the least likely possible world we will need to generate, use this as a threshold and carry out a depth-first search pruning any partial description with probability less that this threshold. If the answer is not accurate enough, we decrease the threshold and try again. This is reminiscent of iterative deepening A · [Korf, 1985], but we can decrease the bound in larger ratios as we do not have to find the most likely possible world.

We can use conflicts [de Kleer, 1991] to form a heuris­ tic function for a multiplicative version of A* [Poole, 1993].

See [Poole, 1993] for a Prolog implementation that incorporates these refinements.

8 Comparison with other systems

The branch and bound search is very similar to the candidate enumeration of de Kleer's focusing mecha­ nism [de Kleer, 1991]. This similarity to a single step in de Kleer's efficient method indicates the potential of the search method. He has also been considering cir­ cuits with thousands of components, which correspond to Bayesian networks with thousands of nodes. It seems to be very promising to combine the pragmatic efficiency issues confronted by de Kleer, with the Bayesian network representation, and the error bounds obtained used in this paper.

Poole [1992a] has proposed a Prolog-like search ap­ proach that can be seen as a top-down variant of the bottom-up algorithm presented here. It is not as efficient as the one here. Even if we consider finding the single most normal world, the algorithm here corresponds to forward chaining on definite clauses (see [Poole, 1992b]), which can be done in linear time, but backward chain­ ing has to search and takes potentially exponential time. The backward chaining approach seems more suitable however when we have a richer language [Poole, 1992b].

D'Ambrosio [1992] has a backward chaining search al­ gorithm for "incremental term computation", where he has concentrated on saving and not recomputing shared structure in the search. This seems to be a very promis­ ing approach for when we do not have as extreme prob­ abilities as we have assumed in this paper.

Shimony and Charniak [1990] have an algorithm that is a backward chaining approach to finding the most likely possible world. The algorithm is not as simple as the one presented here, and has worse asymptotic be­ haviour (as it is a top-down approach- see above). It has not been used to find prior or posterior probabilities, nor has the average-case complexity been investigated.

This paper should be seen as a dual to the TOP-N algorithm of Henrion [1991]. We have a different niche. We take no account of the noisy-OR distribution that Henrion concentrates on.

9 Conclusion

This paper has considered a simple search strategy for computing prior and posterior probabilities in Bayesian networks. It is a general purpose algorithm, that is al­ ways correct, and has a niche where it works very well. We have characterised this niche, and have given bounds on how badly it can be expected to perform on average. How common this niche is, is, of course, an open ques­ tion, but the work in diagnosis and nonmonotonic rea­ soning would suggest that reasoning about normality is a common task.

A Proofs

Theorem 4.3

Proof: Consider what happens to the elements of the queue. Let a be the proportion of the possible worlds that are descendents of elements of the priority queue in which obs 1\ -,g is true. Let (3 be the proportion in which obs 1\ g is true. Then a + (3 is the proportion in which obs is true.

As all of the possible worlds are either in W or are descendents of elements of the priority queue, we have

We want to maximise this formula under the constraints that 0 ~ a ~ 1, 0 ~ (3 ~ 1 and 0 ~ a + (3 ~ 1. There are no internal extrema in this formula, and so the maxima occur at the extremes. There are a = (3 = 0, a = 1/\ (3 = 0 and a = 1/\ (3 = 0. These correspond to the three values in Lemma 4.2, and the theorem follows directly from Lemma 4.2. 0

Theorem 6.4 We can obtain an accuracy off in time

Proof: If we require an accuracy off, we ensure (b:( < f. This can be obtained if we ensure (b6)k = f. Solving for k, we get

This requires 0( nk+ 1 log n) time. Substituting the value for k, the theorem follows. 0

Theorem 6.5 In the space of all systems, to compute the posterior probability of any proposition (of bounded size) given observation obs, we can guarantee an error of less than f ( f < ! ), at least 1 - '1/J of the cases in time

Proof: We can assume that the P(obs) ~ '1/J. This will be wrong less than '1/J of the cases (by definition). If we make sure that PQ < f, then

To achieve error f, we make sure PQ < £'1/J. Then, as £'1/J < 2£Pfj 8 , we will have PQ < 2£Pfj 8 , which, as described above implies that the error will be less than f. Thus we have to ensure that PQ < min(f,t:t/J) = £'1/J. By Theorem 6.4, this can be done in time

0

Acknowledgements

Thanks to Andrew Csinger, Michael Horsch, Runping Qi and Nevin Zhang for valuable comments on this pa­ per. This research was supported under NSERC grant OGP0044121, and under Project B5 of the Institute for Robotics and Intelligent Systems.