9. Exact Inference

Inference tasks

A graphical model (chapter 8) is useful only if questions can be answered from it. For a distribution \(P(x)=\frac1Z\prod_af_a(x_a)\) given as a product of factors, which covers Bayesian networks (\(Z=1\), one factor per CPT) and Markov networks alike, the standard tasks are:

• Marginals and conditionals: \(P(X\mid e)\) for a query variable \(X\) and evidence \(e\). Diagnosis, filtering, and prediction are all of this form.
• The partition function \(Z\), or the probability of the evidence \(P(e)\), which is the same computation with the evidence variables clamped. It is needed for learning and for comparing models.
• MAP inference, or the most probable explanation: \(\arg\max_xP(x\mid e)\), the single most likely joint assignment of all non-evidence variables. A marginal MAP query maximizes over some variables and sums over others, and is harder than both.

All three are hard in general. Computing marginals in Bayesian networks is #P-hard (Cooper, 1990), since counting the satisfying assignments of a formula reduces to it (Appendix B), and even approximating a conditional probability to within a constant relative error is NP-hard (Dagum and Luby, 1993). MAP inference is NP-hard, as it contains satisfiability. What makes inference tractable is structure: the algorithms of this chapter run in time exponential only in a measure of how far the graph is from a tree, and linear in its size. When that measure is too large, chapter 10 gives up exactness.

Variable elimination

Pushing sums inward

Inference by enumeration sums the full joint over the hidden variables. For the burglary network,

\[P(B\mid j,m)=\alpha\sum_e\sum_aP(B)\,P(e)\,P(a\mid B,e)\,P(j\mid a)\,P(m\mid a),\]

a sum of \(2\times2\) terms, each a product of five numbers, for each value of \(B\). Enumeration recomputes the same subexpressions many times: \(P(j\mid a)P(m\mid a)\) is the same for every value of \(e\). Moving each sum as far right as it goes,

\[P(B\mid j,m)=\alpha\,P(B)\sum_eP(e)\sum_aP(a\mid B,e)\,P(j\mid a)\,P(m\mid a),\]

and evaluating from the inside out, storing each intermediate result, removes the repetition. This is dynamic programming, and on a network with \(n\) variables it can turn an \(O(2^n)\) sum into one that is linear in \(n\).

Factors and their operations

Variable elimination carries out this idea on tables. A factor is a function of a set of variables, stored as a table with one entry per joint value. Three operations suffice:

• restriction fixes an evidence variable to its observed value, removing it from the factor;
• product of two factors is a factor on the union of their variables, \(h(x,y,z)=f(x,y)\,g(y,z)\), pointwise, like a database join;
• summing out a variable from a factor, \(h(x)=\sum_yf(x,y)\), marginalizes it.

The algorithm restricts every factor to the evidence, then eliminates the hidden variables one at a time in some elimination order: to eliminate \(Y\), multiply all factors that mention \(Y\) and sum \(Y\) out of the product, replacing those factors by the result. When only the query variable remains, the product of the remaining factors, normalized, is the answer. Two refinements matter in practice. Variables that are neither ancestors of the query nor of the evidence (barren nodes) can be removed before starting, since they sum to one. And with several queries on the same evidence, the work of elimination can be shared, which leads to the message-passing algorithms below.

import string
from itertools import product

import numpy as np


class Factor:
    """A table over named binary variables: vars is a tuple of names, table has one axis per variable."""

    def __init__(self, vars, table):
        self.vars, self.table = tuple(vars), np.asarray(table, float)

    def __mul__(self, other):
        vs = self.vars + tuple(v for v in other.vars if v not in self.vars)
        L = {v: string.ascii_letters[i] for i, v in enumerate(vs)}
        spec = "".join(L[v] for v in self.vars) + "," + "".join(L[v] for v in other.vars) + "->" + "".join(L[v] for v in vs)
        return Factor(vs, np.einsum(spec, self.table, other.table))

    def sum_out(self, v):
        return Factor([u for u in self.vars if u != v], self.table.sum(axis=self.vars.index(v)))

    def restrict(self, v, value):
        return Factor([u for u in self.vars if u != v], np.take(self.table, int(value), axis=self.vars.index(v)))


