IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 47, NO.

2, FRBRUARY 2001 549

Dynamic Programming and the Graphical

Representation of Error-Correcting Codes
Stuart Geman, Member, IEEE, and Kevin Kochanek

Abstract—Graphical representations of codes facilitate the We will concentrate here on two particular kinds of compu-
design of computationally efficient decoding algorithms. This is tations: computation of the most likely configuration of a col-
an example of a general connection between dependency graphs, lection of random variables, and this configuration’s associated
as arise in the representations of Markov random fields, and the
dynamic programming principle. We concentrate on two com- probability. In the context of maximum-likelihood (“soft”) de-
putational tasks: finding the maximum-likelihood codeword and coding, we calculate these functionals under the posterior dis-
finding its posterior probability, given a signal received through tribution—the conditional distribution on the transmitted code-
a noisy channel. These two computations lend themselves to a words given the received signal. What codeword was most likely
particularly elegant version of dynamic programming, whereby transmitted? What is the conditional probability that this max-
the decoding complexity is particularly transparent. We explore
some codes and some graphical representations designed specifi- imum-likelihood codeword was in fact transmitted? In other
cally to facilitate computation. We further explore a coarse-to-fine words, we will compute the maximum a posteriori (MAP) esti-
version of dynamic programming that can produce an exact mator of the transmitted codeword and the (conditional) proba-
maximum-likelihood decoding many orders of magnitude faster bility that the estimator is actually correct. MAP decoding gives
than ordinary dynamic programming. the most likely codeword, and therefore minimizes the prob-
Index Terms—Dynamic programming, graphical models, max- ability of a decoding error. If, alternatively, we seek to mini-
imum-likelihood decoding, soft decoding. mize the number of information-bit errors, then we would max-
imize the posterior marginal at each information bit. In general,
I. INTRODUCTION these are not the same thing, and the better choice is clearly
problem-dependent. MAP makes sense when the codeword rep-

B Y now, nearly everyone in the coding community knows

about the tight relationship between the dependency struc-
ture on a set of random variables and the costs of computing
resents a logical unit, perhaps a product name, a price, a cate-
gory, or an English word. On the other hand, a stream of infor-
mation bits with no a priori structure calls for minimizing the
certain functionals on their joint distribution. One way to de- bit-error rate.
duce this relationship is to generalize, more-or-less straightfor- This paper is about using graphical representations of codes
wardly, the principle of dynamic programming (Bellman [1]), and probability models of channels to calculate the exact MAP
but it can be arrived at from many other directions as well. In decoding and its probability. The emphasis is on finding rep-
fact, an amazing history of discovery and rediscovery emerges resentations and algorithms that make these calculations com-
from a diverse range of applications. Indeed, in one way or putationally feasible. It should perhaps be pointed out that, al-
another, the “forward–backward” algorithm [2], “peeling” [3], though the connections between graphs, codes, and computa-
“parsing” [4], the Viterbi algorithm [5], the “sum–product” al- tion are currently the subjects of intense research, most of the
gorithm [6]–[9], “bucket elimination” [10], and “evidence prop- effort is in a somewhat different direction. The emphasis, in-
agation” [11]–[13], all rest on this very same dynamic-program- stead, is on iterative decoding algorithms that are generally not
ming principle. Oftentimes there is a good reason for using exact and whose convergence properties are generally not well
one variation over another. With one approach, there may be understood, but often exhibit, nevertheless, spectacular perfor-
a book-keeping scheme that avoids duplicating computations; mance. Many codes which would be impossible to decode ex-
with another, the computational flow may be particularly trans- actly lend themselves effectively and efficiently to iterative de-
parent or convenient. The best choice is typically dictated by the coding, in some cases even approaching channel capacity. Gal-
computational constraints and the functionals of interest. lager’s [14] low-density parity-check codes and associated iter-
ative decoding algorithm are probably the first examples, but
few if any recognized the power of Gallager’s construction.
Manuscript received January 1, 2000; revised September 5, 2000. Tanner [15] found clean graphical representations of Gallager
This work was supported in part by the Army Research Office under
Contract DAAH04-96-1-0445, the National Science Foundation under and other related codes, and utilized these representations to
Grant DMS-9217655, and the Office of Naval Research under Contract address a host of design issues concerning computational com-
N00014-97-0249. plexity and minimum distance. Tanner’s graphical representa-
S. Geman is with the Division of Applied Mathematics, Brown University,
Providence, RI 02912 USA (e-mail: geman@[Link]). tions helped to set the stage for the rediscovery of Gallager
K. Kochanek was with the Division of Applied Mathematics, Brown Univer- codes, and for the many extensions that have evolved over the
sity, Providence, RI 02912 USA. He is now with Mathematical Technologies past few years. (See also Bahl et al. [16] for an early paper an-
Inc., Providence, RI 02906 USA (e-mail: kochanek@[Link]).
Communicated by B. J. Frey, Guest Editor. ticipating, rather remarkably, the modern viewpoint.) Of course,
Publisher Item Identifier S 0018-9448(01)00727-1. the turbo codes of Berrou et al. [17] also spurred this line of re-
0018–9448/01$10.00 © 2001 IEEE
550 IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 47, NO. 2, FRBRUARY 2001

search, since these codes, as well, have natural graphical rep- V, we introduce two (not necessarily exclusive) methods for
resentations, come equipped with an iterative decoding algo- reducing the computational demands of maximum-likelihood
rithm, and appear to perform near optimally in certain regimes. decoding. The first is a “thinning” algorithm which controls
Wiberg and collaborators (see [6] and [7]) picked up on Tanner computational costs by reducing information density. The
graphs and turbo codes, and within a general treatment of codes second, coarse-to-fine dynamic programming [22], is a kind
on graphs made connections to soft decoding, general channel of multiscale version of dynamic programming that produces
models (with memory), and iterative and noniterative decoding a provably optimal configuration, often with greatly reduced
algorithms. MacKay and Neal (see MacKay [18] for a thorough computation. We present an exact coarse-to-fine algorithm
discussion) developed some of these same themes, and, addi- for some “context-free” codes, including the RM codes,
tionally, introduced related codes and decoding methods that ap- and demonstrate thousandsfold improvements in decoding
pear to be among the current best performers. The connections efficiency.
and common themes among all of these approaches, as well as
to a broader framework including Bayesian inference, Markov II. DEPENDENCY GRAPHS AND COMPUTING
random fields, belief propagation, and the modern theory of
Given a collection of random variables and a
expert systems, seem to have been first understood and most
probability distribution
clearly formulated by Kschischang and Frey [19].
Our goal in this paper is twofold. We first present a brief tu-
torial on Markov random fields (MRFs), dependency graphs,
and the noniterative computation of various functionals on mar- how difficult is it to compute things like the marginal distri-
ginal and posterior distributions. In the remainder of the paper, bution on , or the most likely configuration ?
we present some new results on how to apply these computa- More than anything else, the dependency relationships among
tional strategies to facilitate the maximum-likelihood decoding the random variables dictate the complexity of computing
of graphical codes. As we have said, our focus is somewhat dif- these and other functionals of . Dependency relationships
ferent from the current trend, in that we concentrate on repre- can be conveniently represented with a graphical structure,
sentations that allow exact computation, and on new methods from which the complexity of various computations can be
for reducing computational complexity. But in light of the often more or less “read off.” These graphical representations, and
good performance of iterative methods, it would be of great in- their connection to computing, are the foundation of modern
terest to systematically compare iterative and exact computa- expert systems [11]–[13] as well as of speech-recognition
tions in a decoding problem that lends itself to both approaches. technologies [2], [23]. In fact, the connection is quite general,
We have not yet made these comparisons. having emerged and re-emerged in these and many other
Section II is about generic computational issues on graphs. application areas, such as genetics [3], coding theory [5], [17],
There is an especially transparent connection between graphs computational linguistics [4], and image analysis [24].
that represent dependency structures and the efficient calcula- This connection between graphs and computing is funda-
tion of a most likely configuration and its associated probability. mental to the modern treatment of soft decoding, and we begin
This is essentially the “bucket elimination” algorithm, and it is here with a brief tutorial. (See also Frey [25] for an excellent
perhaps the most straightforward of the various generalizations introduction to some of the same material.) Our approach to
of dynamic programming. There is no need to triangulate, no graphs is through Markov random fields, and our computational
need to construct junction trees, and no need to worry about cy- focus is on computing most likely configurations and their
cles. Conveniently, the computational cost is readily calculated corresponding probabilities.
once a site-visitation schedule has been established. Further-
more, it is often the case that the most efficient site-visitation A. Markov Random Fields and Their Gibbs Representations
schedule is immediately apparent from the general structure of
the dependency graph. For the purpose of constructing dependency graphs, it is con-
In Section III, we consider the simplest possible case: venient to index random variables by a general finite index set
graphs with linear structure. We review their connection to , rather than the more traditional set . If state
convolutional codes and revisit Viterbi decoding from the MRF spaces are finite, then is a finite vector of random
viewpoint. We point out that the (posterior) probability that the variables with finite range. In most applications, different com-
Viterbi-decoded signal is correct can be calculated as easily ponents have different state spaces, but merely for the sake of
as the decoding itself, and we discuss extensions to channels cleaner notation, we will assume a common (and finite) range
with Markov memory, designed to model burst noise [20], . Thus, is a probability on .
[7]. In Section IV, we generalize to the case of tree-structured is a Markov random field (MRF) with respect to a graph
graphs and make a connection to production systems and if is strictly positive
context-free grammars. We reformulate Forney’s squaring and if
construction [21] to obtain a grammatical representation of
(1)
the Reed–Muller (RM) and other related codes. Moreover,
we discuss the computational complexity of soft decoding
for all and all , where
or of evaluating posterior probabilities for both memoryless
and Markov communications channels. Finally, in Section • indexes the nodes of ;
GEMAN AND KOCHANEK: PROGRAMMING AND GRAPHICAL REPRESENTATION OF ERROR-CORRECTING CODES 551

