5.2 Variable Elimination (VE; the Elimination Algorithm)
It might seem that this is merely the price of inference, but a moment’s reflection will show how naïve our first approach was. The distributive property of multiplication (over addition) licenses swapping the order of multiplication and summation, which obviates the need for a single table for the joint distribution. For example, to compute in Fig. 5.1a, we need to compute
In code, this would amount to
# non-naive way to compute p(X8)
rhostar123 = p_X3givenX1X2*pX1*pX2
assert rhostar123.shape == (2, 2, 2, 1, 1, 1, 1, 1)
phistar12345 = rhostar123.sum(axis=(0,1), keepdims=True)
rhostar345 = phistar12345*p_X4givenX3*p_X5givenX3
assert rhostar345.shape == (1, 1, 2, 2, 2, 1, 1, 1)
phistar3456 = rhostar345.sum(axis=2, keepdims=True)
rhostar456 = phistar3456*p_X6givenX4
assert rhostar456.shape == (1, 1, 1, 2, 2, 2, 1, 1)
phistar4567 = rhostar456.sum(axis=3, keepdims=True)
rhostar567 = phistar4567*p_X7givenX5
assert rhostar567.shape == (1, 1, 1, 1, 2, 2, 2, 1)
phistar5678 = rhostar567.sum(axis=4, keepdims=True)
rhostar678 = phistar5678*p_X8givenX6X7
assert rhostar678.shape == (1, 1, 1, 1, 1, 2, 2, 2)
marginal = rhostar678.sum(axis=(5,6), keepdims=True)
Evidently, it was not necessary to create—or collapse—any tables with more than three ways, so the computational complexity will be of order rather than .
Does the use of the distributive property alone guarantee the cheapest computational complexity? No: in the computation just given, we implicitly made a choice about the order in which variables would be eliminated from the problem, which turned out to be optimal. Suboptimal orderings exist; for example, choosing to eliminate first would result in a five-way table:
Worse still, finding the optimal ordering for an arbitrary graph is NP hard.
On the other hand, we do have strong intuitions about good elimination orderings. To remove a variable, we must first multiply together all “potential” functions that make reference to it, whether they are the conditional and marginal distributions parameterizing the graph, or the intermediate potentials created when eliminating variables (referred to as rhostar* in the code snippets above). It is best to remove nodes whose potentials make the fewest references to other variables. That is why it is better to remove and before removing .
Making observations.
Suppose that we are after, not , but . (Think of this as nevertheless a single row or column of a table, i.e. with just two entries if is binary; whereas would be a full two-by-two table.) Evidently, we could compute this conditional simply by refraining from marginalizing out until the end. That is, we could use the elimination ordering , yielding . At this point we could marginalize out , and use to normalize . Finally, we would simply extract the row of this table corresponding to .
This is almost the right approach, but notice that refraining from marginalizing out has increased the size of the largest table from 3 (in our example of unconditional inference over ) to 4 dimensions. This is particularly vexing since, in the final table, we ultimately throw out all but one of the rows in the dimension.
In fact, we can reduce the maximum table size back down to 3 if we simply throw away the irrelevant parts of the table first. That this is possible should be somewhat intuitive. More formally, it requires two fact: (1) Extraction commutes with normalization: we can extract out the relevant row of the pairwise marginal , namely , and then normalize. (2) Extraction commutes with the sums and products of the rest of the algorithm. This can be seen by dressing extraction up like marginalization—and indeed this will also allow us to write the algorithm in a convenient form:
or, in code, as
Under this formulation, acts just like any other potential function—in particular, it looks like a single row—and extracting out the appropriate row of the table looks just like marginalization. The same considerations (discussed above) about the distributive property of multiplication and about elimination orderings apply just as well to the delta function. In the present example, it would behoove us to remove first, i.e. to use the elimination ordering , rather than to wait until the end.
Thus the only additions to our initial version of variable elimination††margin: variable elimination are the pre-processing step of augmenting our list of potentials with delta functions for all observations we wish to make; and a post-processing step of normalization. And in fact, normalizing singleton marginals is not expensive and is harmless (has no effect) when computing unconditional probabilities, so we could simply make this a standard post-processing step of the algorithm. Then a single algorithm can be used for both conditional and unconditional inference (see the pseudo-code in Algorithm 1).
I emphasize again that what we have assimilated to marginalization is just one of the two parts of the process “conditioning on an observation”—namely, extracting a row from the table. We could call this making an observation (rather than conditioning on one.) The second step is renormalization, which transforms a joint distribution into a conditional distribution. We have relegated this to a post-processing step. But note that if we were interested in the full two-by-two table , we would only normalize and skip row extraction altogether.
Finally, note that VE as described works only for singleton marginals and conditionals, although extending it to marginals over neighbors in the graph is straightforward.
def VariableElimination(graph, observations, target_id):
initialize(graph, target_id)
observe(observations)
update(graph)
target_prob = normalize(target_id)
def initialize(graph, target_id):
// choose ordering with target_id last
// put each vertex’s potential on active list
active_list = graph.vertices
def observe(observations):
// the observations have delta potentials associated w/them
active_list += observations
def update(graph):
for id in ordering do
rhostar = 1
for entry in active_list do
if id in entry.variables then
rhostar *= entry.potential
active_list.remove(entry)
phistar = rhostar.sum(axis=id)
active_list.append(phistar)
def normalize(target_id):
target_vertex = find(graph.vertices, target_id)
Z = target_vertex.potential.sum(axis=target_id)
return target_vertex.potential/Z
Re-using computations.
A related question is: if we can’t compute arbitrary (say) pairwise marginals, like , can we at least reuse some of the computations for computing one singleton marginal, , in order to compute another, ?
The short answer is that yes, we can, but how to do so optimally on an arbitrary graph is not obvious. We will therefore devise, first, such a method only for trees; and then we will derive the full-blown method for arbitrary graphs. The first algorithm, belief propagation, is not merely a diversion; we will use ideas from it in order to derive the second algorithm, the junction-tree algorithm.
Before deriving either, however, let us generate some intuition by considering how we might approach this problem for small graphs—first, Fig. 5.1a, the directed acyclic graph we have been working with. Consider again our computation of above under the elimination ordering , and in particular the intermediate potentials in the code snippet:
Evidently, each intermediate potential is also the marginal probability of the elimination clique††margin: elimination clique , i.e. the set of variables that must be represented together in one table in order to eliminate node/variable . (In order to conclude this, we used some of the independence properties of the graph.) Therefore, to compute the marginal probability of any node in this graph, we could simply eliminate in the order , making sure to retain the “potential functions” over the elimination cliques. Then, finally, we (further) marginalize these clique potentials as needed to get singleton marginals.
The seeds of the junction-tree algorithm.
We can generate more intuition by considering the clique graph††margin: clique graph in Fig. 5.1b. In this graph, each elimination clique created has been assigned its own node. Between each clique node††margin: clique node , a separator node††margin: separator node has been inserted, corresponding to the variables that are shared between its two neighbors in the graph. Now let us assign each to a clique node, and each to a separator node. Then the sequence of computations (5.2) can be re-described as the following operations on the clique graph.
First, assign each of the potentials (conditional distributions) of the original graph, Fig. 5.1a, to exactly one clique node. Any assignment is allowed so long as all of the variables referenced by the potential are among the clique variables. Then each clique potential††margin: clique potential is defined to be the product of the potentials assigned to it. For example, we will assign
Second, assign a separator potential††margin: separator potential to each of the separator nodes; in particular, let for all .22 2 Two (or more) separator nodes can have the same separator sets. Therefore, to avoid ambiguity, we subscript each with all of its neighbors’ indices.
Note well that under this assignment, the product of all the potentials on this clique tree is, as for a directed graphical model, the joint distribution! More precisely, we note that
This is a desirable property and we will try to design our updates so that Eq. 5.2 is invariant to them. In particular, we update the chain in a single “inward” pass as follows. Starting at the “leaf” (1-2-3), marginalize out all variables that are not among the separator variables—in this case, both and , since only is in the separator. Update the separator potential to this marginalized potential:
Next, to account for this change in the separator potential, i.e. to cancel it out of the denominator in Eq. 5.2, multiply it onto the next clique potential.
The second equality follows from the fact that the separator potentials are initialized to unity, , and establishes a symmetry with the updates given below.
Continuing according to this update rule evidently corresponds to the sequence of computations in (5.2), i.e. to a single run of VE on Fig. 5.1a. At the end of this run both the clique potentials and the separator potentials are the marginal distributions over their respective local variables, and furthermore the joint distribution is still given by the ratio of products in Eq. 5.2.
Unnormmalized parameterizations.
Having said that, this method without modification clearly will not suffice in the presence of observations. For example, if we have observed , then computing the conditional probability will still require the “inward” or forward pass through the graph, but computing will require a full outward or backward pass, in order for the observation to propagate from back to .
Let us work through this example. The final step in the forward pass would become
This intermediate potential deserves the double- as well as the single-starred name because it is the first computation of the outward pass as well as the last computation of the inward pass. Now we update the separator potential in the same way as in the inward pass, namely: marginalize out of the clique potential the variables not in the separator node, and set the separator potential equal to the result:
Note that this is a two-way table, in particular, the unnormalized conditional probability over conditioned on . This can be used to enforce consistency with the previous clique potential, by updating it to
The updated potential remains a three-way table, now the unnormalized probability of conditioned on . (Make sure to convince yourself that the second equality holds.)
It is probably clear at this point that this “outward pass” can be continued in the same fashion back to , transforming all the clique marginals into the corresponding unnormalized conditional probabilities. And the outward pass is essentially identical to the inward pass; cf. Eqs. 5.3 and 5.5 and Eqs. 5.4 and 5.6. Indeed, we can summarize the updates in either direction as
(Eq. 5.8 should once again be thought of as broadcast multiplication and division.) What may be less obvious is that, with the addition of a message-passing protocol, this procedure generalizes to clique trees, rather than just the “clique chain” of this example. This is the junction-tree algorithm, or more precisely an implementation known as the Hugin algorithm††margin: Hugin algorithm . It also works without further modification for undirected graphs. Indeed, an undirected graph is essentially like a directed graph with observations, since both are parameterized by an unnormalized joint distribution. Unnormalized joints should therefore be thought of as the generic parameterization, with normalized joint distributions, as in our original problem of computing in Fig. 5.1a, a special case. Corresponding to its specialness is the fact that only a single pass through the network is required to compute all the clique marginals.