def cpt(child, parents, p_true):
    """Factor for P(child | parents) from P(child = 1 | parents), given as an array over the parents' values."""
    p = np.asarray(p_true, float)
    return Factor(tuple(parents) + (child,), np.stack([1 - p, p], axis=-1))


# The burglary network; index 1 means true.
factors = [cpt("B", [], 0.001), cpt("E", [], 0.002),
           cpt("A", ["B", "E"], [[0.001, 0.29], [0.94, 0.95]]),
           cpt("J", ["A"], [0.05, 0.90]), cpt("M", ["A"], [0.01, 0.70])]


def eliminate(factors, query, evidence, order):
    """Variable elimination: restrict to the evidence, then sum out the hidden variables one at a time."""
    fs = []
    for f in factors:
        for v, val in evidence.items():
            if v in f.vars:
                f = f.restrict(v, val)
        fs.append(f)
    largest = 0
    for v in order:
        involved = [f for f in fs if v in f.vars]
        prod = involved[0]
        for f in involved[1:]:
            prod = prod * f
        largest = max(largest, prod.table.size)
        fs = [f for f in fs if v not in f.vars] + [prod.sum_out(v)]
    result = fs[0]
    for f in fs[1:]:
        result = result * f
    table = result.table / result.table.sum()
    return table[1], largest


p, largest = eliminate(factors, "B", {"J": 1, "M": 1}, ["E", "A"])
print(f"P(burglary | John and Mary call) = {p:.4f}; largest intermediate factor has {largest} entries")

# The same query on a chain X1 -> X2 -> ... -> Xn with evidence at the end: elimination versus enumeration.
rng = np.random.default_rng(0)
n = 20
chain = [cpt("X1", [], 0.3)] + [cpt(f"X{i}", [f"X{i - 1}"], [rng.uniform(0.01, 0.05), rng.uniform(0.95, 0.99)])
                                for i in range(2, n + 1)]
p_ve, largest = eliminate(chain, "X1", {f"X{n}": 1}, [f"X{i}" for i in range(2, n)])
num = den = 0.0
for x in product([0, 1], repeat=n):                         # brute force over all 2^20 assignments
    w = 0.3 if x[0] else 0.7
    for i in range(1, n):
        q = chain[i].table[x[i - 1], 1]
        w *= q if x[i] else 1 - q
    if x[-1] == 1:
        den += w
        num += w * x[0]
print(f"chain of {n}: P(X1 | X{n} = 1) = {p_ve:.6f} by elimination (largest factor {largest} entries),"
      f" {num / den:.6f} by summing {2 ** n:,} terms")
# P(burglary | John and Mary call) = 0.2842; largest intermediate factor has 8 entries
# chain of 20: P(X1 | X20 = 1) = 0.428987 by elimination (largest factor 8 entries), 0.428987 by summing 1,048,576 terms

Elimination reproduces the burglary posterior of chapter 8 with no factor larger than eight entries. On a chain of twenty variables it agrees with the brute-force sum over a million assignments while never building a table with more than three variables; its cost grows linearly with the length of the chain, the brute force doubles with each variable.

Elimination order and treewidth

The cost of variable elimination is dominated by the largest factor it creates, and that depends on the order. Eliminating a variable connects all its neighbors, because their factors are multiplied together, so the process can be simulated on the undirected graph of the model (the moral graph, for a Bayesian network): eliminating a node joins its neighbors into a clique, adding fill-in edges, and removes it. The induced width of an order is the largest number of neighbors a node has when it is eliminated; the largest factor has one more variable than that, and with domains of size \(d\) the total cost is \(O(n\,d^{w+1})\) for induced width \(w\). The minimum induced width over all orders is the treewidth of the graph, the same quantity that governs the tree decompositions of chapter 3 (Appendix C).

• Trees have treewidth 1, so inference on them is linear. Singly connected networks, polytrees, in which there is at most one undirected path between any two nodes, have treewidth equal to the largest number of parents, and inference is linear in the total size of the CPTs.
• An \(n\times n\) grid has treewidth \(n\): inference on it costs \(O(d^{n+1})\), exponential in the side length and not in the number of variables \(n^2\).
• Finding an optimal order is NP-hard, and greedy heuristics are used: min-degree eliminates the node with fewest neighbors, and min-fill the one whose elimination adds the fewest edges. They are fast and usually good, but not optimal.