(turn arrows into edges and then connect all parents of each
daughter).1
If our only interest is in graphical representations of Gibbs
distributions, then strict positivity can be relaxed. In particular,
even if we drop the condition in the definition of
Gibbs distributions, we still have the Markov property (1) wrt ,
provided that we avoid conditioning on events with probability
Fig. 1. Dependency graph for a simple Markov random field.
zero. The computational analysis of graphical models, as it turns
out, relies only on the implication Gibbs MRF. Therefore, we
shall proceed without the positivity constraint , which
• is the neighborhood structure of , mean- would otherwise be troublesome in some of our applications.
ing that if and only if and there is an edge
connecting and in ; and B. Marginal and Posterior Distributions
• is shorthand for . In just about every application of MRFs, including coding,
The distribution on the random variable associated with any site we are interested in making inferences about a subset of
, given the values of all other random variables, depends the variables given observations of the remaining variables.
only on the values of the neighbors. This generalizes the familiar What makes MRFs useful in these applications is the fortu-
one-sided Markov property. nate fact that the conditional distribution on the unobserved
An example with is in Fig. 1. The graph variables given the observed variables (i.e., the posterior
summarizes many (conditional) independence relationships. distribution) is an MRF on the subgraph of obtained by
For instance, restricting to “unobserved” sites. The conditional distri-
bution on given
and —see Fig. 1—is Markov wrt the subgraph at sites
which in this case happens to be linear.
Up to a constant, conditional distributions are just joint distri-
is Gibbs with respect to a graph if can be butions with some of the variables fixed, so the statement about
represented as their dependency structure follows immediately from the Gibbs
representation.
(2) Hidden Markov models (speech recognition [23], Kalman fil-
ters [28], etc.), and hidden Markov random fields [29], have de-
pendency graphs like those in Fig. 2, where we have labeled the
where
observable variables of the model using ’s, and the unobserv-
• is the set of cliques in , i.e., the set of fully connected able ones using ’s, in order to distinguish them. In Fig. 2(a) and
subsets of (including singletons); (b), the posterior distribution, , is Markov wrt the linear
• , ; and graph; in Fig. 2(c) it is Markov wrt the nearest neighbor lattice
graph.
• each is a positive function.
The goal is usually to make an inference about , given .
Such representations are certainly not unique: constants can be Therefore, it is significant that the graphical structure repre-
“shifted” (multiply one term, divide another) and can be ab- senting the distribution on given is no more complicated
sorbed into whenever . than the original structure, since, as we shall see shortly, this
The main tool for working with MRFs is the representation structure determines computational complexity. In this regard, it
theorem of Hammersley and Clifford [26], [27]: is MRF with is reassuring that models such as those in Fig. 2 form a very rich
respect to (wrt) if and only if is Gibbs wrt . One direction class: the set of marginal distributions on , obtainable from
is easy: it is straightforward to verify that (2) has the required finite-state space hidden nearest neighbor Markov models is es-
Markov property. But the other direction (MRF wrt Gibbs sentially everything (up to arbitrary approximation, see [29]).
wrt ) is not easy. It does, though, make easy the proof that One way to understand this is to examine the graphical structure
MRFs wrt linear graphs are Markov processes (i.e., that (1) im- of the marginal distribution on . Consider the Gibbs represen-
plies the more familiar “one-sided” Markov property). tation: when the variables are summed (integrated) out, new
Referring to Fig. 1, the cliques are the singletons, the pairs cliques are introduced. Two sites, and , in the subgraph will
, , , , , , , , be connected under the marginal distribution on (after inte-
, , , and the triple . MRF wrt grating out ) if and are already connected under the joint
the graph in Fig. 1 means that can be factored into terms distribution, or if there exists a path, traveling strictly
each of which depends only on the components represented in through the variables, that connects and . So the marginal
one of these cliques. An analogous relationship appears in the on , in the cases depicted in Fig. 2, will in general define a
study of Bayes nets, which use directed acyclic graphs (DAGs). fully connected graph! (Just check that there is always a path
If “respects” a DAG (if factors into conditional probabili- 1Another variant is the Tanner graph (cf. [15], [25]), in which clique func-
ties of individual daughter nodes given their parent nodes) then tions are represented explicitly as specially designated nodes in the dependency
is Markov wrt the corresponding undirected “moral” graph graph.
552 IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 47, NO. 2, FRBRUARY 2001

. In general, has the factorization

where, for convenience, we have absorbed sub-cliques into

super-cliques. Following the site ordering, address first:
compute

and

Notice that “isolates” from the other components, and

is the value of that participates in the most likely
configuration, if that configuration involves . is
the corresponding contribution from terms involving , eval-
uated at . Since is next on the list, and isolates
and from the other components, compute

and

to get, respectively, the value of participating in the optimal

configuration (given ) as well as the contribution from terms
Fig. 2. (a) and (b) Hidden Markov models, and (c) a hidden Markov random involving and . As for
field .

“through ” connecting any two components of .) This ob- and

servation comes (easily—again!) from the Gibbs representation.
(The condition for creating neighbors in a marginal distribution
is necessary, but not quite sufficient—hence the disclaimer “in So far this is standard dynamic programming, but calls for a
general.” Consider, for example, slight generalization. Since it takes both and to isolate ,
, , and from the remaining variables, compute

Then and are marginally independent and

Proceeding in this way, always isolating what has been done

yet they are neighbors in the joint distribution and they are con- from what has not been done, we compute
nected by a path through .)

C. Computation
1) Most Likely Configurations: Computing a most likely se-
quence of words given an acoustic signal, a most likely image
restoration given a corrupted picture, or a most likely code-
word given a channel output, means computing a most likely
configuration under a posterior distribution. In light of our re-
marks about the dependency structure of conditional distribu-
tions (Section II-B), it is evident that the generic problem is
to maximize a probability distribution that is Gibbs relative to
some given graph .
Consider again the simple example in Fig. 1. If, for instance,
then there are possible configurations
and exhaustive search for the most likely one is impractical. On
the other hand, we could use a kind of (generalized) dynamic
programming: choose first an ordering of the nodes in , say
GEMAN AND KOCHANEK: PROGRAMMING AND GRAPHICAL REPRESENTATION OF ERROR-CORRECTING CODES 553

is the value of that participates in the most likely con-

