The Max-Product Algorithm

Massachusetts Institute of Technology Department of Electrical Engineering and Computer Science 6.438 Algorithms for Inference Fall 2014 11 The Max-Product Algorithm In the previous lecture, we introduced the MAP elimination algorithm for computing an MAP configuration for any undirected graph. Today, we specialize the algorithm to the case of trees. The resulting algorithm is the Max-Product algorithm, which recovers the max-marginal at each node of a graph over random variables x1 , . . . , xN : pxi (xi ) = max px (x). xj :j�=i If each x∗i ∈ argmaxxi pxi (xi ) is uniquely attained, i.e., there are no ties at each node, then the vector x∗ = (x1∗ , x2∗ , . . . , x∗N ) is the unique global MAP configuration. Oth­ erwise, keeping track of “backpointers” indicating what maximizing values neighbors took is needed to recover a global MAP configuration. Note that computing x∗i is computationally not the same as first computing the marginal distribution of xi and then looking at what value xi maximizes pxi (xi ) because we are actually maximizing over all other xj for j = i as well! Max-Product relates to MAP elimination in the same way Sum-Product relates to elimination. In particular, for trees, we can just root the tree at an arbitrary node and eliminate from the leaves up to the root at which point we can obtain the max-marginal at the root node. We pass messages from the root back to the leaves to obtain max-marginals at all the nodes. 11.1 Max-Product for Undirected Trees As usual, for undirected graph G = (V, E) with V = {1, 2, . . . , N }, we define a distri­ bution over x1 , . . . , xN via factorization: px (x) ∝ φi (xi ) i∈V ψij (xi , xj ). (i,j)∈E We work through an example to illustrate what happens when G is a tree. The in­ tuition is nearly identical to that of Sum-Product. Consider an undirected graph­ ical model shown below. We run the MAP elimination algorithm with ordering (5, 4, 3, 2, 1), i.e., effectively we’re rooting the tree at node 1: x1 x3 x2 x4 x1 x5 x3 x2 x4 m5 Eliminating nodes 5, 4, and 3, we obtain messages: m5 (x2 ) = max φ5 (x5 )ψ25 (x2 , x5 ) x5 ∈X m4 (x2 ) = max φ4 (x4 )ψ24 (x2 , x4 ) x4 ∈X m3 (x1 ) = max φ3 (x3 )ψ13 (x1 , x3 ) x3 ∈X x1 x3 x2 m4 x1 x2 m4 m5 m3 m5 Eliminating node 2 yields message: m2 (x1 ) = max φ2 (x2 )ψ12 (x1 , x2 )m4 (x2 )m5 (x2 ). x2 Finally, we obtain the max-marginal at node 1: px1 (x1 ) = φ1 (x1 )m2 (x1 )m3 (x1 ), with maximal point x∗1 ∈ argmax px1 (x1 ) = argmax φ1 (x1 )m2 (x1 )m3 (x1 ). x1 ∈X x1 ∈X To obtain max-marginals for all nodes, we could imagine that we instead choose a different node to act as the root and eliminate from leaves to the root. As with the Sum-Product case, the key insight is that messages computed previously could be reused. The bottom line is that it suffices to use the same serial schedule as for Sum-Product: pass messages from leaves to the root and then from the root back to the leaves. Using the above serial schedule for message passing, in general the messages are computed as: Y mi→j (xj ) = max φi (xi )ψij (xi , xj ) mk→i (xi ). (1) xi ∈X k∈N (i)\j Once all messages have been computed in both passes, then we compute max-marginals at each node: Y pxi (xi ) = φi (xi ) mk→i (xi ). (2) k∈N (i) 2 As mentioned previously, in general we need to backtrack to obtain a global configu­ ration, which requires us to store the argmax for each message: Y δi→j (xj ) = argmax φi (xi )ψij (xi , xj ) mk→i (xi ). xi ∈X k∈N (i)\j Then the backtracking works as follows: Let x∗i ∈ argmaxxi pxi (xi ) for an arbitrary node i. For each j ∈ N (i), assign x∗j = δj→i (x∗i ). Recurse until all nodes have been assigned a configuration. Parallel Max-Product for Undirected Trees As with Sum-Product, we have a parallel or flood schedule for message passing: 1. Initialize all messages m0i→j (xj ) = 1 for all (i, j) ∈ E. 2. Iteratively apply the update for t = 0, 1, 2, . . . until convergence: Y mt+1 (x ) = max φ (x )ψ (x , x ) mtk→i (xi ). j i i ij i j i→j xi ∈X k∈N (i)\j 3. Compute max-marginals: pxi (xi ) = φi (xi ) Y mt+1 k→i (xi ). k∈N (i) 4. Backtrack to recover an MAP configuration. As with parallel Sum-Product, we can reduce the amount of computation by precom­ puting at the beginning of each iteration of step 2: Y m m ti (xi ) = mtk→i (xi ). k∈N (i) Then the update equation in step 2 becomes mt+1 i→j (xj ) = max φi (xi )ψij (xi , xj ) xi ∈X m m ti (xi ) . t mj→i (xi ) Numerical stability, Max-Sum, and Min-Sum In the message passing and max-marginal computation equations, we are multiplying many potential function values, which can lead to numerical issues. For example, if the potential functions represent probabilities, e.g., if our undirected graph is actually obtained from a DAG so the potentials are conditional probability distributions, then 3 multiplying many probabilities can result in extremely small values below machine precision. A remedy for this is to work in the log domain. Denote mbi→j (xj ) = log mi→j (xj ). Then taking the log of equations (1) and (2) results in updates for the Max-Sum algorithm, summarized below. Message-passing: ⎫ ⎧ ⎬ ⎨ � mbk→i (xi ) . mbi→j (xj ) = max log φi (xi ) + log ψij (xi , xj ) + xi ∈X ⎩ ⎭ k∈N (i)\j Max-marginals: ⎧ ⎫ ⎨ ⎬ � mbk→i (xi ) . pxi (xi ) = exp log φi (xi ) + ⎩ ⎭ k∈N (i) We could also take the negative log instead of the log, which results in the MinSum algorithm. A rationale for why this might make sense is that if the potentials represent probabilities, then log probabilities take on negative values. In this case, if we wish to deal with strictly non-negative entries, then we we would use Min-Sum. Denote bi→j (xj ) = − log mi→j (xj ). Taking the negative log of equations (1) and (2), we obtain the update rules below. Message-passing: ⎧ ⎫ ⎨ ⎬ � bk→i (xi ) . bi→j (xj ) = min − log φi (xi ) − log ψij (xi , xj ) + xi ∈X ⎩ ⎭ k∈N (i)\j Max-marginals: ⎧ ⎫ ⎨ ⎬ � bk→i (xi ) . pxi (xi ) = exp log φi (xi ) − ⎩ ⎭ k∈N (i) Lastly, we remark that if there are ties in MAP configurations, one way to remove the ties is to perturb the potential functions’ values with a little bit of noise. If the model is numerically well-behaved, perturbing the potential functions’ values with a little bit of noise should not drastically change the model but may buy us node max­ marginals that all have unique optima, implying that the global MAP configuration is unique as well. 11.2 Max-Product for Factor Trees The Max-Product algorithm for factor trees is extremely similar to that of undirected trees, so we just present the parallel schedule version of Max-Product on factor tree G = (V, E, F): 4 1. Initialize all messages m0i→a (xi ) = m0a→i (xi ) = 1 for all i ∈ V, a ∈ F. 2. Iteratively apply the update for t = 0, 1, 2, . . . until convergence: Y t mt+1 max fa ({xi } ∪ {xj : j ∈ N (a) \ i}) mj→a (xj ) a→i (xi ) = xj :j∈N (a)\i mt+1 i→a (xi ) = Y j∈N (a)\i mt+1 b→i (xi ) b∈N (i)\a 3. Compute max-marginals: Y pxi (xi ) = mt+1 a→i (xi ). a∈N (i) 4. Backtrack to recover an MAP configuration. 11.3 Max-Product for Hidden Markov Models We revisit our convolution code HMM example from Lecture 9, truncating it to have four hidden states: m1m2 m2m3 m3m4 m4 y1y2 y3y4 y5y6 y7 To ease the notation a bit for what’s to follow, let’s denote mi mi+1 / mi,i+1 and yj yj+1 / yj,j+1 . We thus obtain: m12 m23 m34 m4 y12 y34 y56 y7 Fix ε ∈ (0, 1/2). Let wε = log 1−ε ε > 0. Then the potential functions are given by: ψmi,i+1 ,mi+1,i+2 (“ab”, “cd”) = 1(b = c) for i ∈ {1, 2}, ψm34 ,m4 (“ab”, “c”) = 1(b = c), ψmi,i+1 ,y2i−1,2i (“ab”, “uv”) = exp {1(a = u)wε + 1(a ⊕ b = v)wε } for i ∈ {1, 2, 3}, ψm4 ,y7 (“a”, “u”) = exp {1(a = u)wε } , 5 where a, b, c, d, u, v ∈ {0, 1}, “ab” denotes a length-2 bit-string, and 1(A) is the indi­ cator function that is 1 when event A is true and 0 otherwise. Suppose we observe y = (0, 0, 0, 1, 0, 0, 0) and want to infer the MAP configuration for m = (m1 , m2 , m3 , m4 ). Let’s solve this using Max-Sum. The first thing to do is to incorporate the data into the model by plugging in the observed values; this will cause the resulting graphical model to just be a line: m12 m23 m34 m4 ϕm12 ϕm23 ϕm34 ϕm4 The new singleton potentials are: φm12 (“ab”) = exp {1(a = 0)wε + 1(a ⊕ b = 0}wε } φm23 (“ab”) = exp {1(a = 0)wε + 1(a ⊕ b = 1}wε } φm34 (“ab”) = exp {1(a = 0)wε + 1(a ⊕ b = 0}wε } φm4 (“a”) = exp {1(a = 0)wε } The pairwise potentials remain the same. Since we’re going to use Max-Sum, we write out the log-potentials: log φm12 (“ab”) = 1(a = 0)wε + 1(a ⊕ b = 0}wε log φm23 (“ab”) = 1(a = 0)wε + 1(a ⊕ b = 1}wε log φm34 (“ab”) = 1(a = 0)wε + 1(a ⊕ b = 0}wε log φm4 (“a”) = 1(a = 0)wε log ψmi,i+1 ,mi+1,i+2 (“ab”, “cd”) = 0 if b = c −∞ otherwise log ψm34 ,m4 (“ab”, “c”) = 0 if b = c −∞ otherwise for i ∈ {1, 2} Using the Max-Sum message update equation, log m(i,i+1)→(i+1,i+2) (“cd”) ( = max 2 log φmi,i+1 (“ab”) + log ψmi,i+1 ,mi+1,i+2 (“ab”, “cd”) + log m(i−1,i)→(i,i+1) (“ab”) “ab”∈{0,1} ( = max log φmi,i+1 (“ac”) + log m(i−1,i)→(i,i+1) (“ac”) , a∈{0,1} where the last equality uses the fact that what the pairwise potentials really are saying is that we only allow transitions where the last bit of the previous state is the first 6 bit of the next state. This means that if we know the next state is “cd”, then the previous state must end in “c”, so we shouldn’t bother optimizing over all bit-strings of length 2. Note that if we pass messages from the left-most node to the right-most node, then what log m(i,i+1)→(i+1,i+2) (“cd”) represents is the highest “score” achievable that ends in mi+1,i+2 = “cd” but that does not include the contribution of the observation associated with mi+1,i+2 . The goal of Max-Sum would then be to find the MAP con­ figuration with the highest overall score, with contributions from all nodes’ singleton potentials and pairwise potentials. We can visualize this setup using a trellis diagram: Hidden states m12 m23 m34 m4 00 00 00 0 logϕm12(“00”) = 2wε logϕm23(“00”) = 1wε logϕm34(“00”) = 2wε logϕm4(“0”) = 1wε 01 01 01 logϕm12(“01”) = 1wε logϕm23(“01”) = 2wε logϕm34(“01”) = 1wε 10 10 10 1 logϕm12(“10”) = 0wε logϕm23(“10”) = 1wε logϕm34(“10”) = 0wε logϕm4(“1”) = 0wε 11 11 11 logϕm12(“11”) = 1wε logϕm23(“11”) = 0wε logϕm34(“11”) = 1wε Possible values hidden states can take For each hidden state’s possible value, incoming edges correspond to which possible previous states are possible. The MAP configuration is obtained by summing across scores (i.e., the log singleton potential values in this case) along a path that goes from one of the left-most nodes to the right-most node of the trellis diagram. If we keep track of which edges we take that maximizes the score along the path we traverse, then we can backtrack once we’re at the end to figure out what the best configuration is, as was done in the MAP elimination algorithm from last lecture. The key idea for the Max-Sum algorithm is that we basically iterate through the vertical layers (i.e., hidden states) in the trellis diagram one at a time from left to right. At each vertical layer, we store the highest possible score achievable that arrives at each possible value that the hidden state can take. Note that for this example, the edges do not have scores associated with them; in general, edges may have scores as well. 7 With this intuition, we run the Max-Sum with the modification that we keep track of which optimizing values give rise to the “highest score so far” at each hidden state’s possible value. At the end, we backtrack to obtain the global MAP configuration (00, 00, 00, 0), which achieves a overall score of 6wε . Note that this algorithm of essentially using Max-Product on an HMM while keeping track of optimizing values to enable backtracking at the end is called the Viterbi algorithm, developed in 1967. 8