Forward-Backward Algorithm, Sum-Product on Factor Graphs

Massachusetts Institute of Technology Department of Electrical Engineering and Computer Science 6.438 Algorithms For Inference Fall 2014 9 Forward-backward algorithm, sum-product on factor graphs The previous lecture introduced belief propagation (sum-product), an efficient infer­ ence algorithm for tree structured graphical models. In this lecture, we specialize it further to the so-called hidden Markov model (HMM), a model which is very useful in practice for problems with temporal structure. 9.1 Example: convolution codes We motivate our discussion of HMMs with a kind of code for communication called a convolution code. In general, the problem of communication is that the sender would like to send a message m, represented as a bit string, to a receiver. The message may be corrupted along the way, so we need to introduce redundancy into the message so that it can be reconstructed accurately even in the presence of noise. To do this, the sender sends a coded message b over a noisy channel. The channel introduces some noise (e.g. by flipping random bits). The receiver receives the “received message” y and then applies a decoding procedure to get the decoded message m̂. A schematic is shown in Figure 1. Clearly, we desire a coding scheme where m̂ = m with high probability, b is not much larger than m, and m̂ can be efficiently computed from y . We now discuss one example of a coding scheme, called a convolution code. Sup­ pose the message m consists of N bits. The coded message b will consist of 2N − 1 bits, alternating between the following: • The odd-numbered bits b2i−1 repeat the message bits mi exactly. • The even-numbered bits b2i are the XOR of message bits mi and mi+1 , denoted mi ⊕ mi+1 . The ratio of the lengths of m and b is called the rate of the code, so this convolution code is a rate 12 code, i.e. for every coded message bit, it can convey 12 message bit. We assume an error model called a binary symmetric channel : each of the bits of the coded message is independently flipped with probability ε. We can represent this as a directed graphical model as shown in Figure 2. Note that from the receiver’s perspective, only the yi ’s are observed, and the task is to infer the mi ’s. In order to perform inference, we must convert this graph into an undirected graphical model. Unfortunately, the straightforward construction, where we moralize the graph, does not result in a tree structure, because of the cliques over mi , mi+1 , and b2i . Instead, we coarsen the representation by combining nodes into “supernodes.” In particular, we will combine all of the adjacent message bits into variables mi mi+1 , and 1 error b1 , . . . , b2N m1 , . . . , mN message y1 , . . . , y2N 1 coded message corruption (binary symmetric channel) 1 m̂1 , . . . , m̂N received message decoded message Figure 1: A schematic representation of the problem setup for convolution codes. m1 m3 m2 b1 b2 b3 b4 b5 y1 y2 y3 y4 y5 … Figure 2: Our convolution code can be represented as a directed graphical model. we will combine pairs of adjacent received message bits y2i−1 y2i , as shown in Figure 3. This results in a tree-structured directed graph, and therefore an undirected tree graph — now we can perform sum-product. 9.2 Hidden Markov models Observe that the graph in Figure 3 is Markov in its hidden states. More generally, a hidden Markov model (HMM) is a graphical model with the structure shown in Figure 4. Intuitively, the variables xi represent a state which evolves over time and which we don’t get to observe, so we refer to them as the hidden state. The variables yi are signals which depend on the state at the same time step, and in most applications are observed, so we refer to them as observations. From the definition of directed graphical models, we see that the HMM represents the factorization property P(x1 , . . . , xN , y1 , . . . , yN ) = P(x1 ) N N i=2 P(xi |xi−1 ) N N j=1 P(yj |xj ). (1) Observe that we can convert this to the undirected representation shown in Figure 4 (b) by taking each of the terms in this product to be a potential. This allows us to 2 m1 m2 m2 m3 m3 m4 y 1 y2 y 3 y4 y 5 y6 m1 m2 m2 m3 m3 m4 y 1 y2 y 3 y4 y 5 y6 ... mN y2N 1 mN mN 3 y2N 2 y2N 1 (a) ... mN y2N 1 mN mN 3 y2N 2 y2N 1 (b) Figure 3: (a) The directed graph from figure 2 can be converted to a tree-structured graph with combined variables. (b) The equivalent undirected graph. x1 x2 x3 y1 y2 y3 ... xN Hidden state yN Observations (a) x1 x2 x3 y1 y2 y3 ... xN yN (b) Figure 4: (a) a hidden Markov model, and (b) its undirected equivalent. 3 1 (x1 ) N (x1 ) x1 ... x3 x2 xN 12 (x1 , x2 ) Figure 5: Representation of an HMM as a chain graph. perform inference using sum product on trees. In particular, our goal is typically to infer the marginal distribution for each of the hidden states given all of the observa­ tions. By plugging in the values of the observations, we can convert the HMM to a chain graph, as shown in Figure 5. This graphical model has two sets of potentials: φ1 (x1 ) = P(x1 , y1 ) φi (xi ) = P(yi |xi ), ∀i ∈ {2, . . . , N } ψi,i+1 (xi , xi+1 ) = P(xi+1 |xi ), ∀i ∈ {1, . . . , N − 1} (2a) (2b) (2c) The resulting sum-product messages are: m1→2 (x2 ) = φ1 (x1 )ψ12 (x1 , x2 ) (3) x1 mi→i+1 (xi+1 ) = φi (xi )ψi,i+1 (xi , xi+1 )mi−1→i (xi ) xi mN →N −1 (xN −1 ) = 2≤i≤N φN (xN )ψN −1,N (xN −1 , xN ) (4) (5) xN mi→i−1 (xi−1 ) = φi (xi )ψi−1,i (xi−1 , xi )mi+1→i (xi ) xi 1 ≤ i ≤ N − 1 (6) Belief propagation on HMMs is also known as the forward-backward algorithm. You can easily show that m1→2 (x2 ) = φ1 (x1 )ψ12 (x1 , x2 ) = x1 x1 P(x1 , y1 )P(x2 |x1 ) P(x1 , y1 , x2 ) = P(y1 , x2 ). = x1 m2→3 (x3 ) = φ2 (x2 )ψ23 (x2 , x3 )m1→2 (x2 ) = x2 = x2 P(x3 , x2 , y2 , y1 ) = P(y1 , y2 , x3 ). x2 4 P(y2 |x2 )P(x3 |x2 )P(y1 , x2 ) Continuing in this fashion, we can show that: mi−1→i (xi ) = P(y1 , y2 , . . . , yi−1 , xi ). Similarly: mN →N −1 (xN −1 ) = X φN (xN )ψN −1,N (xN −1 , xN ) = xN = X mN −1→N −2 (xN −2 ) = X xN X xN P(yN |xN )P(xN |xN −1 ) P(xN , yN |xN −1 ) = P(yN |xN −1 ). φN −1 (xN −1 )ψN −2,N −1 (xN −2 , xN −1 ) xN −1 = X xN −1 = X xN −1 P(yN −1 |xN −1 )P(xN −1 |xN −2 )P(yN |xN −1 ) P(yN −1 , yN , xN −1 |xN −2 ) = P(yN −1 , yN |xN −2 ). Continuing in this fashion, we get: mi+1→i (xi ) = P(yi+1 , yi+2 , . . . , yN |xi ). 9.2.1 α, β Forward-Backward Algorithms and Probabilistic interpretation of messages As we have seen, belief propagation on HMMs takes the form of a two-pass algorithm, consisting of a forward pass, and a backward pass. In fact there is considerably flexibility on how computation is structured in the algorithm, even with the two-pass structure. As a result, we refer to this algorithm and its variants collectively as the forward-backward algorithm. To illustrate how the computation can be rearranged in useful ways, in this section we highlight one variant of the forward-backward algorithm, termed the α, β version for reasons that will become apparent. The α, β forward-backward algorithm was also among the earliest versions developed. To see the rearrangement of interest, first note that each posterior marginal of interest pxi |y1 ,...,yN is proportional to pxi ,y1 ,...,yN and that from the graphical model in Figure 4, (y1 , . . . , yi ) are d-separated from (yi+1 , . . . , yN ) given xi . Therefore, the joint distribution factorizes as: py1 ,...,yN ,xi = py1 ,...,yi ,xi pyi+1 ,...,yN |xi ,y1 ,...,yi = py1 ,...,yi ,xi pyi+1 ,...,yN |xi . 5 (7) (8) Now we derive recursive update rules for computing each of the two parts. First, the forward messages α: αi (xi ) = P(y1 , . . . , yi , xi ) X = P(y1 , . . . , yi , xi , xi−1 ) (9) (10) xi−1 X = xi−1 X = xi−1 P(y1 , . . . , yi−1 , xi−1 )P(xi |xi−1 )P(yi |xi ) (11) αi−1 (xi−1 )P(xi |xi−1 )P(yi |xi ), ∀i ∈ {2, . . . , N }, (12) where α1 (x1 ) = P(x1 , y1 ). Then, the backward messages β: βi (xi ) = P(yi+1 , . . . , yN |xi ) X = P(xi+1 , yi+1 , . . . , yN |xi ) (13) (14) xi+1 = X xi+1 = X xi+1 P(yi+1 |xi+1 )P(xi+1 |xi )P(yi+2 , . . . , yN |xi+1 ) (15) P(yi+1 |xi+1 )P(xi+1 |xi )βi+1 (xi+1 ), ∀i ∈ {1, . . . , N − 1}, (16) where βN (xN ) = 1. Now consider the relationship between α and β and our sum-product messages. The forward message corresponds to part of the formula for α, X mi→i+1 (xi+1 ) = P(y1 , . . . , yi , xi+1 ) = P(y1 , . . . , yi , xi , xi+1 ) xi = X xi P(xi+1 |xi )P(y1 , . . . , yi , xi ) = X xi P(xi+1 |xi )αi (xi ). (17) Observe that the α messages and the forward messages in belief propagation perform exactly the same series of sums and products, but they divide up the steps differently. Following an analogous line of reasoning, we find that the backwards messages are identical, i.e. βi (xi ) = mi+1→i (xi ). This gives a useful probabilistic interpretation to the sum-product messages. We can see directly that the marginal computation corresponding to γi (xi ) / pxi (xi |y1 , . . . , yN ) ∝ pxi (xi , y1 , . . . , yN ) = φi (xi )mi−1→i (xi )mi+1→i (xi ) in the sum-product algorithm is, in terms of the and α, β messages, γi (xi ) = αi (xi )βi (xi ) . P(y1 , . . . , yN ) 6 x1 fa fb x3 x2 fc Figure 6: An example of factor graph belief propagation. In turn, we note that the likelihood P(y1 , . . . , yN ) is obtained from the messages at any node: X P(y1 , . . . , yN ) = αi (xi )βi (xi ), ∀i. xi One appealing feature of this form of the forward-backward algorithm is that the α messages are, themselves, marginals that are often directly useful in applications. In particular, αi (xi ) ∝ P(xi |y1 , . . . , yi ), i.e., αi (xi ) represents the marginal at state node i using data up to observation node i, and as such represents a “causal” marginal of particular interest in real-time appli­ cations. Sometimes these are referred to as “filtered” marginals, for reasons that we will come back to when we revisit the forward-backward algorithm in the context of Gaussian HMMs. Finally, it should be noted that several other closely related variants of the forwardbackward algorithm are also possible. One is the so-called α, γ version, whereby the recursion for β is replaced with a recursion for γ, which has the benefit of producing the desired marginals directly as messages, and the additional benefit that each piece of data yi is used only once, in the forward pass. Specifically, the recursion for the γ messages can be shown to be (see, e.g., Jordan’s notes, Chapter 12): � � X P(xi+1 |xi )αi (xi ) � γi+1 (xi+1 ). γi (xi ) = � � � P(xi+1 |xi )αi (xi ) x x i+1 9.3 i Sum-product on factor graphs We now consider how to perform the sum product algorithm on a slightly more general class of graphical models, tree-structured factor graphs. Observe that this is a strictly larger set of graphs than undirected or directed trees. In either of these two cases, there is one factor for each edge in the original graph, so their equivalent factor graph representation is still a tree. However, some non-tree-structured directed or undirected graphs may have tree structured factor graph representations. One 7 important example is polytrees: recall that these are directed graphs which are treestructured when we ignore the directions of the edges. For instance, observe that the convolution code directed graph shown in Figure 2 is a polytree, even though its undirected representation (obtained by moralizing) is not a tree. Using factor graph belief propagation, it is possible to perform inference in this graph without resorting to the supernode representation of Section 9.1. In factor graph belief propagation, messages are sent between variable nodes and factor nodes. As in undirected belief propagation, each node sends messages to one neighbor by multiplying and/or summing messages from its other neighbors. Factor nodes a multiply incoming messages by their factor and sum out all but the relevant variable i: N X ma→i (xi ) = fa (xi , xN (a)\{i} ) mj→a (xj ) (18) xN (a)\{i} j∈N (a)\{i} The variable nodes simply multiply together their incoming messages: N mb→i (xi ) mi→a (xi ) = (19) b∈N (i)\{a} For instance, consider Figure 6. The messages required to compute px1 (x1 ) are: ma→1 (x1 ) m2→b (x2 ) mc→3 (x3 ) m3→b (x3 ) = = = = fa (x1 ) 1 fc (x3 ) mc→3 (x3 ) X mb→1 (x1 ) = fb (x1 , x2 , x3 )m2→b (x2 )m3→b m(x3 ) x2 ,x3 8 (20) (21) (22) (23) (24)