figuration (with the corresponding cost ). Ev-
idently, then, is the corresponding value of , and
the corresponding value of , and so on,
back through the graph.
The most likely configuration is thereby computed with
many fewer operations than required in a systematic (brute-
force) search. How many operations are needed? The worst of
it is in the computation of , , , and , each of which
involves a triple loop (e.g., with regard to : for every
and every find the best ) and hence order opera- Fig. 3. Star-shaped dependency graph.
tions. Since there are 10 sites, and no site requires more oper-
ations than an upper bound on the computational cost 2) Sometimes, an will produce a tie. These can be
is . decided arbitrarily, in which case one of possibly many
This procedure generalizes. For any graph with sites maximum probability configurations will be found.
let be a “site visitation” schedule. How 3) There is a more-or-less obvious modification that finds
much work is involved in “visiting” a site? When visiting , the most likely configurations, for any .
we need to compute the best for each possible configura- 4) Often an optimal ordering, in the sense of achieving the
tion on the set of sites that isolate from the re- minimum , is transparent, more or less by inspec-
maining sites. In other words, if is the number of neigh- tion. But the problem in general, for arbitrary graphs, is
bors of among in the graph (the NP-hard. See Arnborg et al. [30].
“size of the boundary set of ”), then there are
operations associated with the visit to ( neigh- 2) Marginal and Conditional Probabilities: Often, the ob-
bors plus a loop through the states of ). Therefore, if ject of interest is a marginal or conditional probability distribution
, then the maximum-likelihood config- on a subset of the variables. What is the probability of conges-
uration can be found in tive heart failure in a 50-year-old male with swollen ankles and
pneumonia? What is the probability that a particular decoding,
for example, the maximum-likelihood decoding, is actually cor-
rect given the output of a noisy channel and given that codewords
are (say) a priori equally likely? These conditional probabilities
operations. are quotients of marginal probabilities—probabilities on config-
Actually, things can be better than this. Suppose urations of subsets of the variables. The coding example requires
is not connected (in the graph structure of ). The the marginal probability of the observed channel output, and this
maximization on is not necessarily dependent upon all involves a summation of probabilities over all possible inputs; for
of the bounding variables of the set . In fact, a medical application, we may need to compute marginal prob-
one need only fix the values of the variables at those sites abilities on small subsets of variables, associated with diseases
that isolate the particular connected component containing like congestive heart failure and pneumonia, attributes like age
. Thus, should be interpreted, more favorably, as the and sex, and signs and symptoms like swollen ankles.
largest boundary of a connected component created by the site The general problem is, therefore, to compute the probability
visitation schedule . of a configuration of states for a subset of variables, given a
This makes a substantial difference. In Fig. 3, with visita- Gibbs distribution on an associated graph . There is again a
tion schedule no connected set is gen- dynamic programming principle, operating in much the same
erated with boundary size greater than one: and only way as in the calculation of a most likely configuration. This is
operations are needed. Of course order matters: any probably best illustrated by example.
visitation schedule that starts with will immediately incur Consider again Fig. 3, and suppose we wish to calculate the
computations. marginal distribution on . Then for each value of and
we will need to sum out the other seven variables
Remarks:
1) None of this would really work on a big problem, say with
, even if were small. The functions rep-
resent probabilities, and when is large, the probability
of any configuration is exponentially
small and would generate an underflow. Therefore, in
practice, we maximize the logarithm of instead of it-
self. Products (such as above) are re- The apparent computational cost is , but if we pay attention
placed by sums (such as ), to the order of summation then this can be reduced to less
but otherwise the procedure and the reasoning behind it than . We again define a site visitation schedule, say
are the same. and t his a gain dictates the sequence of
554 IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 47, NO. 2, FRBRUARY 2001

calculations: define The inverse, , is then

Our only interest is in efficiently summing over , so let us

make things transparent by fixing and and writing

in which case the problem becomes the evaluation of

and then finally The clique structure is the same as we started with, and therefore
so is the dynamic programming principle and the number of op-
erations. This time, however, there are no numerical problems.
Consequently, we will take this approach when, in Section IV-B,
we compute some posterior probabilities of maximum-likeli-
In essence, the two schemes, for maximizing and for sum- hood decodings.
ming, are the same. Pick a visitation schedule, fix the variables
on the boundary of the connected set containing the current site, III. LINEAR GRAPHS
and then either maximize or sum. In either case, the number of
elementary computations is no worse than Linear dependency graphs come up in many applications:
speech recognition, convolutional coding, filtering, and control,
among others. In general, there is an observation vector and a
As we have said, conditional probabilities are quotients of mar-
“state vector” , and a joint dependency structure like the one
ginal probabilities, so conditional probabilities are also amenable
in Fig. 2(a) or 2(b).
to dynamic programming. But this will not always work! At least
In a speech recognition system, might represent a portion
not when there are a large number of “observed” variables—vari-
of a phoneme uttered as part of a word, or perhaps even a pair
ables upon which we condition. The problem is again numerical:
of words, so that the state space is potentially quite large, repre-
the a priori probability of any one configuration of the observable
senting the word or pair of words in addition to the phoneme and
variables is exponentially small, but this very same probability
phoneme fraction. The observable is some representation or
is the denominator of the quotient representing the desired con-
encoding of the associated acoustic signal, or more precisely, the
ditional probability. If, for instance, a block code transmits 1024
signal as it has been recorded by the microphone. In the speech
bits, then the unconditioned probability of receiving any partic-
application, this particular dependency graph comes out of the
ular 1024-bit word is hopelessly small—much too small to be
much-used hidden Markov model, under which
computed with a summation scheme like the one recommended
above. On the other hand, the conditional probability of, say, the
most likely transmitted word, given the received word, will typi-
cally be order one. What we are after is a ratio of two very small
numbers, and we need to use caution to avoid gross numerical and
error.
One way around this is to mix the computations of the numer-
ator and denominator in such a way as to avoid exponentially
small terms. It turns out that this is easiest to do if we com- or
pute the inverse of the conditional probability, rather than the
conditional probability itself. To illustrate, let us write for the
“unobservable” components and for the “observable” compo-
nents and for the complete vector of variables. The Gibbs The joint distribution is clearly Gibbs with respect to the graph
distribution has the form in Fig. 2(a) or 2(b), depending on which model is used for
.
The object of interest is, of course, the configuration , and as
and, therefore, we have already noticed, its posterior distribution, given , cor-
responds to a simple linear graph. So the computational com-
plexity, whether computing a MAP sequence or
the posterior probability of such a sequence, is no worse than
(using a left-to-right or right-to-left visitation schedule).
where is the dimension of . Furthermore, in many applications, most transitions
GEMAN AND KOCHANEK: PROGRAMMING AND GRAPHICAL REPRESENTATION OF ERROR-CORRECTING CODES 555

are ruled out, a priori, meaning that for most