Left: the induced width of four elimination orders on \(n\times n\) grids. Eliminating row by row achieves the treewidth, \(n\); the greedy min-degree and min-fill heuristics, with ties broken by node name, are optimal up to side 6 but reach width 21 on the \(15\times15\) grid; a random order reaches 73, which would need tables of \(2^{74}\) entries for binary variables. Right: the total size of the factors created by variable elimination with the min-fill order, for random Bayesian networks in which each node has zero to three parents chosen at random among earlier nodes, against the \(n2^n\) operations of enumeration. The median induced width grows from 3 at 10 variables to 42 at 320, so elimination also becomes exponential, but at 320 variables it needs about \(2\times10^{13}\) entries against \(7\times10^{98}\).

Random networks with long-range connections have large treewidth, and exact inference on them is hopeless beyond a few hundred variables. Real networks built by experts, and the grids, chains, and trees of vision and speech, usually have more structure, and exact methods handle many of them.

Message passing on trees

Sum-product

Variable elimination answers one query. On a tree, all marginals can be computed at twice the cost of one, by keeping the intermediate factors as messages. For a pairwise Markov network \(P(x)\propto\prod_i\phi_i(x_i)\prod_{(i,j)}\psi_{ij}(x_i,x_j)\) on a tree, the message from node \(i\) to its neighbor \(j\) is

\[m_{i\to j}(x_j)=\sum_{x_i}\phi_i(x_i)\,\psi_{ij}(x_i,x_j)\prod_{k\in N(i)\setminus j}m_{k\to i}(x_i),\]

the result of eliminating the whole subtree on \(i\)’s side of the edge, summarized as a function of \(x_j\). The marginal of each node is proportional to its own potential times all incoming messages:

\[P(x_i)\propto\phi_i(x_i)\prod_{k\in N(i)}m_{k\to i}(x_i).\]

A message can be sent once its sender has heard from all its other neighbors. Choosing a root, an upward pass from the leaves and a downward pass from the root compute all \(2(n-1)\) messages, after which every marginal is available (Appendix A). On a factor graph (chapter 8), the same algorithm alternates two kinds of messages, variable-to-factor (the product of the other incoming messages) and factor-to-variable (the factor times the incoming messages, summed over all its other variables), and it is known as the sum-product algorithm (Kschischang, Frey, and Loeliger, 2001). The forward–backward algorithm for hidden Markov models (chapter 11) and the Kalman smoother are sum-product on a chain.

Max-product and MAP inference

Replacing the sum in the message by a maximum gives max-product, which computes max-marginals: for each node and value, the probability of the best joint assignment in which that node takes that value. With ties broken consistently, taking the best value of each node’s max-marginal gives the most probable joint state; alternatively, the upward pass followed by backtracking from the root, with each node choosing the value that achieved the maximum given its parent’s choice, recovers the MAP state as in dynamic programming. In log space, max-product becomes max-sum, and on a chain it is the Viterbi algorithm of chapter 11.

The MAP state is not the state of individually most probable values, because the variables interact:

from itertools import product

import numpy as np

rng = np.random.default_rng(1)
n, k = 10, 3                                           # a tree of 10 variables with 3 states each
parent = [None] + [int(rng.integers(i)) for i in range(1, n)]
edges = [(parent[i], i) for i in range(1, n)]
unary = rng.uniform(0.5, 2.0, (n, k))                   # node potentials
pair = {e: rng.uniform(0.2, 3.0, (k, k)) for e in edges}   # edge potentials, indexed [x_parent, x_child]
nbrs = {i: [] for i in range(n)}
for a, b in edges:
    nbrs[a].append(b)
    nbrs[b].append(a)


def psi(i, j):
    """Edge potential as a matrix indexed [x_i, x_j]."""
    return pair[(i, j)] if (i, j) in pair else pair[(j, i)].T


