Exact inference in a probabilistic graphical model is, in general, #P-hard#P is the class of counting problems attached to NP decision problems. Roth (1996) showed exact inference is complete for it — so a general polynomial-time algorithm is not expected. — but "in general" hides a lot. For many graphs we can compute marginals efficiently, and the trick is almost embarrassingly simple: exploit the distributive law to push summations inside products. That idea, made systematic, is variable elimination (VE).
Factors and the joint distribution
A discrete graphical model over variables $X_1, \dots, X_n$ factorizes the joint as a normalized product of non-negative factors $\phi_c$, each defined over a subset (scope) $\mathbf{X}_c$ of the variables:
$$ p(x_1, \dots, x_n) = \frac{1}{Z} \prod_{c} \phi_c(\mathbf{x}_c), \qquad Z = \sum_{x_1, \dots, x_n} \prod_c \phi_c(\mathbf{x}_c). $$
A factor is just a table: it maps each joint assignment of its scope to a non-negative number. Two operations on factors are all we need. The product $\phi \cdot \psi$ has scope $\mathrm{scope}(\phi) \cup \mathrm{scope}(\psi)$ and multiplies entries that agree on shared variables. Marginalization (summing out) a variable $B$ from $\phi$ produces a new factor $\sum_B \phi$ whose scope drops $B$.
Swapping in a real image. The figure above is an inline SVG so
the post stays self-contained. To use a raster image instead, drop the file in
blog/images/ and replace the <svg> with
<img src="../images/your-figure.png" alt="…"> — the
<figcaption> styling is identical either way.
The inference task
Suppose we want the marginal $p(C)$ in the chain of Figure 1. Naively we sum the joint over everything else:
$$ p(C) \;\propto\; \sum_{A}\sum_{B} \phi(A,B)\,\phi(B,C). $$
If each variable takes $k$ values, this double sum costs $O(k^3)$ — and for a chain of length $n$ it would be $O(k^n)$. That exponential blow-up is the thing we want to avoid.
Pushing sums inside products
The factor $\phi(A,B)$ does not depend on $C$, and $\phi(B,C)$ does not depend on $A$. So we are free to rearrange the sums:
$$ p(C) \;\propto\; \sum_{B} \phi(B,C) \underbrace{\sum_{A} \phi(A,B)}_{\tau(B)} \;=\; \sum_B \phi(B,C)\,\tau(B). $$
We computed an intermediate factor $\tau(B) = \sum_A \phi(A,B)$ once, in $O(k^2)$, and reused it.This reuse of intermediate results is exactly dynamic programming: we solve each subproblem once and cache the answer instead of recomputing it inside the outer sum. The outer sum is another $O(k^2)$. The exponential collapsed to linear-in-$n$. This is the whole idea; everything else is bookkeeping about which variable to sum out next.
Proof Correctness of one elimination step
Claim. Let $\Phi$ be a set of factors and let $Y$ be a variable appearing in the subset $\Phi_Y \subseteq \Phi$. Define the new factor $\tau = \sum_Y \prod_{\phi \in \Phi_Y} \phi$. Then replacing $\Phi_Y$ by the single factor $\tau$ leaves the marginal over the remaining variables unchanged.
Proof. Partition the factors into those that mention $Y$ (call the set $\Phi_Y$) and those that do not ($\Phi_{\bar Y}$). Because a factor in $\Phi_{\bar Y}$ is constant in $Y$, it commutes past the sum $\sum_Y$:
$$ \sum_{Y} \prod_{\phi \in \Phi} \phi \;=\; \sum_{Y} \Big(\prod_{\psi \in \Phi_{\bar Y}} \psi\Big) \Big(\prod_{\phi \in \Phi_{Y}} \phi\Big) \;=\; \Big(\prod_{\psi \in \Phi_{\bar Y}} \psi\Big) \underbrace{\sum_{Y} \prod_{\phi \in \Phi_{Y}} \phi}_{\displaystyle \tau}. $$
The pull-out step is exactly the distributive law $a(b+c) = ab + ac$ applied to the entries of the tables. The right-hand side is the product of $\Phi_{\bar Y} \cup \{\tau\}$, which is the model after elimination. Summing either side over the remaining variables gives identical results, so all downstream marginals are preserved. $\blacksquare$
The algorithm
Generalizing: fix an ordering of the variables to eliminate. For each variable in turn, gather every factor that mentions it, multiply them together, sum the variable out, and put the resulting factor back into the pool. "Summing out" and "marginalizing" mean the same operation — the two names are used interchangeably.
Algorithm Sum-product variable elimination
| Input: factors $\Phi$, elimination ordering $(Y_1, \dots, Y_m)$, query variables $Q$ | ||
| Output: an unnormalized factor over $Q$ | ||
| 1 | for $i = 1$ to $m$ do | eliminate one variable at a time |
| 2 | $\Phi_i \leftarrow \{\, \phi \in \Phi : Y_i \in \operatorname{scope}(\phi) \,\}$ | factors mentioning $Y_i$ |
| 3 | $\psi \leftarrow \prod_{\phi \in \Phi_i} \phi$ | multiply them together |
| 4 | $\tau \leftarrow \sum_{Y_i} \psi$ | sum $Y_i$ out |
| 5 | $\Phi \leftarrow (\Phi \setminus \Phi_i) \cup \{\tau\}$ | replace with the new factor |
| 6 | end | |
| 7 | return $\prod_{\phi \in \Phi} \phi$ | product of surviving factors |
Normalize the returned factor (divide by its sum) to get an actual probability distribution over $Q$.
Cost and elimination order
Each step's cost is dominated by the size of the largest intermediate factor it builds. If $\psi$ in the algorithm above has scope of size $w$, forming and summing it costs $O(k^{\,w+1})$. The largest $w+1$ encountered over the whole run — the induced width of the ordering — is what governs tractability.The minimum induced width over all orderings is the graph's treewidth. Trees have treewidth 1, which is why message passing on trees is cheap.
The catch: finding the ordering that minimizes induced width is itself NP-hard. In practice we lean on cheap heuristics — min-degree and min-fill being the usual suspects.
A reference implementation
Here is a compact, dependency-light implementation. A factor is a pair
(vars, table) where table is a NumPy array whose axes
line up with vars. Fold it open:
Code variable_elimination.py
import numpy as np
from functools import reduce
class Factor:
"""A factor: an array `table` whose axes are labeled by `vars`."""
def __init__(self, vars, table):
self.vars = list(vars)
self.table = np.asarray(table, dtype=float)
def multiply(self, other):
# align both factors onto the union of their variables, then broadcast-multiply
new_vars = self.vars + [v for v in other.vars if v not in self.vars]
a = self._reshape_to(new_vars)
b = other._reshape_to(new_vars)
return Factor(new_vars, a * b)
def sum_out(self, var):
axis = self.vars.index(var)
table = self.table.sum(axis=axis)
return Factor([v for v in self.vars if v != var], table)
def _reshape_to(self, target_vars):
shape = [self.table.shape[self.vars.index(v)] if v in self.vars else 1
for v in target_vars]
# move existing axes into target order, then pad missing ones with size-1 axes
perm = [self.vars.index(v) for v in target_vars if v in self.vars]
return np.transpose(self.table, perm).reshape(shape)
def variable_elimination(factors, order, query):
factors = list(factors)
for y in order: # eliminate one variable at a time
touching = [f for f in factors if y in f.vars]
factors = [f for f in factors if y not in f.vars]
psi = reduce(lambda a, b: a.multiply(b), touching)
factors.append(psi.sum_out(y)) # sum y out, put result back
result = reduce(lambda a, b: a.multiply(b), factors)
result.table /= result.table.sum() # normalize
return result
# --- chain A - B - C from Figure 1 -------------------------------------------
phi_AB = Factor(["A", "B"], [[0.8, 0.2], [0.3, 0.7]])
phi_BC = Factor(["B", "C"], [[0.6, 0.4], [0.1, 0.9]])
p_C = variable_elimination([phi_AB, phi_BC], order=["A", "B"], query=["C"])
print(p_C.vars, p_C.table) # -> ['C'] [0.3675 0.6325]
The order=["A", "B"] reproduces the hand derivation above:
eliminate $A$ to build $\tau(B)$, then eliminate $B$. Swap in any ordering
and the answer is unchanged — only the cost differs.
That's variable elimination end to end: a factor algebra, one distributive-law step proved correct, an ordering-driven loop, and the induced width as the single number that decides whether any of it is fast. Junction trees, belief propagation, and the whole message-passing family are, at heart, VE run cleverly and cached.