. This can help substantially: if the software can be
structured to efficiently ignore the zero transitions, then instead
of operations, we can expect something more like ,
where is the “arity” or number of possible transitions from a
state . m = 2).
Fig. 4. Second-order posterior dependency of a convolutional code(

A. Convolutional Codes and Memoryless Channels Thus, the dependency of on has the form
Convolutional codes also give rise to linear dependency
graphs, though the neighborhood structure is generally richer with the understanding that depends only on
than the nearest neighbor system of first-order Markov pro- the last arguments when .
cesses. We shall formulate, here, convolutional decoding as Suppose the codeword goes through a channel and oc-
an instance of dynamic programming for general graphs, and casionally gets corrupted. Let be the th bit of the th code
recover the well-known Viterbi algorithm [5], [31]. There is block and let be the cor-
nothing new in this exercise (in particular, see [7] and [9]), but responding output bit. The usual channel model has the form
our general viewpoint does suggest some extensions that may
be of some practical value. One could, for instance, perform
exact maximum-likelihood decoding even when the channel where are independent and identically distributed (i.i.d.), so
noise is not white. Such Markov dependency within the error that
process might afford a good model of bursting (see [20]). In
this case, the dynamic programming principle still holds and
the maximum-likelihood decoding is still computable, at a
modest increase in computational cost. Furthermore, the exact Hence,
posterior probability of the maximum-likelihood decoding is
also computable, for about as much additional computation as
was used for the decoding itself.
Recall that an convolutional code is defined
through a generator matrix of the form

.. .. .. ..
. . . .

where each is a submatrix. For convenience, we will where

stick to binary codes, so that the elements of are in .
To maintain the connection with Section II, we introduce
the following (somewhat unusual) notation: will be the
upper-left submatrix of , The clique structure, and hence the computational complexity
will represent the th block of information bits , of dynamic programming, would be unchanged by any of a
and will represent the corresponding block of rich collection of prior models on , governing the arrival of
code (output) bits. If and , information bits. Since under the uniform prior, the posterior
then . distribution is proportional to , maximum-likeli-
Since cliques determine computational complexity, it is hood decoding is the same as MAP decoding, and the (maximal)
useful to observe that cliques are of the form . For example,
the dependency graph for the case is depicted in Fig. 4.
The obvious site visitation schedule is which
incurs a maximum boundary of sites and hence a maximum
computational cost of ( since )
.. at any one site. With sites, the number of operations needed to
.
compute the MAP (or maximum-likelihood) decoding is, there-
.. fore, . This is, of course, the well-known Viterbi
. algorithm.
Suppose turns out to be the maximum-likeli-
hood block of information bits. It would be a good idea to com-
.. pute the associated probability that is correct: . As
. we have already seen more generally in Section II-C-2, these
probabilities are computed by a straightforward adaptation of
556 IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 47, NO. 2, FRBRUARY 2001

the same dynamic programming (Viterbi) algorithm that gave IV. PRODUCTION SYSTEMS AND TREE-STRUCTURED GRAPHS
us in the first place.
From the computational viewpoint, the primary virtue
B. Burst Errors of linear graphs is that the computational cost of dynamic
programming grows linearly with the number of variables,
One way to model burst errors is with a Markov noise process even while the configuration space grows exponentially. More
where and represents general lattice graphs, such as with and nearest
a transmission error ( or ). The typical state is neighbor interactions, behave differently. A
presumably “ ,” but an occasional transition occurs to “ ,” and sublattice of will achieve a maximum boundary of at least
there is some tendency to stay at “ .” If , no matter what the site visitation schedule.
Computation, therefore, grows exponentially (in ) whenever
.
then the situation can be modeled by making (the probability In between the linear graph and the lattice graph (with )
of initiating a burst) very small, and , being the are tree-structured graphs, which fortunately also admit site vis-
average burst length. The channel model is then completed by itation schedules with bounded maximum boundaries. As an
introducing an initializing probability , in example, consider the tree-structured (but cyclic) dependency
which case the model for is graph on in Fig. 5.
Label the sites at level , from left to right, by
and consider the “bottom-up, left-to-right” site visitation sched-
ule: The largest
Many elaborations are possible (as developed, for example, boundary encountered for any connected component has size
in [20], and connected to graphical models in [7]), including, , and this is independent of , the depth of the
for instance, state-dependent bursting in which statistics of tree. Since there are nodes, the number of dynamic
depend on the transmitted data , but the Markov model em- programming operations for computing probabilities and
bodied in and is sensible and, in any case, suitable for most likely configurations of associated Gibbs distributions
illustration. is . Thus computation grows linearly with the
What are the implications for computing maximum-likeli- number of variables for Gibbs distributions on tree-structured
hood (or MAP) decodings? Introduce the error indicators graphs.
if Because of their computational advantages, tree-structured
else. dependencies are attractive modeling tools (e.g., [32], [33]).
They also come up naturally when working with production
Then systems, which define the so-called “context-free” grammars
studied in formal linguistics ([4]).
In this section, we will introduce a suite of error-correcting
codes that are based on, or in any case admit representations
in terms of, production systems. (As we shall see, the approach
turns out to be nothing more than a reformulation of Forney’s
“Squaring Construction,” [21]. See also Gore [34], for an ear-
lier but less developed squaring-type construction.) In compu-
tational linguistics, the range of the (vector of) leaf-node vari-
ables is known as the “yield” or “language.” In our application,
the range is a set of permissible codewords rather than a set of
well-formed sentences. In either application, whether to linguis-
tics or coding, the tree structure is exploited to design efficient
computational algorithms.
In way of illustration, let us examine some of the computa-
tional consequences of a formal grammar representation of the
even-parity code. A context-free grammar (in “Chomsky normal
form”—see [4]) consists of a finite set of nonterminal symbols,
, a start symbol , a finite set of terminal symbols ,
in light of the relation between and . and, for every , a finite set of production rules, each of
Evidently, then, the situation is not much different from the the form
simple i.i.d. channel model. This time, the maximal cliques have
the form which is an expansion
over the i.i.d. model by only one site, and therefore the most
likely decoding and its posterior probability can be computed or
with operations, or about times the decoding
cost under a white-noise model.
GEMAN AND KOCHANEK: PROGRAMMING AND GRAPHICAL REPRESENTATION OF ERROR-CORRECTING CODES 557

Fig. 5. Tree-structured dependency graph. The leaf nodes of the X -graph represent a codeword, from a context-free code, and the Y nodes represent the output
of a memoryless channel.

In general is allowed, but a reduction to read and two more productions are applied (to the level
Chomsky normal form is always possible (again, see [4]). Typ- states), thereby determining the states of the four-level
ically, there is a multitude of production rules for each . sites. Encoding continues through the information bits,
The language, or yield, of the grammar is the set of strings of resulting in a specification of states at every site in the graph.
terminals that can be derived from through repeated appli- In this manner, a length even-parity codeword is produced at
cation of the production rules. Probabilistic grammars include the leaf nodes, inducing a one-to-one correspondence between
a collection of probability distributions, one for each , sequences of information bits and length even-parity
that dictates the choice of production rules. This induces a prob- codewords. A simple example, with , is given in Fig. 7,
ability distribution on the yield. where information bits have been placed in the center of the tri-
In linguistics, usually denotes a parsed sentence, where angles, and where the conventions “bit apply first pro-
the nonterminals represent sentence fragments, such as noun duction” and “bit apply second production” have been
phrase, verb phrase, prepositional phrase, article, noun, and so used to translate from information bits into states, via equations
on, and the terminals represent lexical items, typically words. (4) and (5).
But suppose, instead, that , , , The random variables are , , ,
and the production rules are with state space when , when ,
and when . The production rules induce clique
(3) functions—one for each triangle. In most instances, the “natural
prior” is the uniform prior (in which case maximum-likelihood
Then, evidently, the yield is exactly the set of all nonempty, decoding is MAP decoding), which results, evidently, from a se-
finite, binary strings with an even number of ones. quence of i.i.d. “balanced” information bits. The cor-
We can connect this to Markov random fields, and block responding clique functions are if
codes, by fixing a binary graph structure, assigning to and otherwise ( is always ), for the apical
the root node, assigning terminal states to the leaf nodes, and (root) node triangle; and if is
assigning nonterminal states to the nonleaf nodes. Consider, for the parity of and otherwise, for the remaining
example, the balanced binary tree in Fig. 6. Instead of (3), we triangles.
adopt the production rules In summary, Fig. 6 depicts an MRF (“ ”) in which the mar-
(4) ginal distribution on leaves concentrates on the even-parity code
of length ; the marginal distribution depends on the distribu-
at levels , and tion on information bits; and in particular, the marginal distri-
(5) bution is uniform when the distribution on information bits is
uniform.
at level . A sequence of information bits can be turned
Imagine now a channel, let us say memoryless for the
into a state configuration on the graph by assigning each bit to
time being, through which the leaf-node variables are trans-
one of the triangles in the tree. The bit associated with the
mitted and possibly corrupted, and eventually received as
apical (root) triangle is used to choose between the productions
. If, for example,
and . This fixes the states of the daughter
sites (at level ), and then two more information bits are
558 IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 47, NO. 2, FRBRUARY 2001

Fig. 6. Even-parity code viewed as the yield of a context-free grammar.

transmission across a memoryless channel. In Section IV-C, we

will examine the extensions to Markov models for bursty chan-
nels.

A. Context-Free Codes
There are many ways to generalize. Here we will stick to
binary codes, balanced, binary, tree-structured graphs, and a
symbol set (state space) that depends only on level . We will
assume that the number of productions per symbol also de-
pends only on the level, and furthermore we will restrict our-
Fig. 7. Even-parity code. The information sequence 1100010 is coded as
01001101.
selves to block codes with minimum distance , for some
. If the number of productions per symbol is al-
ways a power of two, then the encoding scheme that we used
with i.i.d. , for some variance , then for the even-parity code generalizes directly: just “peel off” the
number of bits needed to specify a production (this depends only
on the level—not on the state), working down from the top. Ob-
viously, there are other, less restrictive ways to map informa-
and the joint dependency graph has the structure depicted in tion bits into productions, but we will confine ourselves to this
Fig. 5. simple case by devising production rules in which the number
It is immediately evident that the posterior distribution has the of choices depends only on the level and is always a power of
same structure as the prior, and hence, in light of our earlier dis- two.
cussion about dynamic programming on trees, the most likely Look again at Fig. 6. At every (nonterminal) node, each
(or MAP) even-parity word can be computed in “linear possible state can itself be thought of as representing a code;
time”—the number of operations grows linearly with the di- namely, the set of derivable terminal sequences of the subtree
mension of the code. Of course, there are other effi- defined by the chosen node. Site , for example, is the root
cient ways to “soft-decode” the even-parity code, but bear in node for the terminal sequence of sites ,
mind that the “reliability,” , is also calculated for about and (resp., ) if and only if has
the same cost. even (resp., odd) parity. In this way, an at level represents
Obviously, it would be desirable to extend this representation the even-parity code of length , and an represents the
beyond the even-parity code. One way to do this is to generalize corresponding odd-parity code.
the production system ((4), (5)) that generates the even-parity More generally, let be the non-
code. In Section IV-A, we will formulate production systems terminal symbols (states) of a level- site. (Later, will depend
that yield other (usually familiar) error-correcting codes. Then, on , .) If we want a code with minimum distance ,
in Section IV-B, we will study the computational implications then, evidently, the yield of , as expressed at the terminal sites
for computing maximum-likelihood (or MAP) decodings, and must itself be a code with minimum dis-
for computing the probability that the decoding is correct, given tance no smaller than . Taking a hint from the even-parity
GEMAN AND KOCHANEK: PROGRAMMING AND GRAPHICAL REPRESENTATION OF ERROR-CORRECTING CODES 559

code, suppose that each represents a length , distance , The correspondence is built inductively. Start by representing
code (just as and at level represent length , distance with the symbols
, codes), and that the distance between any two of these
codes, is at least (just as the even- and
odd-parity codes are distance apart). Then the pro- Now build the representation for from an already-built rep-
ductions resentation for , through the formula

define level symbols , each (6)

of which represents a length , distance , subtree code,
and, furthermore, whenever . where integer and
Each code is just a union of Cartesian products, so the . In way of example, let us work through the
construction could be written more traditionally as particulars when and

where represents Cartesian product (concatenation)

and equals the matrix shown at the bottom of this page.
The representation is useful because it makes explicit a hier-
What we have is just one example of Forney’s (twisted) archy of distance properties, as can be verified inductively
“squaring construction” [21]. The even-parity code is a special
case with , , and . There is nothing (7)
special about the particular “twisting” , and, in whenever . An immediate consequence is that
fact, other ways of combining and are possible and often
desirable, as we shall see shortly. The point for now is that such
constructions can be iterated, “bottom up,” thereby defining since each is evidently distinct and there are just the
symbols and productions all the way up to the root node at right number of them. Beyond this, the distance properties lead
level (where any of the states more-or-less directly to a hierarchy of codes and separations (as
could be used as a start symbol). These productions define we shall see momentarily), and this hierarchy, in turn, leads to
a tree-structured dependency graph, and through this, a the representation of various codes in terms of production rules.
ready-made computational engine for soft decoding and other The reader may recognize the representation as a version of
decoding calculations. Forney’s iterated squaring construction [21], albeit with a new
Of course, the construction assumes that we start with a set notation. The combinatorial function is not partic-
of level- codes, , with certain distance properties and certain ularly special, and many others could be used (corresponding
separations. One way to get at this is to generalize even further: to different “twistings”). In [21], Forney introduced an espe-
replace the set of separated codes by cially canonical twisting (see Appendix A) for which we will
a hierarchy of codes with a hierarchy of separations. This more later develop an exact coarse-to-fine version of the dynamic pro-
elaborate representation falls out of a hierarchical representation gramming algorithm (in Section V). But let us first develop the
sought-after grammatical representation, since this is indepen-
for itself, which we will now construct.
dent of the particular twisting (and, hence, combinatorial func-
The elements of will be represented by symbols tion) used.
. The integer coefficients are restricted Right away we get a hierarchy of codes and separations, just
by , where , for each . by taking unions over indices: for each and each
There are, then, vectors ; we will , define
set up a one-to-one correspondence between these and the (8)
elements of .
560 IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 47, NO. 2, FRBRUARY 2001

Evidently, in light of (7), is a distance code (with the simpler model corrupts bits independently (memoryless
codewords), and furthermore channel)

(9)

whenever . What is more, the construction of where represents the received signal and
from (see (6)) translates into a construction of is an independent noise process. In this case
from . If then

and we get the joint dependency graph depicted in

Fig. 5.
We ran some experiments with the RM code
, transmitted at low signal-to-noise ratio through a
Gaussian channel: . Following
tradition, we have assumed BPSK modulation
We have, therefore, for any , a production system . For the memoryless channel, the dependency graph
of the posterior is the same as the dependency graph
of the prior , and is tree structured under the grammatical
representation of RM . See Fig. 5, where for the purpose
of organizing the computation we can simply ignore the sites
belonging to the received vector . In each of 50 trials, we
computed both the maximum-likelihood decoding , and
for the associated “reliability” . As suggested earlier, we
used the “bottom-up” site visitation schedule ,
, , , which makes for
about 80 000 operations per decoding cycle and per posterior
probability calculation. In computing , we were careful
to avoid underflow, by following the recipe given at the end of
Section II-C.2.
for productions are ``singletons'' The results are summarized in Table I, where represents
the correct decoding and is the Hamming distance to
(10)
the maximum-likelihood decoding. Notice that the single de-
which, along with the conventions and , coding error is well anticipated by the low posterior probability
yields a distance code for every and every .
, . The meaning of the
case is that, from level down, no more information C. Burst Errors
bits are encoded; the codeword is already determined by the What are the computational implications of a channel model
configuration at level . with Markov memory? Let us take another look at the model
We started with the even-parity code. To recover the even- for burst errors discussed earlier in Section III-B, but this time
parity code, as a special case, just take and identify in the computational context of a tree-structured prior.
with and with for all . The channel model introduced in Section III-B represents
RM Codes: A Canonical Class of Context-Free bursts by a Markov error process such that
Codes: Since, as Forney has shown [21], the RM codes . In terms of the initial probability distribu-
can be derived from an iterated squaring construction, it is tion on , and the transition probability ma-
natural to seek a grammatical representation for them. It turns trix for Prob , the channel
out that the RM grammar is identical to the grammar presented model can be written
above, with the sole exception being that a different choice of
combinatorial function (or twisting) is required. We present
a brief derivation of this fact in Appendix A. In Section V,
we will exploit the grammatical structure of the RM codes to where if , and otherwise. Of course, there
develop a more efficient maximum-likelihood decoder based are other models, that perhaps make more sense in a particular
on the notion of coarse-to-fine dynamic programming. application ([20]), but the computational costs will be the same
or similar for many of these variations.
B. Memoryless Channels Since , we get the joint dependency
The random variables situated at the leaf nodes, graph, for , by starting with the (tree-structured) depen-
, make up the codeword—the bits that are actually dency graph for , and appending with sites for and with
transmitted. Returning to the channel models of Section III, cliques introduced in the channel model, . Since is a
GEMAN AND KOCHANEK: PROGRAMMING AND GRAPHICAL REPRESENTATION OF ERROR-CORRECTING CODES 561

TABLE I
RM (2; 6) WITH BPSK MODULATION AND ADDED WHITE GAUSSIAN NOISE. TYPICAL DECODINGS AND THEIR POSTERIOR PROBABILITIES. FIRST COLUMN IS
TRIAL NUMBER; SECOND COLUMN IS HAMMING DISTANCE BETWEEN THE MAXIMUM-LIKELIHOOD DECODING AND THE CORRECT CODEWORD; THIRD COLUMN
IS THE PROBABILITY ASSIGNED TO THE DECODING UNDER THE POSTERIOR DISTRIBUTION. AT  = 1, MAXIMUM LIKELIHOOD YIELDED ONE DECODING
ERROR (IN TRIAL #3) IN 50 TRIALS. NOTICE THAT THE LOW POSTERIOR PROBABILITY (0:3790) SIGNALS AN ERROR. THE NEXT
LOWEST PROBABILITY WAS 0:5866, BELONGING TO A CORRECT DECODING

Fig. 8. Tree-structured code through a Markov channel. (a) Joint (X; Y ) dependency structure. (b) Dependency structure on X , under P (xjy ).

function of and , and since is a function of and nodes. An optimal ordering of site visits is no longer obvious,
, introduces a clique made up of the four sites asso- and may, in fact, depend upon the various sizes of state spaces
ciated with the four variables and . Taking at the various sites. But there is a particular ordering that
this into account, and taking into account the tree-structured does the job at a cost of no more than 16 times the cost of
prior on , we arrive at the joint dependency struc- decoding the simpler memoryless channel analyzed earlier.
ture depicted in Fig. 8(a). As we have noted before, conditioning Before writing down this particular visitation schedule, which
does not change the graph structure. The posterior dependency may appear at first glance to be rather complex and arbitrary,
graph comes from simply removing the sites representing the it might be better to first describe, less formally, the resulting
received signal as depicted in Fig. 8(b). dynamic programming process.
Compare the dependency graph on in Fig. 5 to Fig. 8(b). The optimization (soft decoding) is based on subtrees. The
Channel memory introduces new neighbors among the leaf subtree rooted at ( and ) is made up
562 IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 47, NO. 2, FRBRUARY 2001

of and the two “daughter” subtrees rooted at and . have at least operations, no matter in what order sites
Suppose that for each of the two daughter subtrees the optimal are visited, and interesting lattices begin with or even
“interior” has been calculated, which is to say the optimal as- . But even a good graph structure does not guarantee
signments of states conditioned upon all possible assignments at a feasible computational problem. Sometimes the state spaces
the triangles’ respective corners—the root ( or ) and, are too large. The codes of Sections III and IV have good graph
for each root, the corresponding pair of corners sitting at level structures but in some instances very large state spaces.
( or ). Then How large are the state spaces under the grammatical repre-
the optimal interior of the level- subtree can be computed by sentations of the RM-like codes introduced in Section IV? Fix
“merging” the daughter subtrees. The merging involves visiting (code length ) and fix (code distance ). The com-
(maximizing over) the rightmost corner of the left daughter sub- putation of the number of states at a particular site does require
tree, the leftmost corner of the right daughter subtree, the root some care, since not every available symbol is actually used. At
node of the left daughter subtree, and then, finally, the root node level there are symbols, but only one is actu-
of the right daughter subtree. Merging continues, “upwards,” ally used—the designated “start” or “sentence” symbol . Still,
until eventually the optimal interior of the entire graph is calcu- the number of symbols employed at a given site is independent
lated, conditioned on the configuration at the graph corners , of the particular start symbol, and in fact depends only on the
, and . The global optimum is then computed after a max- level, (see [36])
imization at each of these three remaining sites.
states of a level- site
The procedure is summarized more formally through a set of
do-loops defining the visitation schedule (see the bottom of this
page).
In effect, every production from incurs a 16-fold increase
in computation cost: the two remaining interior leaf nodes (right
This leads to large state spaces, even for more-or-less
node of the left daughter and left node of the right daughter),
modest values of . The code RM , with
with four configurations, are visited once for each of the four
(length ) and (distance ), generates 2,
configurations of the two leaf nodes that sit at the corners of the
4, 16, 256, 16 384, 1 048 576, 1 048 576, states, at levels
subtree.
respectively. The worst of it, from
As usual, the ordering (and computational analysis) applies
the dynamic programming viewpoint, would appear to be the
equally well for probability calculations, including the posterior
loops associated with the level – cliques, which are triangles
.
involving one level- site and two level- sites. Naively, there
are operations in this loop, and this is, of
V. THINNING AND EXACT COARSE-TO-FINE course, infeasible.
Dynamic programming can get out of hand, computation- But the calculation is naive since the clique function is zero
ally. The graph structure, for example, may be bad: there may for the vast majority of these triples. Nonzero contributions are
be no site visitation schedule that avoids large boundaries. The made, only, by those pairs of level- symbols that represent an
two-state two-dimensional Ising model, on the square allowed production from a level- symbol. The number of op-
lattice, is the prototypical example (cf. [35]). The worst loops erations associated with a triangle is the number of symbols at

Do = 2; p
l % loop over levels
Do i = 1; 2p0l % loop over sites at level l
% Create a connected component of “interior” sites of the subtree
% rooted at sil by “merging” the subtrees rooted at s2l01 l01
i01 and s2i
0
Visit s(i01)2 +2 % visit the right-most corner
% of the subtree rooted at s2l0
i01
1

0
Visit s(i01)2 +2 +1 % visit the left-most corner
% of the subtree rooted at s2l0 1

l01
i
Visit s2i01 % visit the root node
% of the left subtree
Visit s2l0 i
1
% visit the root node
% of the right subtree
End do
End do
% All that remains are the corners of the tree …
Visit s10
Visit s20
Visit s1p
GEMAN AND KOCHANEK: PROGRAMMING AND GRAPHICAL REPRESENTATION OF ERROR-CORRECTING CODES 563

TABLE II
THE PERFORMANCE OF CTFDP RELATIVE TO DP. FOR SIX REPRESENTATIVE THINNED RM CODES, THE TABLE COMPARES THE NUMBER OF DECODING
OPERATIONS REQUIRED BY EXACT DP WITH THE AVERAGE NUMBER OF OPERATIONS FROM A SERIES OF CTFDP DECODING SIMULATIONS. CODEWORKS WERE
TRANSMITTED ACROSS A MEMORYLESS AWGN CHANNEL (STANDARD DEVIATION  ) WITH BPSK MODULATION

the apex level times the number of productions per symbol. If, carding productions (and associated information bits) to reduce
therefore, we loop first over productions, given a symbol at the decoding complexity. What sort of code would emerge from
apex, and then over symbols, then the number of operations will this process? Alternatively, one could imagine imposing strict
be , where is the number of productions per level- limits on the cardinality of each level’s state space and inquiring
symbol. This latter number is independent of both the partic- whether the resulting system remained a consistent context-free
ular symbol and the particular level- site, as is evident from the grammar. These equivalent approaches yield a family of con-
production formulas. For the RM codes, for text-free codes that we will refer to as thinned codes.
and for , so that in the particular example We present a brief introduction to thinned codes in Ap-
there are 1, 1, 1, 2, 16, 1024, and 1 048 576 productions for levels pendix B (though a fuller treatment is available in [36]).
and respectively. The thinned RM code RM , defined to be a
The product is biggest when : each of the two subcode of RM with the number states (or pro-
level- sites contributes about opera- ductions) at any level not exceeding , is readily decodable
tions. The cost of decoding, or of evaluating a posterior proba- by the exact maximum-likelihood decoding algorithms of
bility, is about Section IV in at most (loosely) operations. For example,
RM , a linear subcode of the un-
decodable code RM , is decodable in
approximately operations. In other words, by discarding
only half the information bits from RM we can de-
since is the number of sites at level , and this is about code the remaining thinned code at least times faster.
operations. This may be feasible, but it would be impractical in Moreover, using the coarse-to-fine approach of Section V-B,
most applications. RM can be decoded an additional three to 30 times
We will suggest a few remedies. One (coarse-to-fine dynamic faster depending on the signal-to-noise ratio (see Table II).
programming) is more or less generic, in that it applies in prin- As a general rule, thinned RM codes are poor codes in terms
ciple to any problem of finding a most likely configuration under of coding gain and other performance measures. However, they
a probability with a given graphical dependency structure. An- are useful in the context of context-free codes, because they
other (thinning) is special to the problem at hand: decoding the allow one to vary decoding complexity (often by orders of mag-
grammatical codes introduced in Section IV. nitude) by simply altering a single parameter. One speculative
direction for future inquiry is the following problem. Given a
A. Thinning thinned RM code with a known and manageable decoding com-
As we have just seen, the computational cost of max- plexity, can one find a set of combinatorial functions (or twist-
imum-likelihood decoding or computing a posterior probability ings) that optimize the coding gain of the iterated squaring con-
for a typical context-free code like RM trans- struction?
mitted across a memoryless channel is a simple function of A far more promising approach than thinning for the max-
the state-space cardinality and the number of productions. imum-likelihood decoding of context-free codes is the method
Although the coarse-to-fine dynamic programming (CTFDP) of CTFDP.
algorithm (see Section V-B) can reduce decoding complexity,
sometimes by several orders of magnitude, the really large con- B. CTFDP Programming
text-free codes are still undecodable, at least from a practical CTFDP is what Pearl [37] would call an “admissible
point of view. heuristic,” meaning that it is a variation on dynamic program-
Consider for example the code RM . With an infor- ming that is meant to save operations in a typical problem (it is
mation rate of and distance , its (maximum-like- a “heuristic”), but, at the same time, it is guaranteed to solve the
lihood) decoding complexity exceeds the level– contribution optimization problem (it is “admissible”). The well known
of operations! RM is patently undecod- is an admissible heuristic, as is the iterated complete path (ICP)
able. But suppose one systematically pruned its grammar, dis- algorithm of Kim and Kopac [38], [39]. CTFDP is a kind of
564 IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 47, NO. 2, FRBRUARY 2001

Fig. 9. CTFDP. Darkened circles indicate a chosen path of super states. States along the chosen path are refined and dynamic programming is repeated. In panel
(d), dynamic programming yields a path of singletons, at which point the optimal path has been computed.

hierarchical, or multiresolution, version of ICP invented by C. 2)

Raphael [22]. We employ it here in order to reduce (sometimes
dramatically) the computational cost of maximum-likelihood
decoding. Unfortunately, the ideas behind the method do remember that and are singletons
not extend in any obvious way to the problem of computing Then we could perform dynamic programming on super states,
posterior probabilities. starting with the partition , and using the costs
The topology of the dependency graph is irrelevant, but the . Since presumably , this first
idea of CTFDP is perhaps best illustrated with a simple one- dynamic programming would be quick. We now refine: each
dimensional lattice. Let have joint probability (super) state along the optimal path (as chosen by the previous
distribution dynamic programming pass) is refined into its subset of super
states from the next lower level of coarsening

Now the state spaces are larger than at the first pass, but
and let have common finite-state space still, presumably, not nearly as large as . Another dynamic
(range) . The idea is to coarsen into a small number of programming pass generates another sequence of super states,
“super states,” and to perform multiple passes of dynamic pro- and the refinement/dynamic-programming cycle continues.
gramming on super states, successively refining super states on Evidently, a path will eventually consist only of single states
each pass. Formally, for each “coarsening” and, evidently, in light of properties (1) and (2) charac-
we define a partition of with elements in such a terizing the heuristic, this single-state path solves the original
way that is recovered at coarsening , and refines optimization problem.
. Fig. 9 describes how this process might proceed on a simple
1) ( partitions ): , and linear graph. Super states are delineated by boundary marks, and
the super states along an optimal path are indicated by darkened
;
circles. Here, , , , and the refinements are
2) ( refines ): for every and all binary
every there exists
such that
3) ( recovers ): .
So , are just the individual elements of and , are the individual states
. and . The path chosen in Fig. 9(d) consists only of
Now suppose that for any , and any and , we could single states, so the process ends here, meaning that
find a “heuristic cost” such that
1)
CTFDP may or may not find the optimal solution efficiently.
There seems to be a rather subtle relationship between the struc-
ture of the problem at hand, and the savings won (or lost!) in
and a coarse-to-fine implementation. As it turns out, in the case of
GEMAN AND KOCHANEK: PROGRAMMING AND GRAPHICAL REPRESENTATION OF ERROR-CORRECTING CODES 565

