The Sum-Product Algorithm

Chapter 12 The sum-product algorithm The sum-product algorithm is the basic “decoding” algorithm for codes on graphs. For finite cycle-free graphs, it is finite and exact. However, because all its operations are local, it may also be applied to graphs with cycles; then it becomes iterative and approximate, but in coding applications it often works very well. It has become the standard decoding algorithm for capacity-approaching codes (e.g., turbo codes, LDPC codes). There are many variants and applications of the sum-product algorithm. The most straightforward application is to a posteriori probability (APP) decoding. When applied to a trellis, it becomes the celebrated BCJR decoding algorithm. In the field of statistical inference, it becomes the even more widely known “belief propagation” (BP) algorithm. For Gaussian statespace models, it becomes the Kalman smoother. There is also a “min-sum” or maximum-likelihood sequence detection (MLSD) version of the sum-product algorithm. When applied to a trellis, the min-sum algorithm gives the same result as the Viterbi algorithm. 12.1 The sum-product algorithm on cycle-free graphs We will develop the sum-product algorithm as an APP decoding algorithm for a code C that has a cycle-free normal graph realization. We then discuss generalizations. The code C is therefore described by a realization involving a certain set of symbol variables {Yi , i ∈ I} represented by half-edges (dongles), a certain set of state variables {Σj , j ∈ J } represented by edges, and a certain set of constraint codes {Ck , k ∈ K} of arbitrary degree, such that the graph of the realization is cycle-free; i.e., every edge (and obviously every half-edge) is by itself a cut set. APP decoding is defined in general as follows. We assume that a set of independent observations are made on all symbol variables {Yi , i ∈ I}, resulting in a set of observations r = {ri , i ∈ I} and likelihood vectors {{p(ri | yi ), yi ∈ Yi }, i ∈ I}, where Yi is the alphabet of Yi . The likelihood
of a codeword y = {yi , i ∈ I} ∈ C is then defined as the componentwise product p(r | y) = i∈I p(ri | yi ). Assuming equiprobable codewords, the a posteriori probabilities {p(y | r), y ∈ C} (APPs) are proportional to the likelihoods {p(r | y), y ∈ C}, since by Bayes’ law, p(y | r) = p(r | y)p(y) ∝ p(r | y), p(r) y ∈ C. Let Ci (yi ) denote the subset of codewords in which the symbol variable Yi has the value yi ∈ Yi . Then, up to a scale factor, the symbol APP vector {p(Yi = yi | r), yi ∈ Yi } is given by     p(Yi = yi | r) = p(y | r) ∝ p(r | y) = p(ri
| yi
), yi ∈ Yi . (12.1) y∈Ci (yi ) y∈Ci (yi ) i
∈I y∈Ci (yi ) Similarly, if Cj (sj ) denotes the subset of codewords that are consistent with the state variable Σj having the value sj in the state alphabet Sj , then, up to a scale factor, the state APP vector {p(Σj = sj | r), sj ∈ Sj } is given by   p(Σj = sj | r) ∝ p(ri | yi ), sj ∈ Sj . (12.2) y∈Cj (sj ) i∈I We see that the components of APP vectors are naturally expressed as sums of products. The sum-product algorithm aims to compute these APP vectors for every state and symbol variable. 12.1.1 Past-future decomposition rule The sum-product algorithm is based on two fundamental principles, which we shall call here the past/future decomposition rule and the sum-product update rule. Both of these rules are based on set-theoretic decompositions that are derived from the code graph. The past/future decomposition rule is based on the Cartesian-product decomposition of the cut-set bound (Chapter 11). In this case every edge Σj is a cut set, so the subset of codewords that are consistent with the state variable Σj having the value sj is the Cartesian product Cj (sj ) = Y|P (sj ) × Y|F (sj ), (12.3) where P and F denote the two components of the disconnected graph which results from deleting the edge representing Σj , and Y|P (sj ) and Y|F (sj ) are the sets of symbol values in each component that are consistent with Σj taking the value sj . We now apply an elementary Cartesian-product lemma: Lemma 12.1 (Cartesian-product distributive law) If X and Y are disjoint discrete sets and f (x) and g(y) are any two functions defined on X and Y, then ⎞   ⎛    f (x) ⎝ g(y)⎠ . (12.4) f (x)g(y) = (x,y)∈X ×Y x∈X y∈Y This lemma may be proved simply by writing the terms on the right in a rectangular array and then identifying them with the terms on the left. It says that rather than computing the sum of |X ||Y| products, we can just compute a single product of independent sums over X and Y. This simple lemma lies at the heart of many “fast” algorithms. Using (12.3) and applying this lemma in (12.2), we obtain the past/future decomposition rule ⎛ ⎞⎛ ⎞     p(ri | yi )⎠ ⎝ p(ri | yi )⎠ p(Σj = sj | r) ∝ ⎝ y|P ∈Y|P (sj ) i∈IP y|F ∈Y|F (sj ) i∈IF ∝ p(Σj = sj | r|P )p(Σj = sj | r|F ), (12.5) in which the two terms p(Σj = sj | r|P ) and p(Σj = sj | r|F ) depend only on the likelihoods of the past symbols y|P = {Yi , i ∈ IP } and future symbols y|F = {Yi , i ∈ IF }, respectively. The sum-product algorithm therefore computes the APP vectors {p(Σj = sj | r|P )} and {p(Σj = sj | r|F )} separately, and multiplies them componentwise to obtain {p(Σj = sj | r)}. This is the past/future decomposition rule for state variables. APP vectors for symbol variables are computed similarly. In this case, since symbol variables have degree 1, one of the two components of the graph induced by a cut is just the symbol variable itself, while the other component is the rest of the graph. The past/future decomposition rule thus reduces to the following simple factorization of (12.1): ⎛ ⎞   p(ri
| yi
)⎠ ∝ p(yi | ri )p(yi | r|i
= (12.6) p(Yi = yi | r) ∝ p(ri | yi ) ⎝  i ). y∈Ci (yi ) i
=i In the turbo code literature, the first term, p(yi | ri ), is called the intrinsic information, while the second term, p(yi | r|i
=i ), is called the extrinsic information. 12.1.2 Sum-product update rule The second fundamental principle of the sum-product algorithm is the sum-product update rule. This is a local rule for the calculation of an APP vector, e.g., {p(Σj = sj | r|P ), sj ∈ Sj }, from APP vectors that lie one step further upstream. The local configuration with respect to the edge corresponding to the state variable Σj is illustrated in Figure 1. The edge must be incident on a unique past vertex corresponding to a constraint code Ck . If the degree of Ck is δk , then there are δk − 1 edges further upstream of Ck , corresponding to further past state or symbol variables. For simplicity, we suppose that these are all state variables {Σj
, j  ∈ Kjk }, where we denote their index set by Kjk ⊆ K|P . {Σj
, j  ∈ Kjk } Y|Pj
 @ @ H H@ H ...     Ck Σj Figure 1. Local configuration for sum-product update rule. Since the graph is cycle-free, each of these past edges has its own independent past Pj