def messages(combine):
    """Messages m[i -> j](x_j) for every directed edge, computed leaves-first then root-first.
    combine = np.sum gives sum-product, np.max gives max-product."""
    m = {}

    def send(i, j):
        incoming = unary[i].copy()
        for u in nbrs[i]:
            if u != j:
                incoming *= m[(u, i)]
        msg = combine(incoming[:, None] * psi(i, j), axis=0)
        m[(i, j)] = msg / msg.sum()                    # normalize for numerical safety

    order = []                                         # a DFS order from the root 0
    stack = [0]
    while stack:
        v = stack.pop()
        order.append(v)
        stack.extend(c for c in nbrs[v] if c != parent[v])
    for v in reversed(order[1:]):                      # upward pass: children before parents
        send(v, parent[v])
    for v in order[1:]:                                # downward pass: parents before children
        send(parent[v], v)
    return m


def beliefs(m):
    b = unary.copy()
    for (i, j), msg in m.items():
        b[j] *= msg
    return b / b.sum(axis=1, keepdims=True)


marg = beliefs(messages(np.sum))

# Brute force over all 3^10 joint states.
states = np.array(list(product(range(k), repeat=n)))
w = np.prod(unary[np.arange(n), states], axis=1)
for a, b in edges:
    w *= pair[(a, b)][states[:, a], states[:, b]]
exact = np.array([[w[states[:, i] == s].sum() for s in range(k)] for i in range(n)]) / w.sum()
print(f"sum-product marginals agree with brute force over {len(states):,} states:", np.allclose(marg, exact, atol=1e-12))
print("marginal of variable 0:", np.round(marg[0], 4))

# Max-product: each variable's max-marginal peaks at its value in the most probable joint state (ties aside).
maxmarg = beliefs(messages(np.max))
x_map = maxmarg.argmax(axis=1)
print("MAP state by max-product:", x_map.tolist())
print("MAP state by brute force:", states[w.argmax()].tolist())
x_marg = marg.argmax(axis=1)
print("most probable value of each variable separately:", x_marg.tolist())
print(f"probability of the MAP state {w.max() / w.sum():.5f}; of the state of separate maxima"
      f" {w[np.all(states == x_marg, axis=1)][0] / w.sum():.5f}")
# sum-product marginals agree with brute force over 59,049 states: True
# marginal of variable 0: [0.2499 0.1875 0.5625]
# MAP state by max-product: [2, 1, 0, 2, 2, 1, 1, 1, 2, 2]
# MAP state by brute force: [2, 1, 0, 2, 2, 1, 1, 1, 2, 2]
# most probable value of each variable separately: [2, 1, 1, 2, 2, 1, 1, 1, 2, 2]
# probability of the MAP state 0.00248; of the state of separate maxima 0.00169

Two passes of 18 messages give all ten marginals exactly. Max-product finds the most probable of the 59,049 joint states, which differs from the state that picks each variable’s most probable value in one variable and is almost 50% more probable. Which answer is right depends on the loss: the MAP state minimizes the probability of getting anything wrong, and the per-variable choice minimizes the expected number of wrong variables, the distinction between sequence-level and symbol-level decoding in speech recognition and error-correcting codes.

The junction tree algorithm

Message passing is exact only on trees, but any graph can be turned into a tree of clusters. A clique tree, or junction tree, is a tree whose nodes are clusters of variables such that every factor lies in some cluster and the running intersection property holds: every variable shared by two clusters is present in all clusters on the path between them. The clusters produced by variable elimination, each eliminated node with its neighbors at that moment, form such a tree when connected appropriately; equivalently, triangulate the moral graph and take its maximal cliques. Neighboring clusters exchange messages over their shared variables, the separators, exactly as in sum-product, and after an upward and a downward pass every cluster holds the joint marginal of its variables (it is calibrated). This is the junction tree algorithm (Lauritzen and Spiegelhalter, 1988), the standard exact method for moderately sized networks. Its cost is exponential in the size of the largest cluster, the treewidth plus one, so it faces the same limits as elimination, and it pays in memory for storing all cluster tables.

Two further ideas extend exact inference. Knowledge compilation turns a network, once, into an arithmetic circuit whose evaluation answers any query in time linear in the circuit (Darwiche, 2003); circuits exploit local structure, such as determinism and context-specific independence, that treewidth ignores, and their derivatives give all marginals at once, backpropagation in another guise. Lifted inference exploits symmetry in relational models, treating interchangeable objects as a group instead of one by one. For Gaussian models, exact inference is linear algebra, the conditioning of Foundations chapter 4, and message passing on sparse Gaussian graphs is a way to solve sparse linear systems.

Appendices