squaring constructions, the relationship is often highly favor- recompute the super state productions, solve the new DP problem,
able, and it is worthwhile, therefore, to look at generalizations and again examine the optimal derivation tree. If not, we stop:
beyond the simple one-dimensional lattice. the current optimal derivation tree represents the minimum-cost
It is clear enough how to proceed for more general graphs. codeword. Since the final derivation tree contains only states
Introduce a hierarchy of super states at every node, and define a from ’s own grammar, it certainly generates a codeword (in ).
heuristic cost for every clique function Moreover, this codeword is, by definition, the minimum-cost
super codeword in the final super code—a code that contains
itself.
where is a vector of super states, with one component super Although this CTFDP algorithm must eventually produce a
state for each site , and where means com- solution to the given optimization problem, it need not neces-
ponent-wise membership. Our first choice for a heuristic cost sarily outperform standard DP. For the procedure to converge
would naturally be rapidly, the number of refinements and subsequent DP compu-
tations must be minimal. This suggests that super states should
(11)
consist of aggregations of “similar” states so that their costs
if this were actually computable at a less than prohibitive cost. more closely reflect those of their constituents. In addition, the
As we shall see shortly (Section V-C), this “ideal heuristic” has a determination of super state productions must not be too com-
simple analytic representation for the (canonical) squaring con- putationally demanding. Remarkably, at least in the case of RM
struction, and hence for codes of the type developed in Section codes, one can, in fact, find a natural choice of super states that
IV as well as the variants of these discussed in Section V-A. addresses these concerns, resulting in a substantially faster max-
This is perhaps a little surprising, and it has the fortunate con- imum-likelihood decoder.
sequence that CTFDP for these codes is particularly convenient Using the partitioning scheme introduced in Appendix C,
and sometimes spectacularly efficient. we can implement a CTFDP version of maximum-likelihood
decoding for the RM codes. Table II presents the average ratio
C. Exact Heuristics of CTFDP to DP operations from 50 trials with each of four RM
Recall from Section IV that with grammar-based (squaring) codes and two thinned RM codes (see Section V-A). Except
constructions, clique functions are (up to a multiplicative con- for the very smallest code RM , the CTFDP algorithm
computes the maximum-likelihood codeword substantially
stant) just – valued, indicating un-allowed or allowed produc-
tions, respectively. The “ideal heuristic” suggested in (11), ap- faster on average with an efficiency increasing as the code size
plied to the clique , is then also binary, or signal-to-noise ratio increases. In the case of RM , the
coarse-to-fine procedure is five orders of magnitude faster than
simply indicating the existence of an allowed production within
the effectively impractical DP approach!
the coarsened states
APPENDIX A
if production RM GRAMMARS
with and In this appendix, we present a brief derivation of the RM
otherwise grammar. We refer the reader to Kochanek [36] for a fuller ac-
count.
for any three super states , , and at the three sites , , Originally in Section IV-A, our indexing system for gram-
and . It is useful to think of this as defining “super-state matical symbols was designed to reflect a geometric hierarchy
productions”: the production is allowed if there is a of minimum distances and separations. In the linear context,
corresponding production among the constituent states. both the indexing scheme and distance properties of symbols
With this choice of super-state productions for a context-free emerge naturally from the algebra of RM codes. Each binary
code , any given set of state-space partitions on the under- string in can be uniquely expressed as a sum of coset rep-
lying graph—possibly consisting of super states of varying de- resentatives
grees of coarseness and possibly differing from node to node
at any level—uniquely determines a super code containing .
Clearly, any codeword in can be derived from this super-state
grammar; given the codeword , there exists a super-state
derivation tree that corresponds (by the definition of super-state where and . We choose our coset
productions) to the codeword’s original state derivation tree and representatives for the quotient group RM RM
has the bits of as its terminal assignments. Thus, each such according to the scheme
super-state grammar generates a super code containing .
The CTFDP algorithm proceeds by progressively refining
the super state grammar. Given the solution of the previous where is the -bit binary representation of the integer
dynamic programming (DP) problem—an optimal derivation and is the generator matrix of -fold
tree corresponding to a minimum-cost super codeword, we Kronnecker products of the form of weight
determine whether the optimal derivation tree contains any ordered lexicographically by label —largest first.
non-singlet super states. If so, we refine these super states, Note that the generators and serve as a
566 IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 47, NO. 2, FRBRUARY 2001