. The corresponding sets Y|Pj
of input symbols must be disjoint, and their union must be Y|P . Thus if Ck (sj ) is the set of codewords in the local constraint code Ck that are consistent with Σj = sj , and Y|Pj
(sj
) is the set of y|Pj
∈ Y|Pj
that are consistent with Σj
= sj
, then we have Y|P (sj ) = Ck (sj ) j
∈Kjk Y|Pj
(sj
), (12.7) where the plus sign indicates a disjoint union, and the product sign indicates a Cartesian product. In other words, for each codeword in Ck for which Σj = sj , the set of possible pasts is the Cartesian product of possible pasts of the other state values {sj
, j  ∈ Kjk }, and the total set of possible pasts is the disjoint union of these Cartesian products. Now, again using the Cartesian-product distributive law, it follows from (12.7) that   p(Σj
= sj
| r|Pj
). p(Σj = sj | r|P ) = Ck (sj ) j
∈K (12.8) jk Thus if we know all the upstream APP vectors {p(Σj
= sj
| r|Pj
), sj
∈ Sj
}, then we can compute the APP vector {p(Σj = sj | r|P ), sj ∈ Sj }. Equation (12.8) is the sum-product update rule. We can see that for each sj ∈ Sj it involves a sum of |Ck | products of δk − 1 terms. Its complexity is thus proportional to the size |Ck | of the constraint code Ck . In a trellis, this is what we call the branch complexity. In the special case where Ck is a repetition code, there is only one codeword corresponding to each sj ∈ Sj , so (12.8) becomes simply the following product update rule:  p(Σj = sj | r|P ) = p(Σj
= sj | r|Pj
); (12.9) j
∈Kjk i.e., the components of the upstream APP vectors are simply multiplied componentwise. When the sum-product algorithm is described for Tanner graphs, the product update rule is often stated as a separate rule for variable nodes, because variable nodes in Tanner graphs correspond to repetition codes in normal graphs. Note that for a repetition code of degree 2, the product update rule of (12.9) simply becomes a pass-through of the APP vector; no computation is required. This seems like a good reason to suppress state nodes of degree 2, as we do in normal graphs. 12.1.3 The sum-product algorithm Now we describe the complete sum-product algorithm for a finite cycle-free normal graph, using the past/future decomposition rule (12.5) and the sum-product update rule (12.8). Because the graph is cycle-free, it is a tree. Symbol variables have degree 1 and correspond to leaves of the tree. State variables have degree 2 and correspond to branches. For each edge, we wish to compute two APP vectors, corresponding to past and future. These two vectors can be thought of as two messages going in opposite directions. Using the sum-product update rule, each message may be computed after all upstream messages have been received at the upstream vertex (see Figure 1). Therefore we can think of each vertex as a processor that computes an outgoing message on each edge after it has received incoming messages on all other edges. Because each edge is the root of a finite past tree and a finite future tree, there is a maximum number d of edges to get from any given edge to the furthest leaf node in either direction, which is called its depth d. If a message has depth d, then the depth of any upstream message can be no greater than d − 1. All symbol half-edges have depth d = 0, and all state edges have depth d ≥ 1. The diameter dmax of the tree is the maximum depth of any message. Initially, incoming messages (intrinsic information) are available at all leaves of the tree. All depth-1 messages can then be computed from these depth-0 messages; all depth-2 messages can then be computed from depth-1 and depth-0 messages; etc. In a synchronous (clocked) system, all messages can therefore be computed in dmax clock cycles. Finally, given the two messages on each edge in both directions, all a posteriori probabilities (APPs) can be computed using the past/future decomposition rule (12.5). In summary, given a finite cycle-free normal graph of diameter dmax and intrinsic information for each symbol variable, the sum-product algorithm computes the APPs of all symbol and state variables in dmax clock cycles. One message (APP vector) of size |Sj | is computed for each state variable Σj in each direction. The computational complexity at a vertex corresponding to a constraint code Ck is of the order of |Ck |. (More precisely, the number of pairwise multiplications required is δk (δk − 2)|Ck |.) The sum-product algorithm does not actually require a clock. In an asynchronous implementation, each vertex processor can continuously generate outgoing messages on all incident edges, using whatever incoming messages are available. Eventually all messages must be correct. An analog asynchronous implementation can be extremely fast. We see that there is a clean separation of functions when the sum-product algorithm is implemented on a normal graph. All computations take place at vertices, and the computational complexity at a vertex is proportional to the vertex (constraint code) complexity. The function of ordinary edges (state variables) is purely message-passing (communications), and the communications complexity (bandwidth) is proportional to the edge complexity (state space size). The function of half-edges (symbol variables) is purely input/output; the inputs are the intrinsic APPs, and the ultimate outputs are the extrinsic APP vectors, which combine with the inputs to form the symbol APPs. In integrated-circuit terminology, the constraint codes, state variables and symbol variables correspond to logic, interconnect, and I/O, respectively. 12.2 The BCJR algorithm The chain graph of a trellis (state-space) representation is the archetype of a cycle-free graph. The sum-product algorithm therefore may be used for exact APP decoding on any trellis (statespace) graph. In coding, the resulting algorithm is known as the Bahl-Cocke-Jelinek-Raviv (BCJR) algorithm. (In statistical inference, it is known as the forward-backward algorithm. If all probability distributions are Gaussian, then it becomes the Kalman smoother.) Figure 2 shows the flow of messages and computations when the sum-product algorithm is applied to a trellis. Y0 ι0 ?6ε0 C0 α1  β1 Y1 ι1 ?6ε1 C1 α2  β2 Y2 ι2 ?6ε2 C2 α3  β3 Y3 ι3 ?6ε3 C3 α4  β4 Y4 ι4 ?6ε4 C4 α5  β5 Y5 ι5 ?6ε5 C5 α6  β6 Y6 ι6 ?6ε6 C6 α7  Y7 ι7 ?6ε7 C7 β7 Figure 2. Flow of messages and computations in the sum-product algorithm on a trellis. The input messages are the intrinsic APP vectors ιi = {p(ri | yi ), yi ∈ Yi }, derived from the observations ri ; the output messages are the extrinsic APP vectors εi = {p(yi | r|i
=i ), yi ∈ Yi }. The intermediate messages are the forward state APP vectors αj = p(sj | r|Pj ), sj ∈ Sj } and the backward state APP vectors βj = p(sj | r|Fj ), sj ∈ Sj }, where r|Pj and r|Fj denote the observations before and after sj , respectively. The algorithm proceeds independently in the forward and backward directions. In the forward direction, the messages αj are computed from left to right; αj may be computed by the sumproduct rule from the previous message αj−1 and the most recent input message ιj−1 . In the backward direction, the messages βj are computed from right to left; βj may be computed by the sum-product rule from βj+1 and ιj . Finally, each output message εi may be computed by the sum-product update rule from the messages αi and βi+1 , giving the extrinsic information for each symbol. To find the APP vector of an input symbol Yi , the intrinsic and extrinsic messages ιi and εi are multiplied componentwise, according to the past/future decomposition rule. (In turbo decoding, the desired output is actually the extrinsic likelihood vector, not the APP vector.) Similarly, to find the APP vector of a state variable Σj , the forward and backward messages αj and βj are multiplied componentwise. Exercise 1. Consider the two-state trellis diagram for the binary (7, 6, 2) SPC code shown in Figure 3 of Chapter 10. Suppose that a codeword is chosen equiprobably at random, that the transmitter maps {0, 1} to {±1} as usual, that the resulting real numbers are sent through a discrete-time AWGN channel with noise variance σ 2 = 1 per symbol, and that the received sequence is r = (0.4, −1.0, −0.1, 0.6, 0.7, −0.5, 0.2). Use the sum-product algorithm to determine the APP that each input bit Yi is a 0 or a 1. 12.3 The min-sum algorithm and ML decoding We now show that with minor modifications the sum-product algorithm may be used to perform a variant of maximum-likelihood (ML) sequence decoding rather than APP decoding. On a trellis, the resulting “min-sum” algorithm becomes a variant of the Viterbi algorithm. With the same notation as in the previous section, the min-sum algorithm is defined as follows. Again, let Ci (yi ) denote the subset of codewords in which the symbol variable Yi has the value yi ∈ Yi . Then the metric mi (yi ) of yi is defined as the maximum likelihood of any codeword y ∈ Ci (Yi ); i.e., mi (yi ) = max p(r | y) = max y∈Ci (Yi ) y∈Ci (Yi )  p(ri
| yi
), yi ∈ Yi . (12.10) i
∈I It is clear that the symbol value yi with the maximum metric mi (yi ) will be the value of yi in the codeword y ∈ C that has the maximum global likelihood. Similarly, if Cj (sj ) denotes the subset of codewords that are consistent with the state variable Σj having the value sj in the state alphabet Sj , then the metric mj (sj ) of sj will be defined as the maximum likelihood of any codeword y ∈ Cj (sj ):  max p(ri | yi ), sj ∈ Sj . (12.11) mj (sj ) = y∈Cj (sj ) i∈I We recognize that (12.10) and (12.11) are almost identical to (12.1) and (12.2), with the exception that the sum operator is replaced by a max operator. This suggests that these metrics could be computed by a version of the sum-product algorithm in which “sum” is replaced by “max” everywhere, giving what is called the “max-product algorithm.” In fact this works. The reason is that the operators “max” and “product” operate on probability vectors defined on sets according to the same rules as “sum” and “product.” In particular, assuming that all quantities are non-negative, we have (a) the distributive law: a max{b, c} = max{ab, ac}; (b) the Cartesian-product distributive law: max (x,y)∈X ×Y f (x)g(y) = max f (x) x∈X max g(y) . y∈Y (12.12) Consequently, the derivation of the previous section goes through with just this one change. From (12.3), we now obtain the past/future decomposition rule ⎛ ⎞⎛ ⎞   p(ri | yi )⎠ ⎝ max mj (sj ) = ⎝ max p(ri | yi )⎠ = mj (sj | r|P )mj (sj | r|F ), y|P ∈Y|P (sj ) i∈IP y|F ∈Y|F (sj ) i∈IF in which the partial metrics mj (sj | r|P ) and mj (sj | r|F ) are the maximum likelihoods over the past symbols y|P = {yi , i ∈ IP } and future symbols y|F = {yi , i ∈ IF }, respectively. Similarly, we obtain the max-product update rule  mj (sj | r|P ) = max mj
(sj
| r|Pj
), (12.13) Ck (sj ) j
∈Kjk where the notation is as in the sum-product update rule (12.8). In practice, likelihoods are usually converted to log likelihoods, which converts products to sums and yields the max-sum algorithm. Or, log likelihoods may be converted to negative log likelihoods, which converts max to min and yields the min-sum algorithm. These variations are all trivially equivalent. On a trellis, the forward part of any of these algorithms is equivalent to the Viterbi algorithm (VA). The update rule (12.13) becomes the add-compare-select operation, which is carried out at each state to determine the new metric mj (sj | r|P ) of each state. The VA avoids the backward part of the algorithm by also remembering the survivor history at each state, and then doing a traceback when it gets to the end of the trellis; this traceback corresponds in some sense to the backward part of the sum-product algorithm. Exercise 2. Repeat Exercise 1, using the min-sum algorithm instead of the sum-product algorithm. Decode the same sequence using the Viterbi algorithm, and show how the two computations correspond. Decode the same sequence using Wagner decoding, and show how Wagner decoding relates to the other two methods. 12.4 The sum-product algorithm on graphs with cycles On a graph with cycles, there are several basic approaches to decoding. One approach is to agglomerate the graph enough to eliminate the cycles, and then apply the sum-product algorithm, which will now be exact. The problem is that the complexity of decoding of a cycle-free graph of C cannot be significantly less than the complexity of decoding some trellis for C, as we saw in Chapter 11. Moreover, as we saw in Chapter 10, the complexity of a minimal trellis for a sequence of codes with positive rates and coding gains must increase exponentially with code length. A second approach is simply to apply the sum-product algorithm to the graph with cycles and hope for the best. Because the sum-product rule is local, it may be implemented at any vertex of the graph, using whatever incoming messages are currently available. In a parallel or “flooding” schedule, the sum-product rule is computed at each vertex at all possible times, converting the incoming messages to a set of outgoing messages on all edges. Other schedules are possible, as we will discuss in the next chapter. There is now no guarantee that the sum-product algorithm will converge. In practice, the sum-product algorithm converges with probability near 1 when the code rate is below some threshold which is below but near the Shannon limit. Convergence is slow when the code rate is near the threshold, but rapid when the code rate is somewhat lower. The identification of fixed points of the sum-product algorithm is a topic of current research. Even if the sum-product algorithm converges, there is no guarantee that it will converge to the correct likelihoods or APPs. In general, the converged APPs will be too optimistic (overconfident), because they assume that all messages are from independent inputs, whereas in fact messages enter repeatedly into sum-product updates because of graph cycles. Consequently, decoding performance is suboptimal. In general, the suboptimality is great when the graph has many short cycles, and becomes negligible as cycles get long and sparse (the graph becomes “locally tree-like”). This is why belief propagation has long been considered to be inapplicable to most graphical models with cycles, which typically are based on physical models with inherently short cycles; in coding, by contrast, cycles can be designed to be very long with high probability. A third approach is to beef up the sum-product algorithm so that it still performs well on certain classes of graphs with cycles. Because the sum-product algorithm already works so well in coding applications, this approach is not really needed for coding. However, this is a current topic of research for more general applications in artificial intelligence, optimization and physics.