basis for as do their corresponding -fold Kronnecker prod- to be the code generated by the start symbol and
ucts for . For a complete discussion of set and group parti-
tions of (upon which much of this discussion is based) see productions
[21].
One can verify by induction [36] that these symbols have a
familiar hierarchical structure. In fact, the iterative rule for con-
structing symbols is

where and . But by an

application of the additivity rule for this linear representation of
with terminals and . Note
that the number of productions is now
(adopting the convention that for or ).
A careful examination of its grammar (see [36]) reveals
the grammatical hierarchy of Section IV-A emerges
RM to be a linear
subcode of RM with symbols (or states) at
level . Of course, if we set the unthinned grammar
for RM is recovered.
The decoding complexity of thinned codes is controlled by
Notice that the combinatorial function has been replaced by the limiting parameter . By definition, the number of produc-
the alternative mod where rep- tions is bounded above by . Surprisingly, this condition in-
resents the bitwise exclusive OR operator. duces a similar restriction on the state-space cardinality. Thus,
Thus, we see that the derivation of the RM grammar is simply the number of operations required to decode the thinned RM
a matter of choosing the appropriate twisting. Instead of we in- code RM is loosely upper-bounded by (
troduce a “canonical” twisting which corresponds to Forney’s sites, states per site, productions per state), a large
“iterated group squaring construction” and essentially serves to but fixed multiple of the code length.
linearize the grammatical hierarchy. To distinguish the two de- We conclude our discussion of thinning by observing
velopments, we label the RM symbols with the letter (instead that one could as easily have thinned a general context-free
of ), and everywhere replace by . Then for the RM grammar code—simply by making the appropriate substitutions ;
1) equation (6) (with in place of and in place of ) such codes would exhibit all of the above features except
constructs (inductively) , a linear representa- linearity.
tion of ;
APPENDIX C
2) equation (8) (with in place of ) produces a collection
of distance codes , representing STATE-SPACE PARTITIONS
To implement a CTFDP version of a given DP problem, one
RM must first partition the problem’s state spaces into hierarchies
of super states. In this appendix, we present one particularly
effective choice of super states for thinned RM codes.
We begin by introducing a more compact notation for the
and its cosets in ; and thinned RM grammar. The set of allowed symbols at level for
3) the productions in (10) (with in place of and in the code RM can be re-expressed as
place of ) define a grammatical realization of the codes
.

APPENDIX B where the integer label is derived from the indexes

THINNING of by
For the class of context-free codes constructed in Section concatenation of binary representations. Specifically, let
IV-A and generated by RM-like grammars, one can bound the be the -bit binary expansion of , and let
size of the state space without disrupting the coherence of the re- be the concatenation of these expansions. Then
sulting context-free grammar. Suppose, for example, we impose is the integer with binary expansion ( ).
an upper bound of on the number of allowed productions in The label denotes the corresponding state of an allowed
an intractably large context-free grammar of this type. The yield symbol. The number of states at level is (see [36])
of this reduced grammar is a subcode of the original one—a whereas the number of productions allowed at level is
thinned code. We define the thinned RM code RM .
GEMAN AND KOCHANEK: PROGRAMMING AND GRAPHICAL REPRESENTATION OF ERROR-CORRECTING CODES 567

Among the advantages of this scheme is that single and mul- is an allowed super-state produc-
tiple productions can be jointly expressed as tion if and only if shares a (suitably sized—see [36]) binary
prefix with and shares a binary prefix with
for some ; for if this condition is met, there
is an allowed state production contained within the postulated
super-state production. We are thus able to compute super-state
where the auxiliary integers and are defined by the productions by inspection!
correspondences
REFERENCES
[1] R. Bellman, Dynamic Programming. Princeton, NJ: Princeton Univ.
Press, 1957.
and [2] L. R. Bahl, F. Jelinek, and R. L. Mercer, “A maximum likelihood ap-
proach to continuous speech recognition,” IEEE Trans. Pattern Anal.
Machine Intell., vol. PAMI-5, pp. 179–190, 1983.
[3] C. Cannings, E. A. Thompson, and H. H. Skolnick, “Probability func-
and the binary operator is the same as itself—bitwise tions on complex pedigrees,” Adv. Appl. Probab., vol. 10, pp. 26–61,
exclusive-OR, introduced for notational convenience. In 1978.
[4] J. E. Hopcroft and J. D. Ullman, Introduction to Automata Theory, Lan-
these expressions, is the -bit binary expansion of guages, and Computation. Reading, MA: Addison-Wesley, 1979.
and is the -bit expansion of [5] A. J. Viterbi, “Error bounds for convolutional codes and an asymptoti-
(or if ). Note that for single cally optimum decoding algorithm,” IEEE Trans. Inform. Theory, vol.
IT-13, pp. 260–269, Apr. 1967.
productions (i.e., ), the strings and have [6] N. Wiberg, H.-A. Loeliger, and R. Kötter, “Codes and iterative decoding
length and can, therefore, be ignored; however, for on general graphs,” European Trans. Telecommun., vol. 6, pp. 513–525,
1995.
multiple productions they cannot be ignored. [7] N. Wiberg, “Codes and decoding on general graphs,” Ph.D. dissertation,
A further distinction of this scheme is that the multitude of Dept. Elec. Eng., Linköping Univ., Linköping, Sweden, 1996.
productions for the RM code RM can be readily [8] G. D. Forney Jr., “On iterative decoding and the two-way algorithm,” in
Proc. Intl. Symp. Turbo Codes and Related Topics, Brest, France, 1997.
computed from a comparatively small set of stored inte- [9] F. Kschischang, B. Frey, and H.-A. Loeliger, “Factor graphs and the
gers—the parameters and the set sum-product algorithm,” Department of Electrical and Computer Engi-
of auxiliary ’s and ’s, one pair of integers for neering, University of Toronto, Tech. Rep., 1998.
[10] R. Dechter, “Bucket elimination: A unifying framework for proba-
each state at each level . But of far greater consequence is bilistic inference,” in Learning in Graphical Models, M. I. Jordan,
the foundation we have established for constructing a simple Ed. Cambridge, MA: MIT Press, 1999, pp. 75–104.
system of state-space partitions for thinned RM grammars. [11] J. Pearl, Probabilistic Reasoning In Intelligent Systems: Networks of
Plausible Inference. San Francisco, CA: Kaufmann, 1988.
The basic partition at level is constructed as a succes- [12] S. L. Lauritzen, Graphical Models. Oxford, U.K.: Oxford Univ. Press,
sion of binary refinements of the set of states. There are 1996.
coarsenings, , each with [13] B. D. Ripley, Pattern Recognition and Neural Networks. Cambridge,
U.K.: Cambridge Univ. Press, 1996.
super states denoted by the pair at [14] R. G. Gallager, “Low-density parity-check codes,” Ph.D. dissertation.
coarseness . Specifically, we define the coarsened symbols [15] R. M. Tanner, “A recursive approach to low complexity codes,” IEEE
Trans. Inform. Theory, vol. IT-27, pp. 533–547, Sept. 1981.
[16] L. R. Bahl, J. Cocke, F. Jelinek, and J. Raviv, “Optimal decoding of linear
codes for minimizing symbol error rate,” IEEE Trans. Inform. Theory,
vol. IT-20, pp. 284–287, Mar. 1974.
for . [17] C. Berrou, A. Glavieux, and P. Thitimajshima, “Near Shannon limit
error-correcting coding and decoding: Turbo-codes,” in Proc. Int. Conf.
The extraordinary feature of our choice of super states is that Communications (ICC’93), Geneva, Switzerland, 1993, pp. 1064–1070.
they inherit the underlying structure of the RM grammar. In fact, [18] D. J. C. MacKay, “Good error-correcting codes based on very sparse
these coarsened symbols obey the recursive set relation (proved matrices,” IEEE Trans. Inform. Theory, vol. 45, pp. 399–431, Mar. 1999.
[19] F. Kschischang and B. Frey, “Iterative decoding of compound codes
in [36]) by probability propagation in graphical models,” IEEE J. Select. Areas
Commun., vol. 16, pp. 219–230, Feb. 1998.
[20] L. N. Kanal and A. R. K. Sastry, “Models for channels with memory
(12) and their applications to error control,” in Proc. IEEE, vol. 66, 1978, pp.
724–744.
[21] G. D. Forney Jr., “Coset codes—Part II: Binary lattices and related
codes,” IEEE Trans. Inform. Theory, vol. 34, pp. 1152–1187, Sept.
where and are simple functions of the auxiliary integers 1988.
and , respectively, while and depend only on the level [22] C. S. Raphael, “Coarse-to-fine dynamic programming,” IEEE Trans.
and coarseness . Although this coarsening scheme is some- Pattern Anal. Machine Intell., to be published.
[23] L. R. Rabiner, “A tutorial on hidden Markov models and selected appli-
what cumbersome to express mathematically, its computational cations in speech recognition,” Proc. IEEE, vol. 77, pp. 257–286, 1989.
implementation is straightforward and facilitates the remarkably [24] Y. Amit and A. Kong, “Graphical templates for model registration,”
fast CTFDP decoding algorithm presented in Section V-C. IEEE Trans. Pattern Anal. Machine Intell., vol. 18, pp. 225–236, 1996.
[25] B. J. Frey, Graphical Models for Machine Learning and Digital Com-
For example, since the hierarchy of coarsened symbols munication. Cambridge, MA: MIT Press, 1998.
retains the underlying RM structure, the run-time compu- [26] J. Hammersley and P. Clifford, “Markov fields on finite graphs and lat-
tation of super-state productions required by the CTFDP tices,” Univ. California, Berkeley, Tech. Rep., 1968.
[27] G. Winkler, Image Analysis, Random Fields, and Dynamic Monte Carlo
procedure is trivial. If , , and are Methods: A Mathematical Introduction. New York: Springer–Verlag,
super states at the respective sites , , and , then 1995.
568 IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 47, NO. 2, FRBRUARY 2001

[28] G. Kallianpur, Stochastic Filtering Theory. New York: Springer- [34] W. C. Gore, “Further results on product codes,” IEEE Trans. In-
Verlag, 1980. [Link], vol. IT-16, pp. 446–451, July 1970.
[29] H. Künsch, S. Geman, and A. Kehagias, “Hidden Markov random [35] R. Kindermann and J. Snell, Markov Random Fields and Their Applica-
fields,” Ann. Appl. Probab., vol. 5, pp. 577–602, 1995. tions. Providence, RI: Amer. Math. Soc., 1980.
[30] S. Arnborg, D. G. Corneil, and A. Proskurowski, “Complexity of finding [36] K. Kochanek, “Dynamic programming algorithms for maximum like-
embeddings in a k-tree,” SIAM J. Alg. Discr. Meth., vol. 8, pp. 277–284, lihood decoding,” Ph.D. dissertation, Div. Appl. Math., Brown Univ.,
1987. Providence, RI, 1998.
[31] G. D. Forney Jr., “The Viterbi algorithm,” Proc. IEEE, vol. 61, pp. [37] J. Pearl, Heuristics: Intelligent Search Strategies for Computer Problem
268–278, 1973. Solving. Reading, MA: Addison-Wesley, 1984.
[32] P. W. Fieguth and A. S. Willsky, “Fractal estimation using models [38] A. Kam and G. Kopec, “Document image decoding by heuristic search,”
on multiscale trees,” IEEE Trans. Signal Processing, vol. 44, pp. IEEE Trans. Pattern Anal. Machine Intell., vol. 18, pp. 945–950, Sept.
1297–1300, May 1996. 1996.
[33] M. Meila, “Learning with mixtures of trees,” Ph.D. dissertation, MIT, [39] , “The iterated complete path algorithm,” Xerox Palo Alto Research
Cambridge, MA, 1999. Center, Tech. Rep., 1995.