Page Cover.html

[Cover] [Contents] [Index] Next Section

[Cover] [Contents] [Index] Next Section


Page 1
1
Introduction
In this book, I explore algorithms for pattern classification, unsupervised learning, data compression, and channel coding. At first, it may be hard to imagine how these different research areas can be brought into focus under a single theme that is both novel and of practical value. My hope is to convince the reader that these problems can be attacked in an interesting and fruitful manner using a recently developed class of algorithms that make use of probabilistic structure. These algorithms take advantage of a graphical description of the dependencies between random variables in order to compute, or estimate, probabilities derived from a joint distribution. As simple and well-known examples, the forward-backward algorithm [Baum and Petrie 1966] and the Viterbi algorithm [Forney 1973] make use of a chain-like Markovian relationship between random variables.
The roots of probabilistic structure reach far back to the beginning of the 20th century. In 1921, Sewall Wright [1921] developed "path analysis" as a means to study statistical relationships in biological data. Few new developments were made until the 1960's when statisticians began using graphs to describe restrictions in statistical models called "log-linear models" [Vorobev 1962; Goodman 1970]. In 1963, the idea of hierarchical probabilistic structure briefly reared its head in the engineering research community when Gallager [1963] invented an error-correcting decoding algorithm based on a graphical description of the probabilistic relationships between variables involved in channel coding. Most likely because of the primitive computers available at the time, his algorithm was quickly overlooked by his peers, only to be rediscovered nearly 35 years later independently by at least three research groups, and to be shown to yield unprecedented performance in error-correcting coding applications [Berrou, Glavieux and Thitimajshima 1993; Wiberg, Loeliger and Kotter 1995; MacKay and Neal 1995]. A simpler chain-type Markovian graphical structure later became popular and very useful in the engineering community, largely due to a good tutorial paper by Forney [1973], in which the notion of a "trellis'' was introduced. Probabilistic structure has been extensively developed in the artificial intelligence literature, with applications ranging from taxonomic hierarchies [Woods 1975; Schubert 1976] to medical diagnosis [Spiegelhalter 1990]. In the late 1980's, Pearl [1986] and Lauritzen and Spiegelhalter [1988] independently published a general algorithm for computing probabilities based on a graphical representation of probabilistic structure. This algorithm is practical and exact for only a special type of probabilistic


Page 2
structure that is characterized by a relatively small "state space". Over the last decade, there has also been a tremendous increase in interest in estimating the parameters of models with fixed graphical structure. In the mid 1980's, Hinton and Sejnowski [1986] introduced a maximum likelihood algorithm for learning the parameters of a graph-based log-linear model called a "Markov random field". More recently, approximate algorithms for general models based on directed graphical models called ''Bayesian networks" have been introduced. These include Markov chain Monte Carlo methods [Pearl 1987; Neal 1992], "Helmholtz machines'' [Hinton et al. 1995; Dayan et al. 1995], and variational techniques [Ghahramani and Jordan 1998; Saul, Jaakkola and Jordan 1996; Jaakkola, Saul and Jordan 1996].
1.1 A probabilistic perspective
Offhand, it is not obvious that sophisticated probability models are needed to solve problems in machine learning and digital communication. Given a segment of speech, a classifier outputs a decision, say, as to whether or not the speaker has security clearance. It appears there are no random variables in this model. The classifier may also output a measure of reliability regarding the decision it makes. In this case, it appears there is just one binary random variable that captures the variability in the decision. The mean of this Bernoulli random variable must somehow be related to the input, and this task can be viewed as some sort of function approximation. Similarly, in response to a series of sensory inputs, an animal makes a single motor action. Given a highly redundant image, a data compression algorithm usually produces a unique sequence of codeword symbols. Given the output of a noisy telephone line, a channel decoder (telephone modem) makes a deterministic decision about the contents of the transmitted data file.
The above modelling approaches either require only very low-dimensional probability models or do not use random variables at all, and so they are clearly not representing the true underlying structure in each problem. For example, in reality each speaker has a unique glottis that interacts in a random way with a unique shape of vocal tract and a unique random style of articulation to produce a speech segment. It seems that a fruitful approach to speaker identification would involve representing these random variables and the probabilistic relationships between them. In the following four sections, I give some general ideas about what is missing in the above approaches to pattern classification, unsupervised learning, data compression, and channel coding.


Page 3
1.1.1 Pattern classification
A soft-decision classifier estimates the probability that a given pattern v belongs to each class j Î {0, . . . , J - 1}. That is, the classifier estimates
where the subscript "r" in Pr indicates a true (real) probability (as opposed to one produced by a model). If these probabilities can be accurately estimated, Bayes' decision theory tells that a minimum error rate is achieved by choosing the class j that maximizes Pr( j½v) [Chow 1957; Duda and Hart 1973].
We could use a logistic regression model to estimate Pr( j½v). For example, regression with a flexible neural network has been successfully used to classify individual digits extracted from handwritten United States ZIP codes [Le Cun et al. 1989]. Alternatively, support vector techniques can be used to obtain conditional probability estimates [Vapnik 1998]. However, these approaches ignore the causal structure in the physical production of the input patterns.
In order to faithfully capture the actual physical process that produces each digit, we first ought to specify an a priori distribution P( j) over the digit classes j Î {0, . . . , 9}—maybe some digits are more common than others. Next, for a given class of digit j, we expect there to be a distribution P(h½j) over a set of digit attributes h. These attributes are called "hidden variables", because they are not part of the classifier inputs or outputs. Each element of h might specify the presence or absence of a particular line segment or flourish. Given a set of features h, we expect there to be a distribution P(v½h) over possible images—this distribution models the way in which features combine to make an image, as well as noise such as ink spots. The joint distribution given by this model of the real world can be written
and the distribution over classes given a pattern can be obtained using Bayes' rule with h marginalized out:
So, it appears that to properly model the structure of the problem, we need a more sophisticated probabilistic description than (1.1). In general, a correct model of this sort, where P( j½v) » Pr( j½v), will perform optimally in terms of classification error.


Page 4
1.1.2 Unsupervised learning
The goal of unsupervised learning is for a machine to process a set of training data and then extract structure which cognitive scientists believe to be relevant to perception and conception. For example, given a training set containing examples of handwritten digits, we would like our machine to learn to represent digit class without being told what the digit classes are. When the machine is given a new test image of, say, a 3, we hope that somewhere within the machine there will be a flashing light saying "I think this is a 3!"
One powerful approach to developing unsupervised learning algorithms is to suppose that in addition to the input sensory variables v there are hidden "concept" variables h, and that all the variables are linked together by a parameterized probability model P(v, h½q ), with parameters q . This probability model has a form that we believe to be appropriate for the task of cognition. For example, we may believe that given the hidden variables, the sensory inputs are independent:
The model is fit to a set of training examples [image] using maximum likelihood parameter estimation:
where the marginal probability of the sensory inputs is computed from
We hope that the maximum likelihood model will represent the true physical process that produced the patterns; i.e., that the hidden variables will come to represent meaningful concepts. For example, in the case of unsupervised learning of digits, when a new pattern v is presented to the trained machine, we hope that in the posterior distribution over the hidden variables,
high probability will be given to a representation that correctly identifies the digit.


Page 5
Figure 1.1 An ideal distortion-free digital communication system.
1.1.3 Data compression
Figure 1.1 shows an ideal distortion-free digital communication system. We assume the data source produces patterns that are independent and identically drawn (i.i.d.) with true information content (entropy) [image] bits per pattern. Each source pattern v is fed into an ideal source encoder, which on average produces a block of [image] information bits that look perfectly random. An ideal channel encoder takes in a block of these random-looking information bits (with a length that is usually different from [image]) and produces a longer block of binary signals. The redundant signals in the longer block identify relationships between the information bits, so that an information bit that is received in error can be corrected using the "redundant" signals. To produce a block of N binary signals, the channel encoder takes in a block of RN bits, where R is called the "rate".
An ideal channel decoder takes in a block of N noise-corrupted signals, which are in general now real-valued instead of binary (e.g., if the channel noise is Gaussian). The decoder produces RN information bits that are error-free. An ideal source decoder on average takes in [image] error-free information bits and produces a distortion-free source pattern. Notice that memory may be needed at each stage in this block diagram since, e.g., several source patterns


Page 6
may be needed to obtain enough information bits for the channel encoder to produce a block of signals.
It is important to recognize that digital communication can be broken down into two stages [Cover and Thomas 1991 ]: compressing data to obtain random-looking information bits and then communicating this compressed data in an error-free fashion.
A source code maps each input pattern v to a source codeword u, such that for each valid u there is a unique pattern. I will consider sources where i.i.d. patterns are generated from Pr(v). The purpose of noiseless source coding, or lossless data compression, is to losslessly represent the source patterns by codewords, so that the expected codeword length is as small as possible. Shannon's noiseless source coding theorem [Shannon 1948] states that the average codeword length per source pattern cannot be less than the source entropy:
where [image] is the length of the codeword for v in bits, and [image] is the entropy of the source in bits:
Shannon's theorem also states that if the input patterns are processed in blocks, there exists a source code with an average codeword length per input pattern that is arbitrarily close to
Arithmetic coding [Rissanen and Langdon 1976; Witten, Neal and Cleary 1987] is a practical algorithm for producing near-optimal codewords when the source distribution Pr(v) is known. Sometimes, e.g., if v is binary-valued, these probabilities can be easily estimated from the source. Often, however, the distribution is too complex, and so a more sophisticated parametric model or flexible model must be used to estimate the probabilities. For example, consider a high-dimensional binary image v that is produced by the physical process described in Section 1.1.1, so that
The probabilities used by the arithmetic encoder are obtained by marginalizing out j and h:
Thus, probabilistic descriptions can also be very useful for source coding.


Page 7
1.1.4 Channel coding
A channel code maps a vector of information bits u to a vector of codeword signals x. This mapping adds redundancy to u in order to protect the block against channel noise. (As a simple example, the codeword might consist of three repetitions of the information vector.) After x is transmitted across the channel, the decoder receives a noise-corrupted version y and produces an estimate of the information block û We say that a block error or a word error has occurred if û ¹ u. In its simplest form, Shannon's channel coding theorem [Shannon 1948] states that for any given channel, there exists1 a channel code that can achieve an arbitrarily low probability of block error when the signal-to-noise ratio is greater than a channel-dependent threshold called the Shannon limit. Roughly speaking, the codewords are kept far apart in codeword symbol space, so that when a moderately noise-corrupted codeword is received, it is still possible to determine with high probability which codeword was transmitted.
From a probabilistic perspective, the decoder can minimize the word error rate by choosing an estimate û that maximizes Pr(û|y), or minimize the symbol error rate by choosing an estimate û that maximizes Pk Pr(ûk|y). A probabilistic model can be constructed by examining the encoding process and the channel. We first specify a (usually uniform) distribution for the information blocks, P(u). Often, the encoder uses a set of state variables, s, in order to produce the codeword. These variables are determined from the information block using a distribution P(s|u)—although this relationship is usually deterministic, this probabilistic description will come in handy later on when we study probabilistic decoding. The transmitted codeword is determined from the information block and state variables by P(x|u,s). Finally, the real-valued channel outputs are related to the transmitted codeword by a probability density function p(y|x) that models the channel. The joint distribution given by the model is
and the distribution over information symbol uk given the channel output can be obtained by marginalizing out s, x and u j for all j ¹ k and using Bayes’ rule:
1 Shannon was quite the tease. He proved the code exists, but did not show us a practical way to encode or decode it.


Page 8
Although this probabilistic formulation may seem strange compared to many of the strongly algebraic traditional approaches, it is this formulation that I view as the foundation of the recently proposed high-performance turbocodes [Berrou and Glavieux 1996].
1.1.5 Probabilistic inference
As presented above, pattern classification, unsupervised learning, data compression, and channel coding are all similar in that some type of marginal (and possibly conditioned) distribution is sought for a given joint distribution. Consider a set of random variables z = (z1, z2, . . . , z N) that covary according to a joint distribution P (z1, z2, . . ., z N ). For any two subsets of variables z1 Í z and z2 Í z, I will refer to the computation or estimation of P(z1|z2), or a decision based on P(z1|z2), as probabilistic inference.
Examples of probabilistic inference include the computation of the class probabilities for pattern classification (1.3), the computation of the distribution over concepts in unsupervised learning (1.7), the computation of the input probability for data compression (1.1 1), and the information symbol decisions based on the information symbol probabilities for channel coding (1.13). Notice that in these different cases of probabilistic inference, the joint distributions can be decomposed in different ways. In fact, if we decompose the joint distributions at the level of individual variables instead of vector variables, we can envision a wide variety of rich structures. In the next section, I describe graphical models, which can be used to describe this structure.
1.2 Graphical models: Factor graphs, Markov random fields and Bayesian belief networks
Often, the joint distribution associated with a probabilistic inference problem can be decomposed into locally interacting factors. For example, the joint distributions involved in the applications of Bayes' rule in (1.3), (1.7), (1.1 1), and (1.13) can be expressed in the forms given in (1.2), (1.4), (1.10), and (1.12). By taking advantage of such probabilistic structure, we can design inference algorithms that are more efficient than the blind application of Bayes' rule.
Probabilistic structure can be characterized by a set of conditional independence relationships. (This structural description does not fix the values of the


Page 9
probabilities.) For example, in the case of channel coding, we can use the chain rule of probability to write out the joint distribution:
The probability density function (the last factor) can be simplified, since the received vector y is conditionally independent of the information vector u and the state vector s, given the transmitted codeword x:
By substituting this conditional independency relationship into (1.14), we obtain the more structured form of the joint distribution given in (1.12).
The general idea is to express the joint distribution as a product of factors, where each factor depends on a subset of the random variables. In the simplest case, each factor depends on a single random variable, making marginalization easy. Most distributions that describe practical problems cannot be broken up in this way, and the subsets overlap. Within this richer set of models, some structures lead to highly efficient exact algorithms (e.g., the forward-backward algorithm for a chain-type structure). Other structures are not tractable and lead to approximate algorithms.
It turns out that graph theory provides a succinct way to represent probabilistic structure. A graphical representation for probabilistic structure, along with functions that can be used to derive the joint distribution, is called a graphical model. Examples of graphical models include Markov random fields, Bayesian networks (a.k.a., belief networks, causal networks, or influence diagrams), chain graphs [Lauritzen and Wermuth 1989], and factor graphs. Not only does the graphical representation concisely capture probabilistic structure, but it forms a framework for computing useful probabilities.
1.2.1 Factor graphs
A factor graph [Frey et al. 1998] explicitly indicates how a joint function of many variables factors into a product of functions of smaller sets of variables. Here, we will consider global functions that are probability distributions, although in general factor graphs can be used to describe the decomposition of any function that evaluates to a semi-ring [Wiberg, Loeliger and Kötter 1995; Aji and McEliece 1997; Frey et al. 1998]. Consider a distribution over a set of variables z that can be written as follows:


Page 10
Figure 1.2
The factor graph for the
global probability distribution,
P(u, s, x, y)=P(u)P(s|u)P(x|u,s) p(y|x).

where zj Í z is the subset of variables that the jth local function ¦ j (·) depends on and a is a normalization constant. The local functions may be marginal distributions or conditional distributions (as in (1.12)), or they may be "potentials" like the ones used for Markov random fields (see Section 1.2.2). A factor graph is a bipartite graph (meaning there are two types of vertex and edges may only connect vertices of differing type), with one set of vertices corresponding to the variables in z and another set of vertices corresponding to the local functions, ¦1(·), . . . , ¦J(·). Each local function is connected to the variables on which it depends. For example, Figure 1.2 shows the factor graph for the channel coding model given by (1.12). Notice that the function P(u) can be included with P(s|u) or P(x|u,s) if we wish.

A factor graph may have directed edges in order to indicate conditional dependence as is done in Bayesian networks (see Section 1.2.3). It turns out that in terms of expressing factorization of the global distribution, factor graphs are more general than Markov random fields and Bayesian networks.
1.2.2 Markov random fields
A Markov random field (MRF) (see Kinderman and Snell [1980] for a tutorial) is an undirected graph with vertices corresponding to variables (there aren't any local function vertices as there are in factor graphs). Let ni be the neighbors of variable z i in the MRF. Then, the MRF defines a local Markov property with respect to the global distribution:
where z\{z i} are all the variables in z, except z i. That is, given its neighbors, each variable is independent of all other variables. The global distribution for an MRF can be written in terms of a product of clique potentials, as long as the distribution is positive. A clique is a fully-connected subgraph that cannot


Page 11
Figure 1.3
(a) The Markov random field (MRF) for the global
probability distribution, P(u, s, x, y) =P(u)P(s|u)P(x|u,s) p(y|x).
(b) The MRF drawn as a factor graph.
remain fully connected if more variables are included. For example, in the MRF shown in Figure 1.3a, there are two cliques: z1 = (u,s,x) and z2 = (x,y) If there are C cliques and zj are the variables in clique j then the distribution can be written
where [image] is the potential for clique j and a is a normalization constant.
Notice that the factorization of P(·) implied by any MRF can also be represented by a factor graph (compare (1.16) and (1.18)). We simply create one local function vertex for each clique and connect it to the clique variables. For example, Figure 1.3b shows the factor graph for the MRF shown in Figure 1.3a.
In contrast, not all factorizations that can be represented by factor graphs can be represented by MRFs. Factor graphs are more explicit about factorization than are MRFs. For example, both the factor graph shown in Figure 1.2 and the factor graph shown in Figure 1.3b correspond to the MRF shown in Figure 1.3a. In fact, it is not possible to construct an MRF that explicitly gives the factorization P(u, s, x, y) = P(u)P(s|u)P(x|u, s) p(y|x).
1.2.3 Bayesian networks
Unlike MRFs, Bayesian networks [Pearl 1988] are specified in terms of directed acyclic graphs, in which all edges are directed and in which there are no closed paths when edge directions are followed. A Bayesian network for a set of random variables z = (z1, z2,. . . ,z N ) consists of a directed acyclic graph with one vertex for each variable, and a set of probability functions P(z k|ak), k = 1, . . . , N, where the parents ak of z k are the variables that have directed edges connecting to z k. If z k has no parents, then ak = Æ. For now, we can think of each function P(z k|ak) as an exhaustive list of probabilities corre-


Page 12
Figure 1.4
(a) The Bayesian network for the global probability
distribution, P(u, s, x, y) =P(u)P(s|u)P(x|u,s) p(y|x)
The Bayesian network drawn as a factor graph.
sponding to the possible configurations of z k and ak. The global probability distribution is written
For example, Figure 1.4a shows the Bayesian network for the coding model discussed above.
A Bayesian network can be drawn as a factor graph by introducing one function vertex for each conditional probability function (effectively producing a moral graph [Lauritzen and Spiegelhalter 1988]). The edges of the factor graph may be directed to indicate that the local functions are conditional probabilities, as shown in Figure 1.4b.
As with MRFs, we can actually define a Bayesian network in terms of its independence properties and then derive (1.19). Several definitions will be useful for making this connection. The children ck of z k are the variables that have directed edges connecting from z k. The descendents dk of z k are its children, its children's children, etc. The nondescendents nk of z k are the variables in {z1, z2, . . . , z k -1, . . . ,z N} that are not in dk i.e., nk = z \ (dk È { z k }). Note that nk ¹ z\ dk, since zk is not included in the nondescendents. From these definitions, it follows that ak Í nk. Figure 1.5 shows an example of a Bayesian network, along with the parents, children, descendents and nondescendents of variable z5.
The defining conditional independence property of a Bayesian network is that given the parents of any variable z k, the distribution over z k will not change if any combination of the nondescendents of z k are also given:


Page 13
Figure 1.5
(a) An example of a Bayesian network. (b) The parents a5 of z5 are shown
by a solid loop; the children c5 are shown by a sparse dotted loop; the
descendents d5 are shown by a dense dotted loop; the nondescendents
n5 are shown by a dashed loop.
In other words, z k is conditionally independent of any combination of its nondescendents, given its parents. To take the family hierarchy (not necessarily a tree) analogy further, given the genetic code of Susan's parents, determining the genes of her siblings, her grandparents, her grandparents' children, her children's other parents or any combination of the above does not influence our prediction of Susan's genetic make-up. This is not true for descendents. For example, determining the genes of Susan's grandchildren does influence our prediction of her genetic make-up, even though determining the genes of those parents of Susan's grandchildren who are not Susan's children does not (notice that the latter are nondescendents).2
The expression for the global distribution (1.19) is derived from (1.20) in the following way. Since a Bayesian network contains no directed cycles, it is always possible to choose an ancestral ordering zp (1), zp (2), . . . , zp (N ), where p(×) is a permutation map, so that the descendents of each variable come later in the ordering: dp (k) Í {zp (k+1), . . . , zp (N )}. For example, the variables in the network shown in Figure 1.5 were assigned so that z1, z2, . . . , z13 is an ancestral ordering. Using the general chain rule of probability applied to the ancestral ordering, we have
2 Interestingly, if we have previously determined the genes of Susan's grandchildren, then determining the genes of those parents of Susan's grandchildren who are not Susan's children does influence our prediction of Susan's genetic make-up. See Section 1.2.5 for more details.


Page 14
Figure 1.6
An example of a chain type Bayesian
network, or a Markov chain.
From the definition of the ancestral ordering, it follows that the set of variables that precede zp (k) is a subset of its nondescendents and that the parents of zp (k) are a subset of the variables that precede zp (k):
For this reason, { zp (1), . . . , zp (k-1)} = ap (k)È({ zp (1), . . . , zp (k-1)} \ ap (k)), and taking w = { zp (1), . . . , zp (k-1)} \ ap (k) in (1.20), we have
Inserting this result into (1.21), we obtain the form given in (1.19).
If the probability functions for a Bayesian network are not specified, the network is meant to represent all distributions that can be written in the form given in (1.19). For the network with ancestral ordering z1, z2, . . . , z13 shown in Figure 1.5, (1.19) gives a joint distribution
This product could have been written in any order, but using an ancestral ordering helps clarify the dependencies. In this equation, each variable z k is conditioned on variables whose distributions appear to the left of the distribution for z k. Note that a Bayesian network may have more than one ancestral ordering. In this case, z10, z6, z3, z7, z1, z2, z4, z5, z8, z11, z9, z13, z12 is also an ancestral ordering. An interesting special case of a Bayesian network is the chain-type network shown in Figure 1.6, also known as a first-order Markov


Page 15
chain. Applying (1.19) to this network, we obtain
This type of structure is frequently used to model time series data, where it is often assumed that the next state of a physical system depends only on the previous state. Comparing this network to the more complex networks that appear later in this book, the Bayesian network can be thought of as a generalization of the Markov chain.
1.2.4 Ancestral simulation in Bayesian network s
It is often practically very difficult to simulate vectors z that are distributed according to P(z). For example, in a Markov random field it can take a very long time for a stochastic simulation to reach equilibrium. However, if the joint distribution can be described by a Bayesian network, and if a value for each z k can be drawn from its conditional probability P(z k|ak) in a practical manner, then the ancestral ordering can be used to draw an entire vector. Starting with k = 1 we draw zp (1) from P(zp (1)) We continue to draw zp (k) from P(zp (k)|a p (k)) for k = 2, . . . , N until an entire vector z has been drawn. In this way, the probabilistic structure implied by the graph allows us to decompose the simulation problem into local pieces.
1.2.5 Dependency separation in Bayesian networks
In Section 1.2.3, we saw that a Bayesian network implies that variable z k is conditionally independent of any subset of its nondescendents, given its parents. This is expressed mathematically by (1.20). A convenient way to describe this scenario is to say that '' z k is dependency-separated from any combination of its nondescendents by its parents''.
Consider the uncountable set of distributions [image] that can be described by a given Bayesian network. In general, I will say that zA is dependency-separated from zB by zS ("S" for separation), if and only if
for all P Î [image] . (See Pearl [1988] for an extensive discussion of dependency-separation.) Notice that dependency-separation is symmetric with respect to zA and zB. The type of dependency-separation that can we used to define a Bayesian network in Section 1.2.3 is special, in that zS was the set of parents of the single variable zA and zBwas a subset of the nondescendents of zA.


Page 16
Figure 1.7
The three ways in which a path may be blocked. If all paths from all
variables in zA to all variables in zB are blocked, then zA and zB are
dependency-separated by zs.
It is possible to ascertain dependency separation in the general case simply by inspecting the Bayesian network. If zA, zB, and zS are three disjoint subsets of z, then zS dependency-separates zA from zB if, in every path connecting any variable in zA to any variable in zB, there is at least one variable z k that satisfies one or more of the following three conditions:
1. z k acts as both a parent and a child in the path and z k Î zS (Figure 1.7a), or
2. z k acts as the parent of two variables in the path and zk Î zS (Figure 1.7b), or
3. z k acts as the child of two variables in the path and neither z k nor any of its descendents are in zS (Figure 1.7c).
(Note that the identification of a path does not depend on edge directions.) A path for which one of these conditions is met is said to be blocked. In order to ascertain dependency-separation, we need to consider only paths that do not intersect themselves, since those conditions that hold for any given path will also hold for the path when extra variables that form a loop are considered.
For example, in Figure 1.5a we have the following dependency separation relationships. z6 is dependency-separated from z10(by nothing), since the sole path z6, z9, z13, z10 is blocked by z13 in condition 3. What if z9 is observed? Then, z6 is still dependency-separated from z10, by condition 1 applied to the sole path z6, z9, z13, z10. In contrast, if only z13 is observed, then z6 is not dependency-separated from z10, since there exists a path z6, z9, z13, z10 for which none of the three conditions can be met. This means that once z13 is observed, z6 and z10 may become dependent. Note that there may exist a distribution in [image] where z6 and z10 are independent given z13 , but there exists at least one distribution in [image] where z6 and z10 are dependent given z13.
Here are some more complicated examples. z2 is dependency-separated from z9 by z5, since path z2, z5, z9 is blocked by z5 in condition 1, paths


Page 17
Figure 1.8
The LFSR for a systematic recursive convolutional
code that has a minimum distance of 10.
z2, z4, z8, z12, z9 and z2, z1, z4, z8, z12, z9, are blocked by z12, in condition 3, paths z2, z4, z8, z5, z9 and z2, zl, z4, z8, z5, z9 are blocked by z5 in condition 2 and by z8 in condition 3, path z2, z5, z8, z12, z9 is blocked by z5 in condition 1 and by z12 in condition 3. z2 is dependency-separated from {z3, z7} (by nothing), since the paths z2, z5, z3 and z2, z5, z3, z7 are blocked by z5 in condition 3. This means that in the absence of observations z2 and {z3, z7} are independent. z2 is not dependency-separated from {z3, z7} by z12, since there exists a path z2, z5, z3 for which none of the conditions can be met. Condition 3 applied to z5 fails, because z12 is a descendent of z5. This means that once z12 is observed, z2 and {z3, z7} may become dependent.
1.2.6 Example 1: Recursive convolutional codes and turbocodes
Recall from Section 1.1.4 that the purpose of channel coding is to communicate over a noisy channel in an error-free (or nearly error-free) fashion. To do this, we encode a given binary information vector u as a longer codeword vector x, which contains extra bits whose purpose is to "protect" the information from the channel noise. (An example is a repetition code, where each information bit is simply transmitted several times.) The codeword is converted to a physical form (e.g., radio waves) and then sent over a channel. A vector of noisy signals y is received at the output of the channel. Given y, the decoder must make a guess û at what the original u was.
One very popular class of channel codes (you probably have one of these in your telephone modem) can be described using Bayesian networks. The encoder for a recursive convolutional code is simply a linear feedback shift register (LFSR) that takes in information bits and generates codeword bits. See Lin and Costello [1983] for an extensive treatment of convolutional codes. Figure 1.8 shows the LFSR for the convolutional code that is described below. Each box represents a 1-bit memory element and D indicates a delay buffer.


Page 18
The discs represent addition modulo 2 (XOR). For this particular convolutional code, every second output is actually just a copy of the input bit. This type of code is called systematic. Notice that for each input bit, two output bits are produced, so this is a rate 1/2 code. If there are K information bits, then there will be N = 2K codeword bits. The device shown in Figure 1.8 is called a linear feedback shift register because the output sequence generated by the sum of two input sequences is equal to the sum of the two output sequences that are generated by the individual input sequences (where summation is modulo 2). The details of how to choose the feedback delay taps and the output taps in order to produce a good code can be found in Lin and Costello [1983; 1996]. However, the operation of an encoder of this type is quite simple. The LFSR is initialized so that all memory elements contain 0's. Then, the information bits xk are fed into the LFSR, producing codeword bits x k. Signals that represent the codeword bits are then transmitted over the channel. For example, on a twisted pair of wires, we might apply +1 volts if x k = 1 and -1 volts if x k = 0.
Figure 1.9a shows the Bayesian network for a recursive systematic convolutional code. Normally, the number of information bits K is much larger than 6 (typical numbers range from 100 to 100,000 bits). s k is the state of the LFSR at time k, extended to include the input bit (this makes the network simpler). To fully specify the Bayesian network, we must also provide the conditional distributions. Assuming the information bits are uniformly distributed,
Let S (s k-1, u k) be a function (determined from the LFSR) that maps the previous state and current input to the next state, and let g(s k) be a function that maps the state to the nonsystematic output bit. Then, the deterministic conditional distributions for the states and state outputs are
where d(a, b) = 1 if a = b and 0 otherwise. Assuming that the channel simply adds independent Gaussian noise with variance s 2 to the +1/-1 signals described above, the conditional distributions for the received channel output signals are


Page 19
Figure 1.9
(a) The Bayesian network for a recursive systematic convolutional code.
(b) The Bayesian network for an example of the recently proposed turbocode.
(c) The bit error rate (BER) performance for instances of these two codes.
Given an information vector u, encoding and channel transmission can be simulated by one sweep of ancestral simulation. For example, we first directly copy u0 into x0, which is then used to draw a noisy channel output value y0. Then, we use u0 to determine so, which is then used to determine x1, which is then used to draw a noisy channel output value y1. Then, we directly copy u1 into x2 and so on until the entire channel output vector y has been obtained.
The decoder sees only the vector y, and ideally would infer the most likely value of each information bit, i.e., determine for each k the u k that maximizes P(u k|y). In general such a probabilistic inference is very difficult, but if we take advantage of the graphical structure of the code it turns out it can be


Page 20
done quite easily. In fact, it is possible to compute P(u k|y) k = 0, . . . , K - 1 exactly using the forward-backward (a.k.a. BCJR) algorithm [Baum and Petrie 1966; Bahl et al. 1974], which is just a special case of the general probability propagation algorithm discussed in Section 2.1. Once the block is decoded, we can compare the decoded information bit values with the true ones to determine the number of bit errors made for the block transmission. If we simulate the transmission of many blocks, we can obtain an estimate of the bit error rate (BER).
This procedure was carried out for the convolutional code shown in Figure 1.8. This recursive systematic convolutional code was designed to maximize the smallest Hamming distance between all pairs of codewords [Viterbi and Omura 1979; Lin and Costello 1983] (dmin = 10). The information vector length was K = 5000 (giving a codeword length of N = 10000), and 5000 vectors (25 ´ 106 information bits in all) were transmitted for a fixed noise variance. It is common practice to give BER results as a function of the noise level measured by a signal-to-noise ratio E b/N0 in decibels. For any system with N = 2K and transmission power (variance) of unity, E b/N0 is related to s 2 by E b/N0 = -10 log10 s 2. Figure 1.9c shows the BER as a function of Eb/N0 for this recursive systematic convolutional code.3 Notice that as E b/N0 increases (s 2 decreases), the BER drops.
In the same figure, I also give the BER curve for a simple repetition code, where each information bit is transmitted twice, maintaining N = 2K. If the information bit is 0, a pair of -1's are sent; if the information bit is 1, a pair of +1's are sent. Each pair of received noisy signals is then averaged before a threshold of 0.0 is applied to detect the information bit. The curve on the far left shows Shannon's limit; for a given Eb/N0, it is impossible to communicate with a BER below this curve using binary signalling. (See Sections. 6.1.6 and 6.1.7 for a derivation of this curve.) So, systems of practical interest give performance points that lie between the Shannon limit curve and the curve for the repetition code. Performance points to the left of the Shannon limit are impossible, and performance points to the right of the curve for the repetition code are not of practical interest.
Recently, a code and decoding algorithm were discovered that give unprecedented BER performance. It turns out that the turbodecoding algorithm for these turbocodes [Berrou and Glavieux 1996] is just the probability propagation algorithm discussed in Section 2.1 applied to a code network like the one shown in Figure 1.9b [Frey and Kschischang 1997; Kschischang and Frey 1998; McEliece, MacKay and Cheng 1998]. This Bayesian network contains two recursive convolutional code networks that are connected to the infor-
3 Trellis termination [Lin and Costello 1983] was used to improve the performance of the code.


Page 21
mation bits in different ways. The information bits feed directly into one of the chains (s1), but feed into the second chain (s2) in a permuted order as shown. In order to produce the same number of codeword bits per codeword as would be produced by the recursive systematic convolutional encoder described above, every second output of each LFSR is alternately not transmitted (a procedure called puncturing).
Figure 1.9c shows the BER performance for a turbocode system with K = 65,536 and N = 131,072. 530 vectors (~ 35 ´ 106 information bits) were transmitted to determine the BER for each noise level. Each of the two LFSRs had 4 bits of memory and used identical feedback and output delay taps. All four delayed bits were fed back to the input of the LFSR. Only the bit entering the first delay element and the most-delayed bit were fed forward to the output. (This block length and these constituent LFSRs were proposed by Berrou, Glavieux and Thitimajshima [1993]). The decoding complexity per information bit for the turbocode was roughly twice that for the convolutional code described above. The information bit permuter was chosen at random. The turbocode system clearly outperforms the computationally comparable single convolutional code system. At a BER of 10-5, the turbocode system is tolerant to 3.3 dB more noise than the single convolutional code system, and is only 0.5 dB from the Shannon limit. Also shown on this graph is the performance of a concatenated Reed-Solomon convolutional code described in Lin and Costello [1983], which had been considered to be the best practical code until the proposal of turbocodes. The turbocode system is tolerant to 1.5 dB more noise than the concatenated system.
In Chapter 6, I explore some of the exciting new applications of Bayesian networks to channel coding problems, with a focus on using the probability propagation algorithm discussed in Section 2.1 for inference.
1.2.7 Parameterized Bayesian networks
It is sometimes convenient to represent the conditional distributions P(z k ½ak) in parametric form. That is, the distribution over z k given its parents ak is specified not by an exhaustive list of probability masses, but by a function of z k, ak, and a set of parameters q k. (The subscript k indicates that q k is a set of parameters associated with z k.) In this case, we write the conditional distribution as P(z k ½ak, q k). The total set of parameters is q = {q 1, . . . , q N }, and the parameterized joint distribution is expressed as P(z½q ). Such a parametric form can be useful in applications such as density estimation, pattern classification, and data compression, where the distribution P(z½q ) is to be estimated from a data set. In this case, the parametric form can act as a regularizer. Since the number of possible configurations of each z k and ak is usually quite large, we would need an extremely large data set to estimate all probabilities


Page 22
Figure 1.10
The logistic function g(x) = 1/(1 + exp[-x]).
accurately. Using the parametric form, however, we need only estimate each parameter. As described in Section 2.3, a parametric form is also useful when formulating variational inference algorithms.
A common parametric Bayesian network is the sigmoidal Bayesian network [Neal 1992; Jordan 1995; Saul, Jaakkola and Jordan 1996], whose random variables are all binary. The conditional probability function P(z k |ak, q k ) can be viewed as a regression model that is meant to predict z k from a set of attributes ak. A standard statistical method for predicting a binary-valued variable is logistic regression [McCullagh and Nelder 1983], in which the conditional probability for z k given ak is
where the parameter q k 0 represents a constant bias in the exponent. The logistic function g(x) = 1/(1 + exp[-x]) is used to restrict the probability to lie between 0 and 1. (This function is shown in Figure 1.10.) In terms of log-odds,
which shows how each parent z j ak independently increases or decreases the log-odds for z k, depending on the sign of q kj.


Page 23
Sometimes, for the sake of notational simplicity, I will assume that the set of parents for each variable z k is specified by some parameter constraints. Assume without loss of generality that for a given network, the random variables z have an ancestral ordering z1, z2, . . . , z N . I take P(z k|ak,q k) = [image], where in the second expression the parameters are constrained so that the function does not depend on nonparents. Also, in order to succinctly account for the bias, I will usually assume that there is a dummy variable z0 that is set to z0 = 1. (Thus the notation q k 0 for the bias in the summations above.) Using these notational simplifications and using g(·) for the logistic function, the sigmoidal model described above can be written
where q k j is set to 0 for each nonparent z j .
1.2.8 Example 2: The bars problem
Bayesian networks provide a useful framework for specifying generative models. A generative model can be used to generate data vectors that exhibit interesting structure. The generative models discussed in this book can also be used for pattern classification and data compression, in the fashion described in Sections. 1.1.1 and 1.1.3. If the Bayesian network is parameterized, we can estimate the parameters of the network from a training set by making the generative distribution "close" (say, in the Kullback-Leibler pseudo-distance) to the training set distribution. We hope that in this fashion, we can extract the "true" underlying generative process, or at least one that is equally efficient at describing the data.
For example, the 4 ´ 4 binary images shown in Figure 1.11a were generated by first selecting an orientation (horizontal or vertical) with equal probability, and then randomly instantiating each of the four possible bars with that orientation with probability 0.5. (The all-on images were removed from the training set, since the orientation of the bars in an all-on image is ambiguous.) Using the Helmholtz machine and the wake-sleep algorithm (described in Sections. 2.4 and 3.4.4), I fit the parameterized network shown in Figure 1.11b to a large training set of 2 ´ 106 images produced in this way. The network has three layers of binary variables: 1 in the top layer, 8 in the middle layer, and 16 in the visible layer (the image). The variables in adjacent layers are fully-connected, and the conditional distributions are modelled using logistic regression, as described in the previous section. After parameter estimation


Page 24
Figure 1.11
(a) Examples of training images typical of the "bars problem". (b) The
graph for a parameterized Bayesian network that was estimated from a
large training set using the wake-sleep algorithm. Edges that terminate on
the box are connected to all vertices within the box. (c) The parameters
of the Bayesian network clearly show that the network has learned the
notion of horizontal and vertical bars (see the main text for a more
complete description).
(see Hinton, et al. [1995] for details), ancestral simulation of the network produces output images that are indistinguishable from the training images.
After learning, the bias for the top-layer variable is nearly zero, so that under the joint distribution it has the value 1 as often as it has the value 0. The values of the other parameters are depicted in Figure 1.11c. The eight large blocks on the left show the parameters associated with the connections that feed into and out of the middle-layer variables. The bias for a variable is shown by the small black or white square on the top right of the block for each middle-layer variable. Positive parameters are white, negative parameters are black, and the area of the square is proportional to the magnitude of the parameter. (The largest parameter shown in the figure is 14.1.) The parameter associated with the connection from the top-layer variable to a middle-layer variable is shown by the small square on the top left of the block for each middle-layer variable. Finally, the parameters associated with the connections from a middle-layer variable to the visible variables are shown by the 4 ´ 4 grid of squares in the


Page 25
block for each middle-layer variable. The biases for the 16 visible variables are shown by the 4 ´ 4 grid of squares on the far right of Figure 1.11c.
It is clear from these parameter values that each middle-layer variable represents the presence (value of 1) or absence (value of 0) of a particular horizontal or vertical bar. If the top-layer variable is 1, the probability that a horizontal bar is present is nearly 0, since the biases for these variables are nearly 0 and the parameters that connect these variables to the top-layer variable are large and negative. On the other hand, if the top-layer variable is 0, the probability that a horizontal bar is present is 0.5. In this way, the network captures the true generative model that produced the training data.
In Chapters 3 and 5, I show how Monte Carlo inference, variational inference, and Helmholtz machines can be used to fit Bayesian networks to training data for the purposes of pattern classification, unsupervised learning, and data compression.
1.3 Organization of this book
In the remainder of this book, I use graphical models (mostly Bayesian networks) as a platform to develop algorithms for pattern classification, unsupervised learning, data compression, and channel coding. The last of these problems is quite different from the former ones, since we will usually design an error-correcting code using a graphical model and then use probabilistic inference to perform decoding. On the other hand, for pattern classification, unsupervised learning, and data compression, we will usually estimate a parameterized Bayesian network from some training data and then use probabilistic inference to classify a new pattern, perform perceptual inference, or produce a source codeword for a new pattern.
In Chapter 2, I discuss different ways to perform probabilistic inference, including probability propagation, Markov chain Monte Carlo, variational optimization, and the Helmholtz machine.
Several types of Bayesian networks that are suitable for pattern classification are presented in Chapter 3. I show how Markov chain Monte Carlo, variational optimization, and the Helmholtz machine wake-sleep algorithm can be used for probabilistic inference and parameter estimation in these networks. Based on a digit classification problem, I compare the performances of these systems with several standard algorithms, including the k-nearest neighbor method and classification and regression trees (CART).
Learning to extract structure from data without using a supervised signal such as class identity is another interesting parameter estimation problem. In Chapter 4, I examine unsupervised learning in Bayesian networks that have


Page 26
binary-valued and real-valued variables. It turns out that simple algorithms in neural-network type Bayesian networks can be used to learn hierarchical structure from data, such as noisy images of horizontal and vertical bars.
In Chapter 5, I consider the problem of how to efficiently compress data using Bayesian networks with hidden variables. When there are hidden variables, a Bayesian network may assign many source codewords of similar length to a particular input pattern. I present the ''bits-back'' coding algorithm that can be used to efficiently communicate patterns, despite this redundancy in the source code.
In Chapter 6, I present several recently invented error-correcting codes in terms of Bayesian networks and show that their corresponding iterative decoding algorithms can be derived as special cases of probability propagation. In particular, the recently proposed turbodecoding algorithm, which brought researchers a leap closer to (and almost up against) the Shannon limit, is an instance of probability propagation. Motivated by these results and the breadth in perspective offered by graphical models, I present a new class of "trellis-constrained codes", which when iteratively decoded are competitive with turbocodes. I also present two approaches for speeding up a popular class of computationally burdensome iterative decoding algorithms.
Finally, in Chapter 7, I discuss some directions which look promising for further research on graphical models in machine learning and digital communication.


Page 27
2
Probabilistic Inference in Graphical Models
In this chapter, I discuss methods for probabilistic inference that make use of the graphical model description of the joint distribution. Many readers may be aware of how probabilistic inference in a Markov chain is simplified by a chain-type graphical structure that leads to the forward-backward algorithm [Baum and Petrie 1966]. A generalized form of this simplification holds for those graphical models that have only a single path (when edge directions are ignored) between any two vertices. In Section 2. 1, I review an algorithm for "probability propagation", which can be used to infer the exact distributions over individual variables or small groups of variables in such models. It turns out that the generalized Viterbi algorithm is just probability propagation with summations replaced by minimizations (a + b ® min (a, b)) and products replaced by summations (a · b ® a + b).
For networks that have multiple paths between one or more pairs of vertices, probability propagation is not exact. Although there are procedures for attempting to convert an original network to one that is appropriate for probability propagation [Spiegelhalter 1986; Lauritzen and Spiegelhalter 1988], these procedures are not practically fruitful when the number of multiple paths is large. In these cases, approximate inference methods must be used. In Section 2.2, I discuss a Monte Carlo approach to inference, where we attempt to produce a sample from the desired distribution. Histograms based on the sample can then be used to approximate the true marginal distributions of interest. In Section 2.3, I present a variational method for approximate inference. Here, we construct a parameterized approximation to the true distribution and then attempt to optimize the parameters of this variational approximation in order to make it as close as possible to the true distribution. This technique requires that the distribution specified by the Bayesian network can be expressed in a form suitable for mathematical analysis. Finally, in Section 2.4, 1 present the Helmholtz machine. This method can be very efficient, and is tailored to inference in Bayesian networks whose parameters are estimated from data.
2.1 Exact inference using probability propagation (the sum-product algorithm)
In the late 1980's, Pearl [1986; 1988] and Lauritzen and Spiegelhalter [Spiegelhalter 1986; Lauritzen and Spiegelhalter 1988; Lauritzen 1996] independently


Page 28
Figure 2.1
A simple factor graph that describes the factorization
P(s, u, v, w, x, y, z)=¦A (s,u,v)¦B (v,w)¦C (u,x)¦D (u,y)¦E (y,z).
published an exact probability propagation algorithm (a.k.a. the "sum-product algorithm" [Wiberg 1996; Frey et al. 1998]) for inferring the distributions over individual variables in singly-connected graphical models. A singly-connected network has only a single path (ignoring edge directions) connecting any two vertices.
It is easy to see how probability propagation works by studying a simple example. Figure 2.1 shows a singly-connected factor graph that describes the following factorization of the global distribution:
Suppose we wish to compute P(s|z = z') for the observed value z' of z. First, notice that P(s|z = z') = P(s, z = z')/Ss P(s, z = z'), so that we can instead compute P(s, z = z') and then normalize to get P(s|z = z'). If we directly compute Su ,v ,w, x, y P(s, u, v, w, x, y, z = z'), we will have to sum over n u ´ n v ´ n w ´ n x ´ n y terms, where n u is the number of values u can take on, and so on.
In fact, the individual summations over u, v, w, x and y distribute across the products in (2.1); e.g., the summation over y can take place over just the last two factors. One distribution of sums gives


Page 29
Each summation can be viewed as a processing stage that removes a variable and produces an output function. The underbraces indicate what the output of each summation depends on. The principle of probability propagation is to decompose the global sum into products of local sums, so that the computation is tractable. Using this procedure, the total number of additions in the example is n u ´ n v + n w + n x + n y. Whereas the number of additions in the brute force marginalization scales exponentially with the total number of variables, the number of additions in probability propagation scales linearly with the exponential of the number of variables in the largest local function. For example, in (2.2), the summation over u and v dominates when the number of values the variables can take on is large (say, n u = n v = n w = n x = n y = 100).
There are essentially two types of computation performed in probability propagation: multiplication of local marginals; and summation over local joint distributions. Within the second set of curly braces in (2.2), the product of the terms in square braces is an example of multiplication of local marginals. The summation over the product of fA(s, u, v) and the two terms in curly braces is an example of a summation over a local joint distribution. Since these operations take place locally in the graph, we can view the edges of the graphical model as message passing channels for inference computations. By passing short real-valued vectors between neighboring vertices in the singly-connected graphical model, the probability propagation algorithm can perform exact probabilistic inference efficiently. One flavor of probability propagation is the generalized forward-backward algorithm, in which messages are passed in a highly regular way.
2.1.1 The generalized forward-backward algorithm
It is easiest to understand the operation of the generalized forward-backward algorithm in factor graphs, although the algorithm can be applied to MRFs and Bayesian networks. In fact, MRFs and Bayesian networks are special cases of factor graphs in this context. So, we will begin with a Bayesian network, convert it to a factor graph, and then see how the generalized forward-backward algorithm works.
To begin with, the singly-connected Bayesian network shown in Figure 2.2a is converted to the factor graph shown in Figure 2.2b. The marginal probabilities for the root variables are simply included in neighboring local functions. Note that the conversion from a Bayesian network to a factor graph can remove local cycles. For example, if there were an edge from z1 to z2 in the original Bayesian network, the singly-connected factor graph shown in Figure 2.2b is still obtained. At this point, we can ignore the edge directions, and simply interpret the parent-child conditional probabilities as local functions in the factor graph.


Page 30
Figure 2.2
(a) shows a singly-connected Bayesian network which is easily
converted to a factor graph, as shown in (b). At this point the
conditional distributions can be interpreted simply as local
functions. A tree is obtained by choosing z 9 as a root vertex
as shown in (c), where the edge directions have been left off.
(d) shows a network fragment with a variable-to-function
message m x®A and a function-to-variable message m A ® x
are used in the generalized forward-backward algorithm.
Next, the factor graph is arranged as a horizontal tree with an arbitrarily chosen "root" vertex on the far right. For example, if the circled vertex z 9 in Figure 2.2b is chosen as the root, we obtain the tree shown in Figure 2.2c. (Imagine the network sits in a viscous fluid and we grasp the root vertex and pull it down and then to the right.) Beginning at the left-most level, messages are passed level by level forward to the root, while being multiplied and marginalized. Each vertex stores its messages for later use. Then, messages are passed level by level backward from the root to the leaves. During these two phases, the observed variables v are held constant and are not marginalized out. The total number of messages passed in this fashion is just twice the number of edges in the graph, since each edge passes a message in both directions. Once both passes are complete, each vertex zi combines its stored incoming messages to obtain P(z i |v).
During each of the forward and backward passes, two types of vector messages are passed: variable-to-function messages consisting of products of local marginals; and function-to-variable messages consisting of sums of distributions over nearby variables. To make this more clear, consider the network fragment shown in Figure 2.2d, with variables x, y and z, and local functions


Page 31
f A(x, y, z), f B(x, . . .) and f C(x, . . . ) (the other variables to which f B (·) and f C (·) are connected are not shown). The variable-to-function message m x®A(x) is computed from
unless x is observed, in which case we set
where x' is the observed value and d (x, x') = 1 if x = x', d (x, x') = 0 if x ¹ x'. The function-to-variable message m x®A(x) is computed from
After the generalized forward and backward phases are complete, the distribution over x conditioned on the observations v is computed from
where b is computed to ensure normalization, i.e., that åx P(x½v) = 1.
2.1.2 The burglar alarm problem
In order to illustrate how the generalized forward-backward algorithm works, consider a variant of the simple "burglar alarm" network described by Pearl [ 1988]. The network describes a shoddy burglar alarm that is sensitive not only to burglars, but also to earthquakes. (Judea Pearl lives in California.) The three binary random variables in the network are b for "burglary", e for "earthquake" and a for ''alarm''. A value of 0 for one of these variables indicates that the corresponding event has not occurred, whereas a value of 1 indicates that the corresponding event has occurred. Figure 2.3a shows the network, which has the following conditional probability relationships:
Suppose that while you are away at a conference, the burglar alarm contacts you by cell phone and informs you that the alarm is ringing (a = 1). We would like to infer the distribution over the two causes to make a well informed decision about whether or not you should be concerned about a bur-


Page 32
Figure 2.3
(a) The Bayesian network for the burglar alarm problem, with variables b
("burglary"), e ("earthquake") and a ("alarm"). (b) a factor graph for the
network in (a). (c) - (e) and (f) - (h) show the messages passed during the
generalized forward and backward passes.
glary. Since this network is quite simple, we can apply Bayes rule P(b, e|a) = [image] to obtain the exact solution,
The most likely explanation for the ringing alarm is that a burglary took place. However, although an earthquake is also a likely explanation, it is relatively unlikely that both a burglar and an earthquake simultaneously caused the alarm.
Now, consider using the generalized forward-backward algorithm for probabilistic inference in this network. We first rewrite the Bayesian network as a factor graph, as shown in Figure 2.3b, where
After we arbitrarily select e as the root, the generalized forward-backward algorithm sends messages in the forward direction as shown in Figures 2.3c-e, and then in the backward direction as shown in Figures 2.3f-h as follows:


Page 33
(c) Since B is connected to b by one edge, the function-to-variable message m B®b is equal to the local function P(b) (see (2.50)): m B®b = (0.9, 0.1).
(d) Since b is connected to B by one edge, m b®A is equal to m B®b (see (2.3)): m b®A = m B®b = (0.9, 0.1). a is observed to be 1, so we set its message to m a®A = (0, 1).
(e) This time, the function-to-variable message m A®e will require a marginalization as indicated in (2.5):
The message m E®e is set equal to the local function P(e): m E®e = (0.9, 0.1).
(f) The message m e®E does not influence the distributions over e and b, so we can neglect computing it. The message m e®A is set equal to m E®e (see (2.5)): m e®A = m E®e = (0.9, 0.1).
(g) Again, the message m A®a is not needed, so we don't compute it. The message m A®b requires a marginalization, as indicated in (2.5):
(h) The message m B®b is not needed so we don't compute it.
Finally, we compute P(b½a = 1) and P(e½a = 1) using (2.6):
These distributions are exactly equal to the marginal posterior distributions computed from (2.8).


Page 34
2.1.3 Probability propagation (the sum-product algorithm)
The highly regular way that messages are passed in the generalized forward-backward algorithm can be relaxed to obtain a more general probability propagation algorithm. It turns out that as long as a few simple rules are followed, messages may be passed in any order (even in parallel) to obtain conditional probabilities. These rules prescribe how the network is to be initialized for propagation, and how messages are created, propagated, absorbed and buffered. Aside from these rules, the formulas for propagating messages are identical to those in (2.3) to (2.6).
Before propagation begins, the network must be initialized. This procedure computes the a priori incoming messages for each vertex, and corresponds to a generalized forward-backward pass without any observations. After initialization, each vertex z i has available its a priori probability P(z i). In some networks (such as those used for channel coding) these probabilities are uniform and so the initialization procedure can be skipped.
Messages are now created in response to observations. If variable y is observed to have the value y', then a message must be sent out on each of the edges connected to y, using (2.4).
Messages are propagated in response to other messages. If variable y receives a message on an edge, y must send out messages on all other edges, and similarly for the function vertices.
Messages are absorbed by vertices that are connected by only a single edge. This rule follows naturally from the propagation rule, since if such a vertex receives a message on its only edge, the vertex is not required to propagate it back.
It is not necessary that messages be propagated without delay. In fact, a vertex may buffer one or more outgoing messages and pass them at any time. (It is usually most convenient to compute them at a later time, too.) For example, if a vertex has just received a message and is about to receive another one, computations can often be saved by waiting for the second message before computing and sending out a set of messages.
At any time during propagation, vertex y can compute a current estimate [image] using (2.6). If the above rules are followed and propagation continues until there are no buffered messages remaining in the network, then the estimates will equal the exact probabilities: .
Instead of a complete initialization, it is possible to simply buffer the initial messages. Since these messages will be propagated eventually, this has the same final result as the initialization procedure described above, although the intermediate probabilities may differ.


Page 35
The generalized forward-backward algorithm described in the previous section can be viewed as a special case of probability propagation. First, the network is arranged as a tree. Then, the initialization messages are buffered, as are the messages that would be created in response to any observations. (At this point, no computations have been performed.) During the forward pass, each right-going message induces a set of buffered left-going messages and a single right-going message. The latter right-going message is passed to the next level, where it too induces a set of buffered left-going messages and a single right-going message. So, once the forward pass is complete, there are no more buffered right-going messages in the network. During the backward pass, each vertex receives a left-going message from its only right-hand edge. Since an incoming message to a vertex never induces an outgoing message on the same edge, the left-going message will induce only a set of buffered left-going messages. So, there will be no buffered messages remaining in the network once the backward pass is complete. Finally, each vertex z i can compute the exact value for P(z i|v).
2.1.4 Grouping and duplicating variables in Bayesian networks
Often, it is possible to convert a multiply-connected Bayesian network to a singly-connected Bayesian network, so that probability propagation can then be applied in a practical manner. To do this, we group variables, until there are no more multiple paths in the network. Graphically, two variables z j and z k are grouped by removing from the graph z j and z k as well as the edges to which they are connected, and then introducing a new vector variable{z j, z k}. The set of parents of the new vector variable is the union of the sets of parents of the two old variables. The set of parents of each child of z j and z k is extended to include both z j and z k . New edges are introduced to reflect these relationships. This grouping operation will produce a valid Bayesian network as long as z k is not an indirect descendent of zj and vice versa. Otherwise, a directed cycle will result from the grouping, violating the requirement that a Bayesian network have no directed cycles. Note that if z k is a child of z j, and at the same time not an indirect descendent, the grouping is still valid, since no directed cycles are produced.
It turns out that this grouping operation also preserves the representational capacity of the network. Any distribution represented by the old network can be represented by the new one. In fact, all of the conditional probabilities P(z k|ak) in the new network will be the same as in the old network, except the ones that involve either of the grouped variables. The latter conditional probabilities can quite easily be derived from the old ones.
Although grouping variables may help to produce a singly-connected network to which probability propagation can be applied, the grouping operation


Page 36
Figure 2.4
Transforming the multiply-connected Bayesian network for a recursive
convolutional code (a) into a singly-connected network (c) involves
variable duplication (a) - (b) and grouping (b) - (c).
also hides the structure that makes probability propagation an attractive inference method in the first place. So, it is important to produce minimal groupings and retain as much of the structure as possible. Not surprisingly, by grouping variables, any network can be made singly-connected—simply group all variables together into a single vertex. However, extreme groupings of this sort usually eliminate too much structure. "Probability propagation" for the single vertex is equivalent to manipulating the full joint distribution, which in most practical cases is unwieldy.
Another useful operation is duplicating variables. A variable zj can be duplicated by adding an extra variable z N +1 to the network, and creating the following new parent-child relationships: . This procedure is especially useful in combination with grouping, since although we may wish to group z j and z k in order to make the network singly-connected, we may also wish to graphically distinguish z j from the vector variable{z j, z k}.
For example, the recursive convolutional code network shown in Figure 1.9a can be derived from the more natural recursive convolutional code network shown in Figure 2.4a. The latter network explicitly shows the dependence of the encoder state variable s k on the previous information symbol u k - 1 and the previous state s k - 1, as well as the dependence of the encoder output x k on u k and s k. This network is multiply-connected, so probability propagation cannot be used to compute P(u k|y) for maximum a posteriori information symbol decoding. To convert the network to a singly-connected one, we first duplicate the information symbols (so that they are graphically distinguished in the final network) as shown in Figure 2.4b. Then, we group pairs of information symbols and state variables as shown by the dashed loop, producing the singly-connected network shown in Figure 2.4c. Note that by grouping variables in


Page 37
this way, the number of values that each new state s'k can take on is increased by a factor of two. Although in many cases grouping can be used to produce a tractable network, there are cases where it is impossible to find an appropriate grouping. In fact, it turns out that inference in multiply-connected networks is a very difficult problem, in general.
2.1.5 Exact inference in multiply-connected networks is NP-hard
Probability propagation is an exact method of inference for singly-connected graphical models. Cooper [1990] has shown that probabilistic inference in Bayesian networks is in general NP-hard. The same is true for MRFs and factor graphs. Summations relevant to inference, such as the ones in (1.3), (1.7), (1.11) and (1.13), contain an exponential number of terms and it appears that in general these summations cannot be simplified. Researchers have thus focused on developing exact inference algorithms for restricted classes of graphical models (e.g., probability propagation for singly-connected networks), and on developing approximate inference algorithms for networks that are intractable (assuming P ¹ NP). In fact, Dagum [1993] (see also Dagum and Chavez [1993]) has shown that for general Bayesian networks, approximate inference to a desired number of digits of precision is NP-hard. (I.e., the time needed to obtain an approximate inference that is accurate to n digits is believed to be exponential in n.)
One obvious approach to approximate inference in a multiply-connected graphical model is to use the probability propagation algorithm while ignoring the fact that the network is multiply-connected. Each vertex propagates messages as if the network were singly-connected. In this case, the propagation procedure will never terminate, because there will be loops in which messages will endlessly circulate. Although this method has provided excellent results in the area of channel coding, it is frowned upon in other areas (such as medical diagnosis) because there is little theoretical understanding of the behavior of this iterative procedure.
Another disadvantage of the probability propagation algorithm is that it is cumbersome for inferring the joint distribution over several variables (e.g., u1 and u4 in Figure 2.4c). This inference is accomplished by first computing the distribution over u1 given the observations, using one forward-backward sweep. Then, the distribution over u4 given the observations and each of the possible values for u1 is computed using one forward-backward sweep for each possible value for u1. (Notice that these sweeps may be partial, since they need only take into account the effects of clamping u1 to different values.) If the variables of interest have n possible configurations, roughly n (possibly


Page 38
partial) forward-backward sweeps are needed. If we cannot afford the time to perform these sweeps, a faster approximate algorithm may be desirable.
In the following sections, I describe several other approaches to approximate probabilistic inference, including Monte Carlo, variational inference and Helmholtz machines. In these sections, I focus on inference in Bayesian networks.
2.2 Monte Carlo inference: Gibbs sampling and slice sampling
The Monte Carlo method [Hammersley and Handscomb 1964; Kalos and Whitlock 1986; Ripley 1987] makes use of pseudo-random numbers in order to perform computations. Monte Carlo inference uses random numbers in order to perform inference in a graphical model that describes a joint distribution P(z). If we can somehow obtain a reasonably large sample from the distribution P(h|v) over some unobserved hidden variables h Í z in a graphical model, given some observed visible variables v Í z, then relative frequencies can be used for approximate inference.
2.2.1 Inference by ancestral simulation in Bayesian networks
One brute force Monte Carlo approach is to simply simulate the network using the ancestral ordering. Then, we extract from the sample all those vectors that have the desired value for the component v. Next, we compile a frequency histogram for the different values that h can take on. Although this approach is sometimes useful (notably, when using a small network to verify that a more sophisticated inference method works), in general it is not computationally efficient. The problem is that the value of v that we wish to condition on may occur extremely rarely, so that an inordinate sample size must be used in order to obtain useful results.
If we are interested in a subset hI Í h of the hidden variables, it so happens that in some cases ancestral simulation can be used to obtain a sample from P(hI|v) in an efficient manner. In general, if the parents of the visible variables are dependency-separated (see Section 1.2.5) from the hidden variables of interest by the visible variables, then ancestral simulation can be used to obtain a sample from P(hI|v). If the visible variables have no parents, then it follows trivially that the variables in this null set are dependency-separated from the hidden variables. Using the ancestral ordering, a value is drawn for each hidden variable given its parents. After one complete sweep, the value for [image] will be an unbiased draw from P(hI|v).
Ancestral simulation can be used for inference in this case because, if the condition holds, then every path connecting each variable in hI to the parents


Page 39
Figure 2.5
An example of a Bayesian network.
of each visible variable is blocked. This means that disconnecting each visible variable from its parents will not change the distribution P(hI v). Since each visible variable will then have no parents, its value can be included as fixed constant in the conditional probability functions for its children. We are then left with a new Bayesian network that describes a distribution P'(h) over the variables h that were not observed in the original network. Although in general P'(h) ¹ P(h|v), as shown above we have P'(hI) = P(hI|v). So, we may simply use ancestral simulation in the new network to obtain samples from P(hI v). Notice that ancestral simulation in the new network is equivalent to ancestral simulation for the unobserved variables in the original network.
For example, suppose that in the multiply-connected network with ancestral ordering z1, z2, z3, z4, z5, z6, z7 shown in Figure 2.5, the set of visible variables is {z1, z4} and we would like to infer the distribution over the subset of hidden variables {z6, z7}. Since {z6, z7} is dependency-separated by {z1, z4} from the parents {z2} of {z1, z4}, we can estimate P(z6, z7 | z1, z4) by ancestral simulation. We draw a value for z2, then for z3 given z1, then for z5, then for z6 given z3 and z4, then for z7 given z4 and z5. We can estimate P(z6, z7 | z1, z4) by repeating this procedure over and over while building up a histogram for {z6, z7}.
Notice that in the above example, it was not really necessary to draw a value for z2, since for the ancestral simulation method to work, it was required that z2 be dependency-separated from {z6, z7} by {z1, z4}; i.e., given {z1, z4}, z2 does not influence {z6, z7}. In general, values need only be drawn for those variables in the ancestral ordering that are not dependency-separated from hI by v. If we can easily identify these, we can save simulation computations.
2.2.2 Gibbs sampling
When inference by ancestral simulation is not possible, Markov chain Monte Carlo is often used (see an excellent review of these methods by Neal [1993]). Given v, a temporal sequence h(1), h(2), . . . of the hidden variable values is


Page 40
produced by simulating a Markov chain whose stationary distribution is carefully constructed (e.g., as described below) to be equal to P(h½v). By collecting these values over time, an approximate sample is obtained. Ideally, the chain is run long enough so that equilibrium is reached. In practice, the Markov chain may be terminated before equilibrium is reached, so that the simulation time can be kept within a reasonable limit. Once collected, the sample can be used to produce a frequency histogram of the variables of interest in h.
The Gibbs sampling algorithm is the simplest of the Markov chain Monte Carlo methods, and has been successfully applied to Bayesian networks [Pearl 1987; Pearl 1988; Neal 1992] as well as other graphical models [Geman and Geman 1984; Hinton and Sejnowski 1986]. In this algorithm, each successive state h(r) is chosen by modifying only a single variable in the previous state [image] . The variables are usually modified in sequence. If at time [image], we have decided to modify z k Î h, then we draw a value [image] from
Usually, we cannot obtain this distribution directly, but instead must first compute the joint probability which is proportional to the conditional probability:
where the constant of proportionality does not depend on [image] . The joint probability can usually be computed easily from (1.19). If [image] is discrete, we compute the joint probability for each value that it can take on, normalize these values, and then randomly draw a value [image] from this normalized distribution. When [image] is a continuous random variable, it can be quite difficult to draw a value from its distribution. Efficient sampling methods for several special types of continuous parametric distribution are given in Devroye [1986] and Ripley [1987]. In order to draw values from other types of distribution, more sophisticated techniques such as adaptive rejection sampling [Gilks and Wild 1992] must be used.
2.2.3 Gibbs sampling for the burglar alarm problem
In order to illustrate how Gibbs sampling works, I use the simple burglar alarm problem presented in Section 2.1.2, whose Bayesian network is shown in Figure 2.6a.
In order to perform Gibbs sampling, we need the probabilities for each of the hidden variables conditioned on all the other variables. Since these conditional probabilities are proportional to the joint probabilities, we can compute them in the following way, using P(b=1½e=0, a = 1) as an example:


Page 41
Figure 2.6
(a) The Bayesian network for the burglar alarm problem, with variables
b (''burglary''), e ("earthquake") and a ("alarm"). (b) 400 steps of Gibbs
sampling for the variables b and e when the alarm is observed to be
ringing (a = 1).
where P(b,e,a) = P(a|b,e)P(b)P(e) is the joint distribution determined from the network specification in (2.7). Similarly,
Gibbs sampling proceeds by alternately visiting b and e, while sampling from P(b|e, a = 1) and P(e|b, a = 1) using the above formulas. Figure 2.6b shows the values of b and e for 400 steps of Gibbs sampling, starting from an initial configuration (b = 0, e = 0). (In each step, one variable is updated.) The Markov chain shows that the configurations (b= 0, e = 0) and (b= 1, e = 1) are unlikely compared to (b= 1, e = 0) and (b= 0, e = 1). The correct probabilities


Page 42
in (2.8) can be approximated using the relative frequencies computed from this chain:
These are quite close to the correct values given in (2.8). Usually, an initial segment of the Markov chain is discarded when computing these statistics. The motivation for this procedure is that we would like to have samples that are typical of the equilibrium distribution, not the initial configuration.
2.2.4 Slice sampling for continuous variables
In networks with continuous variables, it is often not an easy task to sample from the conditional distribution of each hidden variable, as Gibbs sampling requires. Unlike the case for discrete variables, it is usually not possible to compute the joint distribution for every configuration of a hidden variable. There are infinitely many configurations, and it is often practically impossible to determine an effective discretization. Methods for sampling from continuous distributions include the Metropolis algorithm [Metropolis et al. 1953; Neal 1993] and hybrid methods that use "momentum" in order to help search the configuration space [Duane et al. 1987; Neal 1993]. In this section, I review a technique called slice sampling [Neal 1997; Frey 1997a], that can be used for drawing a value z from a univariate probability density p(z)—in the context of inference, p(z) is the conditional distribution [image] Slice sampling does not directly produce values distributed according to p(z), but instead produces a Markov chain that is guaranteed to converge to p(z). At each step in the sequence, the old value zold is used as a guide for where to pick the new value znew. When used in a system with many variables, these updates may be interleaved for greatly improved efficiency.
To perform slice sampling, all that is needed is an efficient way to evaluate a function f(z) that is proportional to p(z)—in this application, the easily computed joint probability [image] is appropriate. Figure 2.7a shows an example of a univariate distribution, p(z). The version of slice sampling discussed here requires that all of the probability mass lies within a bounded interval as shown. To obtain znew from znew, ¦(znew) is first computed and then a uniform random value is drawn from [0, ¦(znew)]. The distribution is then horizontally "sliced" at this value, as shown in Figure 2.7a.


Page 43
Figure 2.7
Slice sampling. After obtaining a random slice from the density (a), random
values are drawn until one is accepted. (b) and (c) show two such sequences.
Any z for which f(z) is greater than this value is considered to be part of the slice, as indicated by the bold line segments in the picture shown at the top of Figure 2.7b. Ideally, znew would now be drawn uniformly from the slice. However, determining the line segments that comprise the slice is not easy, for although it is easy to determine whether a particular z is in the slice, it is much more difficult to determine the line segment boundaries, especially if the distribution is multimodal. Instead, a uniform value is drawn from the original interval as shown in the second picture of Figure 2.7b. If this value is in the slice it is accepted as znew (note that this decision requires an evaluation of f (z)). Otherwise either the left or the right interval boundary is moved to this new value, while keeping znew in the interval. This procedure is repeated until a value is accepted. For the sequence in Figure 2.7b, the new value is in the same mode as the old one, whereas for the sequence in Figure 2.7c, the new value is in the other mode. Once znew is obtained, it is used as znew for the next step. As shown in Neal [1997], this procedure satisfies detailed balance and therefore gives the desired stationary distribution p(z).
2.3 Variational inference
In contrast to both the rather unprincipled approach of applying probability propagation to multiply-connected networks, and the computationally intensive stochastic approach of Monte Carlo, variational inference is a nonstochastic technique that directly addresses the quality of inference. In the Bayesian network literature, variational inference methods [Saul, Jaakkola and Jordan 1996; Ghahramani and Jordan 1998; Jaakkola, Saul and Jordan 1996] were introduced as an alternative variation on the central theme of Helmholtz machines [Hinton et al. 1995; Dayan et al. 1995], which are described in Sec-


Page 44
tion 2.4. However, I will present variational inference first, because it is simpler to understand.
Suppose we are given a set of visible variables v Í z. (This set may includes different variables on different occasions.) In order to solve the inference problem of estimating P(h½v), we introduce a parameterized variational distribution Q(h½x) that is meant to approximate P(h½v). The most appropriate form of this distribution will depend on many factors, including the network specification and the quality of inference desired. Next, the distance between P(h½v) and Q(h½x) (e.g., Euclidean, relative entropy) is minimized with respect to x, either directly or by using an optimization technique such as a Newton-like method or a conjugate gradient method [Fletcher 1987]. Once optimized, the distribution Q(h½x) is used as an approximation to P(h½v).
A main advantage of variational inference over probability propagation in multiply-connected networks is the explicit choice of a distance measure that is minimized. Although probability propagation is optimal for singly-connected networks, there is very little known theoretically about the quality of inference that results when the network is multiply-connected. On the other hand, there is no general guarantee that in multiply-connected networks, variational methods will perform better than probability propagation. An example where probability propagation in multiply-connected networks works very well for practical purposes is the celebrated turbodecoding algorithm for error-correcting coding [Berrou, Glavieux and Thitimajshima 1993; Frey and Kschischang 1997].
Compared to Monte Carlo, variational inference may provide the designer with a more structured approach to choosing a computationally tolerable approximation to P(h½v). However, variational methods do not usually provide a means to obtain exact inference. Also, variational inference can only be applied when the network is well-tailored to a sensible distance measure along with a fruitful form of variational distribution. (For example, the majority of work on variational methods for Bayesian networks to date has focused on networks that are parameterized.) In contrast, Monte Carlo methods can be applied to any Bayesian network, and can be designed so that they are guaranteed to converge to the correct solution.
2.3.1 Choosing the distance measure
Depending on the particular problem, different measures of distance may be appropriate. For example, in the case of hard-decision classification and hard-decision channel coding, a binary distance is ideal. Under this distance, the distributions are identical if they lead to the same decisions. Otherwise, the distance is incremented for each incorrect decision. In practice, this distance must be softened in order to use continuous optimization methods.


Page 45
Figure 2.8
The effect of using (a) Dq||p versus (b) D p||q when fitting a variational
distribution q(z|x) that is unimodel to a bimodal distribution p(z).
As another example, we will see in Chapters 3 and 5 that for pattern classification and data compression, the appropriate "distance" is the Kullback Leibler divergence, or relative entropy, between Q(h|x) and P(h|v)
Notice that this is not a true distance since it is not symmetric: DQ||P ¹DP||Q,where
(For density functions, the summations are replaced by integrals.)
The choice of whether to use D Q||P or D P||Q depends on our objective. The former places emphasis on not inferring unlikely values of h at the cost of not inferring some of the likely values, whereas the latter places emphasis on inferring all likely values of h at the cost of inferring some of the unlikely values. For example, consider a real-valued univariate probability density p(z) over z that has two modes, as shown in Figure 2.8. Suppose the variational distribution q(z|x) is a Gaussian with x consisting of a mean and a variance. Figure 2.8a shows the optimum variational distribution that is obtained by minimizing D q||p, whereas Figure 2.8b shows the optimum variational distribution that is obtained by minimizing D p||q.
Notice that in order to compute DQ||P in (2.18), we need P(h|v), which is what we were after in the first place. So, in practice, we usually minimize the following free energy function:


Page 46
Notice that minimizing F Q ½½P with respect to x gives the same set of parameters as minimizing D Q ½½P, since log P(v) does not depend on x. In order to compute F Q ½½P, we only need P(h,v), which is readily available in Bayesian networks. (P(h, v) is not easy to compute in other types of graphical models, such as Markov random fields.)
2.3.2 Choosing the form of the variational distribution
The form of Q(h½x) will strongly influence the quality of the variational inference as well as the tractability of computing the distance and its derivatives (which may be needed for the optimization procedure). Exact inference can be achieved in principle by associating one parameter xh with each state of the hidden variables h, where xh is meant to be an estimate of P(h|v). However, computing the distance will require an explicit summation over all possible states of the hidden variables. The number of terms in this sum equals the number of possible configurations of the hidden variables, so this approach will only be tractable when there are not many configurations of the hidden variables. In fact, in most cases the above procedure will not be any more computationally efficient than directly computing P(h½v) using Bayes' rule.
We would like to choose Q(h½x) so that the effect of the hidden variables h in the distance measure can be integrated out either analytically or using a reasonably small number of computations. In this way, the distance and its gradients can be determined without having to numerically examine each possible state of the hidden variables h.
2.3.3 Variational inference for the burglar alarm problem
In this section, I illustrate variational inference using the burglar alarm network described in Section 2.2.3. One type of variational distribution that is often used is the product form distribution. Under this variational distribution, the hidden variables are independent. For continuous variables, further assumptions may be needed regarding the distributions for each hidden variable (e.g., see Section 4.2). For the binary burglar b and earthquake e variables in the burglar alarm network, we can specify an arbitrary product-form distribution using the parameters x1 and x2 for the probabilities that b = 1 and e = 1 respectively. That is,
Inserting this variational distribution into (2.20), and using the simple binary entropy function H(x1) = x1 log x1 - (1 - x1) log(1 - x1), we get


Page 47
Notice that the product form of Q(b, e|x) was used to simplify the first term of the second equality.
At this point, without any further restrictions, we have not gained any computational advantage by using the variational approach. To compute F Q||P and its derivatives, we must still examine all possible configurations of the hidden variables to compute the expectation of log P(b, e, a = 1). In order to make profitable use of variational inference, log P(b, e, a = 1) must have a form that makes the computation of F Q||P easy. It turns out that the conditional probabilities (2.7) for the burglar alarm network were obtained from
So, the joint distribution P(b, e, a = 1) can be written
Substituting this into (2.22), we get
Notice that the hidden variables b and e do not appear in this final expression. Because of the product form of Q(b, e|x) and the exponential form of P(b, e, a = 1), we were able to integrate them out.


Page 48
Figure 2.9
The contours of FQ||P for a variational
technique applied to the burglar alarm
problem. The global minimum occurs
at x1 = 0.951, x2 = 0.186.
Figure 2.9 shows a plot of the contours of F Q||P as a function of x1 and x2. The global minimum occurs at x1 = 0.951 x2 = 0.186, which means the inference estimates are [image] (b = 1|a = 1) = 0.951 and [image] (e = 1|a = 1) = 0.186. These estimates clearly favor a burglar as the cause of the alarm. Recall that Gibbs sampling allowed us to estimate covariance statistics between the two hidden variables. Variational inference does not readily produce those estimates. However, compared to the marginal probabilities P(b = 1|a = 1) = 0.751 and P(e = 1|a = 1) = 0.349 produced by the probability propagation algorithm, the variational method places more emphasis on the more likely cause b. In this sense, the variational technique produces a product form distribution that reveals covariance better than the marginals, produced, say, by probability propagation. For example, using only the marginal probabilities produced by propagation, we might conclude that the probability that both a burglar and an earthquake occurred is 0.751 ´ 0.349 = 0.262. In contrast, using the probabilities produced by the variational method gives 0.951 ´ 0.186 = 0.177, which is closer to the correct value of 0.116 given in (2.8).
In this case, because the burglar alarm network is so small, the analytic form of F Q||P in (2.25) is not much simpler than the expression that would be obtained if (2.21) were substituted into (2.22) and explicit summation over all values of b and e were performed. However, for larger networks, exponential computational savings may be achieved by using conditional distributions that lead to simple forms of log P(h, v).


Page 49
2.3.4 Bounds and extended representations
In practice, the form of log P(h,v) is often not simple, so that a straightforward variational approach cannot be attempted. In these cases, it may be possible to derive an upper bound on the distance that does not depend on h, and then try to minimize the bound instead of the distance itself [Saul, Jaakkola and Jordan 1996]. Effectively, we approximate log 1/P(h,v) with an upper bound that can be integrated analytically.
Alternatively, we may express each conditional distribution P(zk|ak) in terms of conditional distributions over an extended set of variables [Jaakkola, Saul and Jordan 1996]. For example, P(zk|ak) might be the marginal distribution of P(zk, yk|ak) where yk is part of the extended representation. Let yH be the extension variables associated with the variables in h. It is sometimes possible to introduce a variational distribution Q(h, yH|x) over the extended representation for which h and yH can be integrated out in the distance measure.
2.4 Helmholtz machines
One of the main drawbacks of Markov chain Monte Carlo inference and variational inference is that for complex networks, each time a set of variables is observed, either a computationally taxing Markov chain must be simulated, or a high-dimensional optimization must be performed to find the best variational distribution. The essential problem, of course, is that the optimal distribution over h is different for different values v of the visible variables. A Helmholtz machine [Dayan et al. 1995; Hinton et al. 1995] tackles this problem by coupling the original generative network with a recognition Bayesian network that is meant to be capable of quickly producing an estimate of, or an approximate sample from, P(h½v). This recognition network essentially replaces the variational optimization needed for variational inference. It is called a "recognition" network because it is meant to recognize the hidden variable values, or "causes", that are responsible for the values of the visible variables.
As described above, the job of the recognition network is to quickly produce an approximation to P(h½v). Obviously, the recognition network must be different from the generative network, or the inference could be done directly on the generative network. I will highlight this difference by labeling the recognition distribution with Q. So, the recognition network is used to compute Q(h½v), which is an approximation to P(h½v) as given by the generative network. Various types of recognition network are described below, but they all share a common property. Since the recognition network is a Bayesian network, we cannot expect to be able to quickly compute Q(h½v) for arbitrary


Page 50
Figure 2.10
(a) An example of a factorial recognition network.
(b) A factorial recognition network for the burglar alarm problem.
sets h and v. In fact, I will usually assume that the set of visible variables is the same for each inference case, although, of course, the values for the visible variables may differ from case to case. This restriction is the main disadvantage of the Helmholtz machine compared to Monte Carlo inference and variational inference, which usually place no restrictions on which variables are observed.
2.4.1 Factorial recognition networks
To ensure that the inference process is fast, we ought to design the recognition network so that the computation of Q(h|v) can be carried out efficiently. The simplest recognition network in this sense is one for which each variable in h is dependency-separated from each other variable in h by the visible variables v. In other words, given the visible variables, the hidden variables are independent. I will refer to such a network as a factorial recognition network, since given the visible variables, the distribution over the hidden variables can be factored into a product of probabilities:
A factorial recognition network with h = {z1, z2, z3} is shown in Figure 2.10a. By condition 2 in Section 1.2.5, variables in h are dependency separated by v = {z4, z5, z6, z7}.
In many cases, the product form approximation given in (2.26) is not very close to P(h|v). However, it is the easiest network to design or estimate, and because the hidden variables are independent given the visible variables, it is computationally efficient for inference.


Page 51
Figure 2.11
(a) An example of a nonfactorial recognition network.
(b) A nonfactorial recognition network for the burglar alarm problem.
Figure 2.10b shows a factorial recognition network for the burglar alarm problem. The recognition distribution is given by Q(b, e|a) = Q(b|a)Q(e|a), and so is limited to the sane inference estimates as the variational technique described in Section 2.3.3. Namely, the factorial recognition network cannot capture the covariance between the two causes. Can we design a recognition network that can give better estimates? The answer is ''yes'', by using a nonfactorial recognition network.
2.4.2 Nonfactorial recognition networks
Although it is easy to imagine situations where a factorial recognition network will suffice, for the burglar alarm problem discussed above we saw that a factorial recognition network could not capture the covariance between the two causes of the alarm. In this section, I describe nonfactorial recognition networks that are more powerful than factorial ones.
A nonfactorial recognition network can represent a distribution where at least one variable in h is not dependency-separated from at least one other variable in h by the visible variables v. Of course, there are many ways to make a network nonfactorial. For example, a nonfactorial recognition network is obtained by making some hidden variables depend on other ones in addition to the visible variables. Figure 2.11a shows a fully-connected nonfactorial recognition network, which can be contrasted with the factorial network in Figure 2.10a.
Another way to produce a nonfactorial recognition network is through the use of auxiliary variables or dangling units [Dayan and Hinton 1996]. These variables do not influence the output of the generative model, but help facilitate inference in the recognition network. For example, an auxiliary variable in the recognition network can be used to choose between two or more modes.


Page 52
2.4.3 The stochastic Helmholtz machine
Suppose we are interested in only one of the hidden variables, and we would like to obtain its distribution given the visible variables, after marginalizing out the other hidden variables. For a factorial recognition network, each hidden variable is independent given the visible variables. So, the marginal distribution is obtained simply by ignoring the other hidden variables. In fact, the marginal distribution for z k Î h in this case is Q(z k|v), which is part of the recognition network specification.
Such a simple procedure for marginalization is not in general available for nonfactorial recognition networks. In these networks, the hidden variables are not independent given the visible variables. However, Monte Carlo provides an easy way to estimate marginal statistics. If we can obtain a sufficiently large sample from the recognition network, the distribution for z k can be approximated by constructing a histogram for z k alone. Of course, we could directly apply Monte Carlo methods such as Gibbs sampling (Section 2.2.2) to the generative network. However, the hope is that we can carefully design the nonfactorial recognition network so that it is better suited to Monte Carlo than the generative network. In fact, we can avoid complicated Markov chain Monte Carlo by using a recognition network for which ancestral simulation (see Section 1.2.4) can be used.
In general, recognition networks can be either factorial or nonfactorial and stochastic or nonstochastic. Here, "nonstochastic" refers to the way the recognition network is used, not to what the network represents. All Bayesian networks represent a stochastic phenomena, but not all networks are used with Monte Carlo. A factorial recognition network can easily be operated stochastically, simply by choosing each hidden variable zk from its distribution Q(z k|h). A nonfactorial recognition network can easily be operated nonstochastically, since the joint distribution over the hidden variables factors and the marginal distribution for each hidden variable is readily available. However, a nonfactorial recognition network usually cannot be operated nonstochastically. As described above, the dependencies between the hidden variables makes this difficult. However, there are special cases where nonfactorial recognition networks can be operated nonstochastically. In particular, recognition networks that can be viewed as a mixture of factorial networks can be operated nonstochastically with relative ease.
2.4.4 A recognition network that solves the burglar alarm problem
In many cases, a simple nonfactorial recognition network can be used to represent covariances between hidden variables. A nonfactorial recognition network for the burglar alarm problem is shown in Figure 2.11b. The difference


Page 53
between this network and the factorial one in Figure 2.10b, is that e now depends on b as well as a. The conditional distributions for a recognition network that performs exact inference are
Sampling hidden variables using ancestral simulation in this network is actually more efficient than using Gibbs sampling in the generative network, as described in Section 2.2.3. (The computational savings are quite low in this case, because there are only two hidden variables.)
The joint distribution over the hidden variables given a = 1 can be computed from Q(b,e½a=1) = Q(e½b,a=1)Q(b½a=1):
These probabilities are identical to the probabilities in (2.8) for exact inference.


Page 55
3
Pattern Classification
Automated methods for making decisions based on inputs play a very important role both in engineering applications and in helping us understand how biological systems respond to their environments. As many engineers and cognitive scientists will attest, the terms "input" and "decision" for this pattern classification problem are not clearly defined in theory. In practice, the problem is usually decomposed through design and analysis. The input to the classifier is provided by a preprocessor that transcribes part of the physical state of the world. Different preprocessors are appropriate for different classifiers, and often an iterative process is used to find the optimal preprocessor-classifier pair for a given problem. In general, the preprocessor uses simple statistical and signal processing techniques, whereas the classifier is left with the "hard" problem of coming up with decisions.
A very simple method for making hard decisions is the nearest neighbor classifier. This classifier keeps a database of labeled training patterns. Given a test pattern, the nearest neighbor classifier outputs the class of the pattern in its database that is "closest" to the test pattern. Any distance metric may be used, but typically Euclidean distance or one of its generalizations are used. Figure 3.1 shows a selection of normalized and quantized 8 ´ 8 binary images of handwritten digits made available by the US Postal Service Office of Advanced Technology. A database with a total of 7000 patterns was constructed with 700 patterns from each digit class. Using nearest neighbor classification, a misclassification rate of 6.7% was obtained on a test set of 4000 patterns. Slightly better results can be achieved by using the k-nearest neighbor method. This method picks the most common class in the k training patterns that are closest to the test pattern.
One interesting property of the k-nearest neighbor method is that it is a consistent classifier. That is, as the number of training cases T tends to infinity, the decisions produced by the k-nearest neighbor method (with, e.g., k = [image])become Bayes optimal. However, although k-nearest neighbor classification works quite well when a large training set is available, it performs poorly when training data is limited. Figure 3.2 shows a training set consisting of two classes with 30 2-dimensional real-valued patterns in each class. Suppose we wish to classify the indicated test point. The nearest neighbor method will choose class A. In fact, just as our intuition tells us, the test point was drawn from class B. If a k-nearest neighbor classifier is used, class A will consistently be erroneously chosen for sensible values of k.


Page 56
Figure 3.1
450 examples of 8 ´ 8 binary images of handwritten digits.
The above example illustrates a fundamental flaw with the nearest neighbor approach—namely, that it does not use global structure. Viewing the data from class B with a local (narrow) "window", the test pattern seems very unlikely. However, a more global examination of the data from class B leads us to believe that the data comes from a roughly sinusoidal manifold, and that just by chance there isn't any training data for this class in the central region of the figure. Under this view, the test pattern is much more likely. An even more global examination indicates that the two classes of data are probably similar, except for the fact that they lie on manifolds that are relatively inverted. As a result, by inverting one class of data, we actually have 60 points available for estimating the prototypical manifold. In this way, we obtain even more evidence that the test point is from class B.
One way to endow methods with the ability to extract global structure is to use parameterized models that can generalize in nontrivial ways. In Bayesian terms, we have prior expectations about certain properties of the data. For example, we expect the probability density function for the data within a given class to be smooth on some scale. The class of distributions that our model can represent should reflect these prior expectations. By fitting the model to the data, this prior knowledge is then modified to obtain a more data-driven set of posterior expectations. In the above example, we decide that a sinusoidal manifold is a reasonable compromise between our prior expectations regarding continuity and the observed data within class B.
Simple parametric models, such as multidimensional Gaussian density functions, can be used to obtain some degree of generalization. However, overly simple models of this sort are inflexible in that they cannot generalize in complex ways. Also, for real-world problems, such inflexible models are often


Page 57
Figure 3.2
Two classes of 2-dimensional training
data and a test point.
inconsistent, since they often cannot represent the complexity in natural data sets. In this chapter, I examine the use of more flexible Bayesian network models for pattern classification.
I begin this chapter with a description of how Bayesian networks can be used for pattern classification. Then, in Section 3.2, I present the "autoregressive" network which is quite simple, but performs surprisingly well as a pattern classifier. In Section 3.3, I describe maximum likelihood estimation and "maximum likelihood-bound" estimation for models with latent (hidden) variables. Latent variables are not part of the input pattern, but are meant to represent higher-order structure in the data (e.g., handwriting style). In Section 3.4, I review three techniques for estimating the parameters of sigmoidal Bayesian networks with latent variables: Gibbs sampling, variational inference, and the wake-sleep algorithm. Then, in Section 3.5, all of these models are compared with the k-nearest neighbor classifier and a tree-based classifier when classifying handwritten digits.
An area which is closely related to estimating probability models for pattern classification is unsupervised learning. I view unsupervised learning as the process of estimating a probability model for a class of data. The hope is that some of the latent variables in the model will come to represent interesting features, and that these features can then be automatically extracted for novel input patterns. In Section 4.1 I present results for the Helmholtz machine, when it is given the task of trying to extract structure from noisy 16 ´ 16 images of horizontal and vertical bars. Finally, in Section 4.2, I present a new type of parameterized Bayesian network that can be used to simultaneously extract continuous and categorical structure in an unsupervised manner.


Page 58
3.1 Bayesian networks for pattern classification
Bayesian networks provide a means of producing structured probabilistic models with arbitrary complexity. In this sense, they are flexible models. The majority of this chapter is devoted to using Bayesian networks to produce one model for each class of training data. A new test pattern is classified by choosing the class of the model that is best suited to the test pattern. In contrast, it is certainly possible to construct a Bayesian network that has one set of pattern variables v, a variable that represents the class j, plus other variables that represent important physical effects. An inference method can then be used to compute P( j½v) using the network. An advantage of this approach is that the model may make efficient use of the similarities and differences between all of the classes. For example, if each class of data in Figure 3.2 is modelled separately, then the similarity between the two classes cannot be exploited as described above. In practice, however, a parameter estimation algorithm may fail to find such similarities and in the process of trying to model both classes fail to properly extract the features from any one class. Another disadvantage of the single-model approach is that a new class of data cannot be introduced without refitting the model. Despite these disadvantages, the single-model approach is seductively interesting. In Chapter 4, I study networks that are estimated from unlabeled data, where the hidden variables automatically come to represent data classes.
The multiple-model approach to pattern classification consists of estimating one model for each of the J classes of data. In this sense, each model is conditioned on a class number. For the sake of generality, I will assume that the jth model has a set of features or hidden attributes hj that help model the pattern variables v. Network j thus represents a distribution P(v, hj½j). Finally, the class probabilities P( j) must be specified; these are simply determined from the relative sizes of the classes of data and any prior knowledge we have at hand. (For example, even though a training set contains 10 patterns from class 0 and 14 patterns from class 1, if we know ahead of time that the classes are equally likely then we set P( j=0)=P( j=1).)
Ideally, the model estimate for class j will yield a distribution that is close to the true distribution Pr(v½j) of the data from class j:
However, even if the approximation is good, the sum in the above expression is exponential in the number of feature variables, and so cannot be computed


Page 59
directly. Instead, for a given test pattern v, one of the inference methods described in Chapter 1 can be used to produce class likelihood estimates,
Finally, Bayes' rule is used to produce soft classification decisions,
and a hard decision j* can be made by choosing the best class,
The technique used to estimate the class models and the inference method used to estimate P( j½v) depend on the structure of the networks. Before examining intractable models for which inference and parameter estimation must be approximated, I discuss an interesting class of tractable systems. For the sake of notational simplicity, the following sections present models and algorithms for estimating P(v), with the class index j left off. It should be kept in mind that one such density model must be estimated for each class.
3.2 Autoregressive networks
There are a variety of Bayesian network architectures for which inference and parameter estimation can be performed exactly within a reasonable amount of time. An architecture of this type that I discuss here is easy to implement and works surprisingly well on some problems. I define an autoregressive network as a fully-connected parameterized Bayesian network without any latent variables. The graph for the network is thus specified completely by an ancestral ordering. Unless I am considering different ancestral orderings of the variables, I will usually assume that the variables are labeled in the ancestral order. Then, the parameterized distribution P(v½q) for the data can be written
where q is the entire set of parameters, and q i is the set of parameters associated with input v i. Each of the conditional probability distributions in this expression is represented using some sort of parametric or flexible model. Figure 3.3 shows an example of an autoregressive network with five variables.


Page 60
Figure 3.3
An autoregressive network with
ancestral ordering v1, v2, v 3, v 4, v5.
3.2.1 The logistic autoregressive network
If the pattern consists of binary variables (v i Î {0, 1}) logistic regression [McCullagh and Nelder 1983] (see Section 1.2.7) may be used:
where g(x) = 1/(1 + exp[-x]) is the logistic function, and a dummy variable v0 = 1 is used to account for a constant in the arguments of the exponent.
For this logistic autoregressive network, P(v|q) can be computed in O(N2) time in the following way. For each variable v i, the sum [image] is determined from the values of vl, . . . , v i-1, and then [image] is determined from the value of vi using (3.6).P(v|q) is then computed using (3.5).
3.2.2 MAP estimation for autoregressive networks
An autoregressive network can be fit to a set of training patterns v(1). . . , v(T ) using MAP parameter estimation. To do so, we need to specify both a prior distribution over the parameters q, and also the training set likelihood given the parameters. Assuming that each training case is independent and identically drawn (i.i.d.), the log-likelihood of the training set is


Page 61
If the parameters are independent under the prior, then given a training set, the ith term in the sum of the last expression in (3.7) depends on a set of parameters qi = {q i l, . . . ,qi,i-1} that are independent of all the other sets of parameters q i', i' ¹ i. So, MAP estimation can be broken down into N subproblems, where subproblem i is to estimate the parameters q i for the model that predicts vi from
Here, I derive a gradient-based MAP estimation method for the logistic regression model used in the logistic autoregressive network. Let the data for subproblem i be denoted [image]. Up to a constant of proportionality that does not depend on qi, the likelihood of the training data for subproblem i is
where the last expression is obtained from (3.6).
I use a prior distribution under which the parameters in q i are independent and normally distributed with mean 0 and a fixed variance [image]:
Up to a constant of proportionality that is independent of q i, the posterior distribution over the model parameters q i given the data [image] is
Taking the logarithm of this expression and leaving out constants that do not affect the optimization procedure, MAP estimation for subproblem i entails maximizing


Page 62
I use the conjugate gradient optimization method [Fletcher 1987] which requires the derivatives of [image] :
Both [image] and its derivatives can be computed in O(iT) time.
3.2.3 Scaled priors in logistic autoregressive networks
In the prior distribution over the parameters (3.9), how should the variance [image] of the parameters for the ith input depend on i? That is, before seeing any training data, how do we expect the variance of the parameters for a variable to depend on how many inputs that variable receives?
Assume we don't have prior knowledge of a preferred ordering of the variables. By symmetry, it makes sense to assume a uniform prior distribution over the variables; i.e., under the prior each variable is equally likely to have each of the values 0 and 1. The dummy variable v0 = 1 is exempt from this prior, of course. Now, consider the prior probability predictions made for vi. This prior distribution has two sources of variability: a Gaussian prior over the parameters q i, and a uniform distribution over the inputs [image] . By symmetry, this prior distribution over the probability predictions made for vi should not depend on i. As shown below, this restriction determines how to set the variance for the parameters q i for each variable v i.
Since the prediction for vi is determined by its total input [image] , the the above restriction can be applied to the total inputs. Averaging over the two sources of variability, we get a mean value of the total input for v i of
We can take E[q ikvk] = E[q ik]E[v k], since the parameters and the inputs are independent under the prior. The final step holds since the parameters have mean 0.


Page 63
Since the mean total input is 0, the variance of the total input for vi is
Now, since q ik and q ij, k ¹ j are independent under the prior, E[q ikq ij] is nonzero only if j = k. So, the variance of the total input for v i is
Under the prior, the probability predictions for vi should not depend on i. So, the variances of the total inputs for v 1 and vi should not differ:
Note that [image] is the variance of the total input for v1, which has no input variables. All of the variances can be set by picking a reasonable value for [image] . In my simulations, I chose [image]. This value allows for probabilities near 0 and 1 at the output of the logistic function, without favoring them too much (see Figure 1.10 on page 22).
It may be a good idea to let the biases in the network have a separate Gaussian prior, although I have not yet explored this possibility experimentally.
3.2.4 Ensembles of autoregressive networks
An autoregressive network is specified by chosing an order for the variables v1, . . . , v N . Leaving computational considerations aside, if the subproblem models [image] are consistent (i.e., they converge to the correct distribution as the number of training examples tends to infinity) and there is a sufficiently large training set, then the particular ordering chosen is not important. The model for subproblem i will correctly represent the real conditional distribution [image], and so the product of the subproblem distributions will give the true joint distribution. However, the data sets considered here are small, and the parametric subproblem models considered here (e.g., logistic regression) are inconsistent for many distributions of data. In this case, the order of the variables is important in two contrasting ways. Certain orderings may give rise to simpler true conditional distributions [image] that can be more accurately represented by the model distributions [image] . In contrast, for a given training set, different orderings may lead to different amounts of overfitting.


Page 64
I do not address here the difficult issue of how to select an ordering that optimally balances these two effects. This problem is difficult both because the discrete ordering cannot be optimized by a gradient-based method and because for the training sets I will consider here, there is not enough data available to get a reliable estimate of which ordering is best. Instead of searching for an optimal ordering, I estimate an ensemble of autoregressive networks, where each network uses a randomly selected ordering of the variables. The probability prediction for a given vector v is then taken to be the average of the predictions over the ensemble of networks.
3.3 Estimating latent variable models using the EM algorithm
The notion of unobserved or hidden variables arises in many model estimation contexts. For example, due to mechanical failure, training data derived from physical measurements may sometimes lack values for some variables in some cases. In contrast, it is often useful to build hidden variables into a model by design. These variables are meant to represent latent causes that influence the visible variables. Several of the Bayesian network models discussed in the remainder of this chapter are latent variable models (e.g., see Section 3.4). For the sake of notational simplicity, I will use v to refer to the observed variables and h to refer to the unobserved variables. This is a slight abuse of notation, since it can happen that some visible variables are unobserved. For example, several of the photo-sensors in a digital camera may be burned out, so that some of the variables in the image pattern v are unobserved.
We would like to estimate a probabilistic model P(z) for a training set consisting of T patterns v(1), v(2), . . . , v(T ), where each pattern specifies the values of an observed subset v of the variables in z. In general, each training case may specify a different subset of visible variables.
Let h(t) = z \ v(t) be the set of hidden (unobserved) variables for training case t. Assuming that the training cases are i.i.d., the log-likelihood of the training data is


Page 65
To maximize this log-likelihood, we set its derivative with respect to each parameter q in q to zero:
The relation [image] was used in the first and third line of the derivation. Even though [image] log [image] is quite often easy to compute, in many cases of practical interest the system of equations obtained by setting [image] log [image] to zero for each q is highly nonlinear and cannot be solved in closed-form. One approach is to perform gradient descent in log P(v(t), h(t)|q ), while sampling from P(v(t), h(t)|q ) using Markov chain Monte Carlo. This gives a Monte Carlo approximation to gradient descent in log [image] as given in (3.18). Another approach is to solve the system of nonlinear equations iteratively. Although in principle any method for solving a nonlinear system of equations can be used (e.g., Newton's method [Fletcher 1987]), the structure of (3.18) gives rise to a particularly simple two-phase iterative method, called the expectation maximization (EM) algorithm [Baum and Petrie 1966; Dempster, Laird and Rubin 1977].
3.3.1 The expectation maximization (EM) algorithm
Often, we have available an efficient method for estimating the model when all of the variables are visible. That is, the system of equations obtained by setting


Page 66
for arbitrary h(t) can be solved quite easily. Notice that it is essentially this system of equations that is obtained if the dependence of P(h(t)|v(t),q ) on q in (3.18) is ignored. The summation over h(t) in (3.18) has the effect of replicating training case t once for each configuration of the hidden variables for that case, and weighting each replication by P(h(t)|v(t), q ). This observation leads to the following iterative two-phase EM algorithm:
1. E-step: Compute P(h(t)|v(t),q )for each configuration h(t) of the hidden variables for each training case, and set Q(h(t)) ¬ P(h(t)|v(t), q ).
2. M-step: Solve the following system of equations for q ,
Stop if a convergence criterion is satisfied; otherwise go to 1.
In practice, the values of Q(h(t)) for each training case are usually not stored during the E-step. Instead, statistics that are sufficient for the M-step are accumulated while processing the training set. There are several proofs that each EM iteration is guaranteed to increase the likelihood of the training data [Baum and Petrie 1966; Dempster, Laird and Rubin 1977; Meng and Rubin 1992; Neal and Hinton 1993]. After presenting a more general algorithm for maximizing lower bounds on the data likelihood [image] , I will show that each iteration of EM is guaranteed to increase the data likelihood.
3.3.2 The generalized expectation maximization algorithm
Neal and Hinton [1993] introduced a new view of the EM algorithm as a method for maximizing a lower bound on the likelihood of a training set. This interpretation opened the door to tractable approximations to EM for models that are clearly intractable. This ''generalized expectation maximization algorithm'' is an approximation to ML estimation that follows from using the wrong distribution Q(h(t)) in the M-step of the EM algorithm; i.e., Q(h(t)) ¹ P(h(t)|v(t), q). There are practical reasons for using a suboptimal distribution Q(h(t)), the most obvious being that in some cases it is computationally infeasible to compute Q(h(t)) for every configuration of the hidden variables h(t) for each training case. For example, some of the Bayesian network models discussed below have over one million configurations per training case.


Page 67
The bound used in generalized EM is obtained using the following form of Jensen's inequality [Cover and Thomas 1991]:
where Si q i = 1, and ai are arbitrary scalars. Applying this inequality to the log-likelihood of the training data (3.17), we get
The goal of generalized EM is to jointly estimate a distribution Q(h(t)) (which may or may not be parameterized) and a distribution P(v(t), h(t)|q ), so as to maximize this lower bound on the likelihood. This leads to the following algorithm:
1. Generalized E-step: Increase the bound B Q||P with respect to a distribution Q(h(t)).
2. Generalized M-step: Increase the bound B Q||P with respect to q .
Note that unlike the E-step of the EM algorithm, the generalized E-step may produce a Q-distribution for which Q(h(t)) ¹ P(h(t)|v(t), q ).
The EM algorithm can be viewed as a special case of the generalized EM algorithm, where we alternately maximize the bound B Q||P with respect to an unconstrained distribution Q(h(t)), and then with respect to P(v(t), h(t)| q ) via q . If the bound is maximized with respect to Q(h(t)) during the generalized E-step, while enforcing S h(t) Q(h(t)) = 1 using a Lagrange multiplier, we obtain Q(h(t)) = P(h(t)|v(t), ). This form of the generalized EM algorithm is identical to the standard EM algorithm presented in the previous section. Also, in this case the inequality in (3.22) becomes an equality: . It follows that the EM algorithm is a maximum likelihood estimation method.


Page 68
Note that in general, generalized EM does not give the same estimates as ML estimation. As a degenerate example, imagine that we use ML estimation to obtain a model P(v(t), h(t)|q ) from a training set, and that we then apply generalized EM with Q(h(t)) fixed at a uniform distribution. In this case, the bound can be increased by moving P(v(t), h(t)|q ) away from the ML estimate (unless P(v(t)|h(t),q ) happens to be uniform, in which case a uniform Q(h(t)) makes the bound tight so that P(v(t), h(t)|q ) will not change). However, as long as we are able produce estimates of Q(h(t)) that are "close" to P(v(t), h(t)|q ), generalized EM will be close to ML estimation. Of course, in most cases, if we have the computational resources available to obtain an ML estimate, generalized EM should not be used. In Section 3.4, I introduce a class of Bayesian networks that have many latent (hidden) variables. For these networks, it is computationally intractable to perform EM and so generalized EM is used.
3.4 Multiple-cause networks
In many cases, it makes sense to postulate that a data vector v naturally arises from the consequences of a set of hidden causes h. For example, an image may be nicely described as a two-dimensional rendition of a combination of objects. If h k is a binary variable indicating the presence of object k in the image, then the model distribution P(v|h, q V) is the distribution over images given which objects are present. (q V is a set of parameters associated with the distribution over v). This distribution is meant to capture the way in which the objects interact to form the image as well as any inexplicable noise.
The model P(v|h, q V) may be simplified by assuming that the K causes dependency-separate the image pixels. That is, once we know which causes are present, each pixel is independent of the others. In this case, [image] [image] . If the visible variables are binary, each conditional distribution can be implemented using, for example, logistic regression. In contrast to the logistic autoregressive network where each visible variable is regressed on a subset of the other visible variables (see (3.5)), in the multiple-cause network each visible variable is regressed on the hidden cause variables h:
where [image] , and we take h0 = 1 in order to account for a constant in the summations. Binary Bayesian networks which use logistic regression for the conditional distributions are often called binary sigmoidal


Page 69
Figure 3.4
A multiple-cause network with five visible
variables v and four hidden cause variables h.
networks [Neal 1992] and are sometimes called stochastic multi-layer perceptrons.
To complete the model, we provide a distribution P(h|q H) over the set of causes. Although it seems natural that in many cases the hidden variables h k might be interdependent, for the sake of simplicity, I will assume for now that they are not:
where [image]. These probabilities can be nicely parameterized using the logistic function:
An example of this type of multiple-cause Bayesian network is shown in Figure 3.4. Notice that the dependency-separation of the variables v by the set of hidden variables h is ensured by condition 2 described in Section 1.2.5.
Supposing that we have somehow obtained an accurate model of the true causal process for each class of data (e.g., using a method described below), in order to perform classification we would like to compute the marginal probability P(v|q ) for each class model. This can be computed exactly using
where q = {q H, q V} is the entire set of parameters. However, this sum is exponential in the number of causes K and so in practice, we must use another approach. It is obvious from Figure 3.4 that probability propagation cannot


Page 70
be used to obtain an exact result, since the Bayesian network contains many cycles. In fact, we must use an approximate inference method.
Since exact probabilistic inference is needed for ML and MAP parameter estimation, these estimation methods are also intractable. For example, in order to perform the E-step of the EM algorithm, we must compute P(h|v,q ), which has an exponential number (2K) of terms.
In the next three sections, I show how Gibbs sampling, variational inference, and the stochastic Helmholtz machine can be used to approximate P(v|q ) and perform generalized EM parameter estimation in multiple-cause networks.
3.4.1 Estimation by iterative probability propagation
For generalized EM by probability propagation, the Q-distribution in the likelihood bound (3.22) can be computed by applying probability propagation in the multiply-connected network. (Since probability propagation is only exact in singly-connected graphical models, propagation in the loopy network will give approximate inference.) Due to the loops in the network, probability propagation will not naturally terminate, but will instead lead to an iterative algorithm. Initially, I frowned upon this approach because modifying the model parameters may cause intermediate stages of iterative inference to become unstable. However, it turns out that variational inference is also susceptible to this effect, so in fact estimation by iterative probability propagation may be an effective learning procedure.
3.4.2 Estimation by Gibbs sampling
In order to estimate a multiple-cause network from a set of training examples v(1), . . . , v(T), we can perform on-line gradient descent in log P(v,h|q ) while sampling from P(h|v,q ) using the Gibbs sampling method described in Section 2.2.2. For the current training case v(t), we simulate a Markov chain to obtain a configuration h(t) of the hidden variables. (Notice that in general h(t) will be different each time [image] is processed—ideally, h(t) will have a distribution P(h|v(t),q ).) In order to perform Gibbs sampling for the logistic multiple-cause network described above, we need to be able to sample from the distribution for each hidden variable hk given the other cause variables and the visible variables. Since the cause variables are binary, we only need to compute a function that is proportional to [image] . The two values can then be normalized to obtain [image] . Since the total joint probability for h and v can be easily computed in O(KN) time using the ancestral ordering, the joint probability can be used to compute the conditional distribution as follows:


Page 71
For a given training case v(t), the latent variables are visited in a specified order while drawing a new value for each variable from its conditional distribution. The entire set of latent variables h is processed in this fashion for a specified number of times before the Markov chain is terminated and some configuration h(t) of the latent variables is produced.
The biases for the hidden variables and the parameters connecting the hidden variables to the visible variables are adjusted by following the derivatives of log P(v(t), h(t)|q using a learning rate [image] as follows:
In order to perform classification, we would like to compute P(v|q) for an input vector:
This problem can be viewed as a form of free energy estimation [Sheykhet and Simkin 1990; Neal 1993]. I use a very simple approximation that is quite fast and works well in practice for classification purposes. Since the number of terms in the above sum is exponential in the number of causes, I approximate it by assuming that the majority of the total probability mass is contributed by a small number of clusters in h-space. These clusters are found by simulating a Markov chain as described above. At multiple points in the chain, the configuration of h and all neighboring configurations (i.e., those configurations within a Hamming distance of 1) are added to a list of "significant terms". Only the neighboring states of h are considered because once P(v, h|q ) has been computed, it is easy to compute the probabilities for configurations that differ from h by only one bit. After a specified number of clusters have been visited in this manner, the above sum is approximated by adding up the terms for the tabulated configurations. This method for estimating P(v|q ) will not work well when there is a large number of clusters in h-space that contribute significantly to the sum. However, I have found that in practice the Gibbs sampling learning algorithm tends to favor a small number of clusters, making this approximation reasonable.


Page 72
3.4.3 Generalized EM using variational inference
In this section, I review the variational method developed by Saul et al. [ 1996] for generalized EM in sigmoidal Bayesian networks. It turns out that a product-form variational distribution leads to an intractable bound, and so the bound itself must be bounded by a tractable function.
For generalized EM by variational inference, the Q-distribution in the likelihood bound (3.22) depends on some variational parameters x. For the sake of simplicity, consider the bound for one training case v:
Generalized EM entails iteratively maximizing this bound, first by varying x (the generalized E-step), and second by adjusting the model parameters q (the generalized M-step). The first term in this bound is the entropy of the variational distribution Q(h|x), and the second term is the expected log-probability of h and v under the variational distribution.
Here, I consider a product-form variational distribution over the K latent variables h1, . . . h k (notice that h0 is not included since it is fixed to h0 = 1):
where xk is the probability under the variational distribution that h k = 1. Using this variational distribution, the first term in B Q||P simplifies to
The second term in (3.31) is


Page 73
Since the conditional probabilities are given by logistic regression, this term contains many expectations of nonlinear functions. The first step to simplifying these expectations is to express the conditional probability [image] given in (3.25) in the following way:
The expectation of log [image] is
Similarly, the conditional probability [image] can be written
The expectation of log [image] is then
The overall bound for v is
Except for the last term, the values of these terms and their derivatives can quite easily be computed. The explicit summation over h in the last term cannot be reduced to a tractable form. However, the last term can be bounded by introducing some extra variational parameters v. (See Saul, Jaakkola and Jordan [1996] for details.) Generalized EM for the new bound B'Q ||P £ B Q ||P entails iteratively maximizing this bound, first by varying x and v (the generalized E-step), and second by adjusting the model parameters q (the generalized M-step).


Page 74
3.4.4 The stochastic Helmholtz machine
A stochastic Helmholtz machine consists of a pair of Bayesian networks that are fit to training data using an algorithm that approximates generalized EM, where the bound on the likelihood may be very complex. In addition to the multiple-cause network that describes P(v, h|q ) (the generative network), there is a recognition network that describes Q(h|v,f). The stochastic Helmholtz machine requires that the recognition network have an ancestral ordering such that it is easy to draw samples from Q(h|v,f). The advantage of the stochastic Helmholtz machine over Markov chain Monte Carlo is that each sample from the recognition network is independent, as opposed to dependent on the last sample. The advantage of the stochastic Helmholtz machine over variational inference is that more complicated (e.g., nonfactorial) distributions can be represented by the inference process used for generalized EM. The main disadvantage of the stochastic Helmholtz machine is that a recognition network that is compatible with the generative network must somehow be estimated, and this can be a very difficult task when a complex recognition network is used. An example of inference in the stochastic Helmholtz machine is described in Section 2.4.4. Here, I describe the wake-sleep algorithm for on-line estimation of both the generative and recognition parameters (q and f ) [Hinton et al. 1995].
Suppose we have a current generative network (which may or may not be a good model of the data) and a current recognition network. For a parameterized recognition network, the likelihood bound in (3.22) is
This bound can be estimated by averaging log P(h,v(t)| q /Q(h,v(t), f ) over multiple recognition sweeps for each input pattern. In each recognition sweep, the recognition network is stochastically simulated to obtain a configuration h of the latent variables.
We would like to maximize B Q ||P with respect to the recognition network parameters f for all v, if possible. As discussed in Section 3.3.2, the unconstrained recognition distribution that maximizes this likelihood bound is
However, except for very simple recognition networks, this optimization is intractable for the same reason that exact inference is intractable. Instead, we optimize a different function whose global maxima give identical recognition


Page 75
networks in certain limits to those produced by maximizing B Q ||P. The limits may not apply in practice, so that the recognition network may be slightly suboptimal.
Assume for the moment that the recognition network is consistent with the distribution P(h|v,q ). In this case, the parameters f that maximize
will also maximize B Q ||P in (3.40). (Note the reversed order of the distributions). So, for a given generative network, the optimum recognition network can be found by maximizing B P ||Q with respect to the recognition parameters f. The derivative of B P ||Q with respect to a recognition network parameter f is
So, the recognition network can be estimated using stochastic gradient descent by sampling h and v from P(h|v,q ) using ancestral simulation, and then adjusting the recognition network parameters so as to increase the log-probability of the hidden variables given the visible variables. This procedure is called sleep-phase learning, since the recognition network is adjusted to be better suited to the "fantasies" produced by the generative network.
In practice, sleep-phase learning is only an approximation to the generalized E-step described in (Section 3.3.2) for one main reason. An ideal recognition network produces a good approximation to P(h|v,q ) even for a vector v that has a very small probability under P(h|v,q ). (This corresponds to a plausible real-world pattern that the generative network has not yet learned.) For sleep-phase learning to produce such a recognition network, an extremely large sample size must be drawn from P(h|v,q ) in order to get an example of the unlikely vector. For the sake of tractability, a relatively small sample size must be used, which implies that the ideal recognition network cannot be found. This means that in practice, maximizing B P ||Q does not give the same recognition network as would be obtained by maximizing B Q ||P . In fact, in order to prevent overfitting of the recognition network, an inconsistent parametric recognition network is used, so that the global maxima of the two functions may not even coincide.
For a given recognition network, the generative network is adjusted in the wake-phase using a Monte Carlo implementation of the generalized M-step. That is, on-line stochastic gradient descent in the likelihood bound B Q ||P is


Page 76
performed with respect to the generative network parameters q. The derivative of the bound with respect to a generative network parameter q is
For a training vector v(t), the recognition network is used to sample values for the latent variables h. Then, the generative parameters are adjusted so as to increase the log-probability of the latent variables and the visible variables.
The two phases of learning are usually applied in alternation. A training pattern is presented; the recognition network is used to obtain a random h; and the generative network is adjusted. Next, the generative network is used to obtain a random h and then a random v; and the recognition network is adjusted. The result of this constant mixing of the two phases is that the generative network becomes better at modelling the training data and at the same time tries to produce causes for the training data that can be properly inferred by the restricted recognition network. This can be seen mathematically by breaking B Q ||P into two pieces:
The first term encourages the generative network to model the data, whereas the second term (a negative Kullback-Leibler pseudo-distance) encourages it to be compatible with the recognition network. As a result of the latter, for a generative network that is estimated using the wake-sleep algorithm, the global maxima of B Q ||P and B P ||Q often do coincide.
Assuming that under the recognition distribution, the latent variables are independent given the visible variables, we have:
Also, consider modelling each of these components using logistic regression:


Page 77
Figure 3.5
A recognition network that implements Q(h|v,f )
for the generative network shown in Figure 3.4.
where we take v0 = 1 in order to account for a constant in the exponents. This recognition network is shown in Figure 3.5.
For this logistic recognition network, the recognition parameters are adjusted as follows during the sleep phase, in order to increase log Q(h|v,f):
It turns out that in many practical cases this recognition network is sufficient for producing good density estimates. However, if it is estimated in conjunction with a fixed generative network that describes the simple burglar alarm problem (see Section 2.2.3), the likelihood bound B Q ||P may actually decrease. Consider how the recognition network is modified for fantasies where the burglar alarm is ringing. We simulate the generative network, obtaining values for b, e, and a, and discard those samples for which a ¹ 1. For the recognition network parameters f BA and f EA that connect the common consequence a to the two causes b and e, the expected learning updates become
where f B 0 and f E 0 are the recognition biases for b and e. Each connection is modified so as to predict as closely as possible the marginal posterior distributions P(b = 1|a = 1) and P(e = 1|a = 1) over the corresponding causes b and e. After training, the recognition distribution over b and e given a = 1 will


Page 78
be the product of the marginals. For the configuration b = 1, a = 1,
where the values for P(b = 1|a = 1) and P(e = 1|a = 1) were computed from (2.8). This value is quite a bit higher than the correct value of P(b = 1,e = 1|a = 1) = 0.116. In fact, if we assume that b and e are independent given a = 1, the recognition distribution that maximizes the likelihood bound B Q ||P has [image]. This is an example where maximizing B P ||Q is a poor approximation to maximizing B Q ||P. Notice, however, that the problem arises because we are using an inconsistent recognition network.
3.4.5 Hierarchical networks
Earlier in this section, I presented multiple-cause networks with the assumption that the causes were independent in the generative network (see (3.24) on page 69). Even with the assumption that the causes are independent, it is still possible to represent quite complex correlations in the visible variables v. However, in many cases the causes are certainly interdependent. For example, if the causes variables represent the presence or absence of various objects in facial images, we expect that both a toque and a top-hat are not present simultaneously.
In order to model interdependent causes, we can simply add another layer of ''meta-causes'' at a higher level in the network. Even if we assume that the "meta-causes" are independent, the network can still represent fairly complex relationships between the causes. Such a hierarchical network is shown in Figure 3.6a. We have already seen examples of hierarchical networks, such as the network used in the bars problem example (Figure 1.11).
The parameter estimation methods already described in this chapter are applicable to hierarchical networks with any reasonable depth. Of particular interest, however, is the recognition model for the Helmholtz machine. Figure 3.6b shows a recognition model that is appropriate for the network in Figure 3.6a. Note that the top layer of hidden variables could receive input from the bottom layer of visible variables, not just the middle layer of hidden variables. However, this introduces extra parameters into the Helmholtz machine, which may worsen the effect of overfitting. In my experiments, I use layered generative and recognition networks like the ones shown in Figure 3.6. In some cases, adding extra connections may help.


Page 79
Figure 3.6
(a) A hierarchical network with three extra "meta-cause" variables which
produce correlations between the cause variables (middle layer).
(b) A recognition network.
3.4.6 Ensembles of networks
According to the Bayesian doctrine for prediction, when using a density model P(v|q ) to estimate the probability of v, we ought to integrate out the model parameters q . For a prior parameter q probability density p(q ) and a training data likelihood [image] , we use the posterior [image] (where the constant of proportionality does not depend on q ) to obtain an estimate
For example, this integral can be approximated using Laplace's approximation [Spiegelhalter and Lauritzen 1990], Markov chain Monte Carlo methods [Neal 1993; Neal 1996] or variational techniques [Jaakkola and Jordan 1997]. Here, I consider an ensemble method that is less sophisticated than the above approaches, but is also easier to implement and and in practice usually gives a significant improvement over MAP parameter estimation.
Suppose we perform MAP parameter estimation using multiple restarts (different random initial parameters) so that we have an ensemble of M models, where model m has parameters q m. Each model may correspond to a different local maximum of the posterior [image], and we assume that each model is equally likely in the posterior. We then approximate the above integral by


Page 80
If P(v|q ) does not change much over the width of each mode in the posterior, then as long as the modes are properly represented by the ensemble of models, [image] will be very close to the correct value [image] given by integration. On the other hand, if there is a mode in the posterior that is so wide that P(v|q ) does vary significantly across the mode, then [image] may be quite different from [image] . This is because only a peak in the posterior is being included in the sum, while the mass surrounding the peak is being ignored, even though the corresponding predictions are quite variable.
Since the Bayesian networks described above are flexible models, we expect that with limited training data they may have multiple data likelihood optima (corresponding to multiple peaks in the posterior, if we assume a uniform prior over network parameters). For this reason, when time permits, a significant classification rate improvement can be obtained by using an ensemble of networks for each class of data. The classification decision is then based on the average probabilities computed from the ensemble for each class.
3.5 Classification of handwritten digits
An interesting and useful pattern classification problem is the classification of handwritten digits. In this section, I present results on the classification of 8 ´ 8 binary images of handwritten digits made available by the US Postal Service Office of Advanced Technology. I compare the following Bayesian network methods: logistic autoregressive classifier (LARC-1), a stochastic Helmholtz machine with one hidden layer (SHM-1), a stochastic Helmholtz machine with two hidden layers (SHM-2), and an ensemble of stochastic Helmholtz machines with one hidden layer (ESHM-1). In order to place the performance of these networks in context, I include the following methods: classification and regression trees (CART-1), the naive Bayes classifier (NBAYESC-1), and the k-nearest neighbor method (KNN-CLASS-1). The performances of these classifiers are assessed using 5 different training set sizes (120, 240, 480, 960 and 1920 cases) so that the effect of the number of training cases on each method can be studied. After describing the classifiers and the methods used to estimate them, I present and discuss the performance results.
3.5.1 Logistic autoregressive classifiers: LARC-1,ELARC-1
LARC-1 models each of the 10 classes of data using a logistic autoregressive network (see Section 3.2), where the variables are ordered in a raster-scan fashion. Once each of the 10 networks have been estimated from the training data, a test pattern is classified by outputing the class corresponding to the network that gives the greatest likelihood to the pattern.


Page 81
Before estimating each network from its respective class of training data, the double precision parameters q were initialized to uniformly random values on [-0.01, 0.01). Overfitting was prevented by using MAP estimation with a scaled Gaussian parameter prior. The prior variance of the first input was set to [image]. A conjugate gradient algorithm was used for MAP estimation.
ELARC-1 uses an ensemble of 8 logistic autoregressive networks, where each element in the ensemble uses a different ordering of the variables in v. One of the elements uses the raster-scan ordering, whereas the other 7 elements use a randomly selected ordering. The probability of a test pattern for a given class is estimated by averaging the probability estimates from each of the 8 networks in the ensemble for that class.
3.5.2 The Gibbs machine: GM-1
GM-1 models the distribution of each of the 10 classes of data using a logistic multiple-cause network of the type shown in Figures 3.4, that is trained using Gibbs sampling. Each network has 64 visible binary (0/1) variables (8 ´ 8) and one hidden layer of 16 binary (0/1) variables. Once the 10 networks have been estimated, classification of a test pattern proceeds by estimating the probability of the pattern under each network, using the method described in Section 3.4.2. The class corresponding to the network that gives the highest probability is output as the prediction.
Before estimating a network using Gibbs sampling, all of its double precision parameters were initialized to uniformly random values on [-0.01, 0.01). For each training pattern, a single configuration of the hidden variables was obtained by performing 10 sweeps of Gibbs sampling, while annealing the network from a temperature of 5.0 to 1.0 using a 1/ T schedule, where T is the sweep count. Then, the network parameters were adjusted using a learning rate of 0.01. For a training set of T patterns, a randomly chosen set of [image] cases were set aside as a "validation" set. By monitoring the probability estimate for this validation data, early-stopping was used to prevent overfitting. After every 10 epochs of learning (1 epoch = one sweep through the remaining [image] training cases), for each validation pattern, 10 sweeps of Gibbs sampling with annealing were performed as described above, and then 20 sweeps of Gibbs sampling at unity temperature were performed to obtain 20 configurations. Then, the probability of the validation pattern was estimated by computing the probability mass associated with each configuration and its 1-nearest neighbors. Each network was trained for a minimum of 100 epochs (the validation probability estimate was still computed every 10 epochs in this interval). Then, learning was stopped after the current epoch n, if the epoch nmax at which the maximum validation probability estimate occurred took place no less than n/3 epochs ago. Also, in order to terminate learning


Page 82
runs where the validation probability estimate continued to increase asymptotically towards a limit, a maximum of 2000 training epochs were performed. In summary, learning was stopped at epoch n if n ³ 2000 or if nmax £ 2n/3 and n ³ 100. (A similar early stopping technique has been used with regression models [Rasmussen 1996].)
3.5.3 The mean field Bayesian network: MFBN-1
MFBN- I models the distribution of each of the 10 classes of data using a logistic multiple-cause network of the type shown in Figures 3.4. Each network is fit using the variational technique described in Section 3.4.3. Each network has 64 visible binary (0/1) variables (8 ´ 8) and one hidden layer of 16 binary (0/1) variables. Once the 10 networks have been estimated, classification of a test pattern proceeds by computing the likelihood bound for each network, using the method described in Section 3.4.3. The class corresponding to the network that gives the highest bound is output as the prediction.
Before estimating a network using the variational method, all of its double precision parameters (q )were initialized to uniformly random values on [—0.01,0.01). For each training case, the variational parameters were initialized to uniformly random values on [—0.01, 0.01). The variational bound was increased at each generalized E-step using the following iterative method [Saul, Jaakkola and Jordan 1996]. After each iteration, if the bound did not increase by more than 1% then no more iterations were performed for the current training case. A maximum of 10 iterations was performed. These algorithm parameters were suggested by Jaakkola (personal communication). The variational bound was increased at each generalized M-step using batch gradient descent with a learning rate of 0.01.
The validation procedure used to train each network was identical to the one used for GM- 1, except that the variational bound was used instead of an estimate of the validation case probability. The validation bound was computed every 5 epochs. No fewer than 100 epochs were performed, and no more than 1000 epochs were performed.
3.5.4 Stochastic Helmholtz machines: SHM-1, SHM-2, ESHM-1
SHM-1 models the distribution of each of the 10 classes of data using a stochastic Helmholtz machine with 64 visible binary (0/1) variables (8 ´ 8) and one hidden layer of 16 binary (0/1) variables. The generative and recognition networks are of the form shown in Figures 3.4 and 3.5, and logistic regression is used to implement the conditional relationships. The likelihood bound for a given input pattern is estimated using 20 recognition sweeps. Once the 10 machines have been estimated, classification of a test pattern proceeds by es-


Page 83
timating the likelihood bound for each machine. The class corresponding to the machine that gives the highest likelihood bound estimate is output as the prediction.
Before estimating a Helmholtz machine using the wake-sleep algorithm, all of its double precision parameters (q and f ) were initialized to uniformly random values on [-0.01, 0.01). A learning rate of 0.01 was used for both phases of learning. The validation procedure used to train each machine was identical to the one used for GM-1, except that instead of obtaining an estimate of the validation set probability as described above, 20 epochs of recognition passes were performed on the validation set to obtain an estimate of the likelihood bound for the validation data.
SHM-2 is similar to SHM-1, except that it uses stochastic Helmholtz machines with a visible layer of 64 binary variables, a middle hidden layer of 16 binary variables, and a top hidden layer of 8 binary variables. The generative and recognition networks are of the form shown in Figure 3.6.
ESHM-1 uses an ensemble of 8 SHM-1 networks to model each class of patterns. Each network in an ensemble is estimated using the above procedure, where a different randomly chosen validation set of [image] patterns is set aside for each network. Also, different initial random parameters are chosen for each network. Once 8 networks have been estimated for each of the 10 classes of data, a test pattern is processed by approximating 8 likelihood bounds for each data class. These are averaged together within each class to obtain 10 average likelihood bounds. The final class decision for the test pattern is based on these averages.
3.5.5 The classification and regression tree: CART-1
This tree-based classifier has previously been run on several classification tasks in DELVE [Rasmussen et al. 1996]. CART-1 uses a binary decision tree to classify the test patterns, where each node in the tree makes a binary decision based on an axis-aligned decision surface in the input space, and each leaf in the tree has a class label. A test pattern is classified by traversing the tree from the root to a leaf, while following the decisions at each node. That is, decision node d j looks at a particular input variable [image] and compares it to a threshold [image] the right child is chosen, and otherwise the left child is chosen. When a leaf is reached, the class of the leaf node is output by the classifier.
The tree is constructed from a training set using 10-fold cross validation. The details of how the tree is produced can be found in Breiman et al. [1984].1
1 I used Version 1.1 of the CART software provided by California Statistical Software Inc., 961 Yorkshire Ct. Lafayette, California 94549.


Page 84
3.5.6 The naive Bayes classifier: NBAYESC-1
The naive Bayes method of modelling can be viewed as a multiple-cause Bayesian network where there aren't any hidden cause variables. That is, we assume that each of the inputs is independent given the class identity. For the binary input case, this model becomes very simple. The naive Bayes model for each class of data is
where q i Î [0, 1] is the probability that v i = 1 under the model. For a given class of training data, I use the Bayesian method to obtain a minimum squared-loss estimate of P(v), assuming a uniform prior for q :
Once one such estimate is obtained for each of the 10 classes of training data, a test pattern is classified by choosing the class that gives the probability to the pattern. Notice that if ¦i = 0 or T , the probability P(v) is not 0 or 1. This prevents overfitting.
3.5.7 The k-nearest neighbor classifier: KNN-CLASS-1
This is the only nonparametric classifier studied in this section. The software I used was contributed to DELVE by Michael Revow [Rasmussen et al. 1996]. In order to guess the class of an input pattern v, the k-nearest neighbor classifier considers the classes of the k training patterns that are nearest to v in Euclidean distance. Let n j j = 0, . . . , J - 1 be the number of such training patterns in class j, so that [image] . Then, the k-nearest neighbor classifier outputs the most frequent class:
If two or more classes have the maximum number of k-nearest neighbor training patterns, then the classifier chooses the class whose training patterns are closest to v on average in Euclidean distance.
k is chosen using leave-one-out cross validation. If there are T training patterns, T new training sets with T - 1 patterns each are produced by leaving each pattern out once. k is set to 1, and the misclassification rate on the left out patterns is computed using the k-nearest neighbor classifier. Then, k is in-


Page 85
creased and this process is repeated until k = T - 1. The value for k that gives the lowest leave-one-out cross-validation error is used to make predictions for the test patterns.
In order to estimate the probability that v comes from each class, the k-nearest neighbor method uses
In this case, the squared difference between the predicted probability vector and the true class identity vector (a vectors of 0's with a single 1) is used as the cross-validation metric to determine k. See the DELVE manual [Rasmussen et al. 1996] for more information.
3.5.8 Results
The performances of the classification methods described above were assessed using the DELVE (data for evaluating learning in valid experiments) system [Rasmussen et al. 1996]. Under this system, the digit classification problem that I am interested in is called a prototask. A particular choice of training set size (e.g., 120 training patterns) is called a task. In order to get an accurate measure of the performances of the methods (with error bars), each method was trained and tested at least 4 times using disjoint training set - test set pairs (each of which is called a task instance). An original data set consisting of 10,960 patterns (1096 of each class) was partitioned into a training collection of 7680 patterns and a test collection of 3280 patterns. For each of the tasks with training set sizes 120, 240, 480, and 960, the training collection was partitioned into 8 disjoint training sets; for the task with training set size 1920, the training collection was partitioned into 4 disjoint training sets. Notice that for the tasks with training set sizes 960 and 1920, all of the training partition cases were used, whereas for the tasks with training set sizes 120, 240, and 480, not all of the training partition cases were used. For each of the tasks with training set sizes 120, 240, 480, and 960, the test collection was partitioned into 8 disjoint test sets with 410 patterns each; for the task with training set size 1920, the test collection was partitioned into 4 disjoint test sets with 820 patterns each. This way of partitioning the data eliminates the dependence between each of the 8 experiments performed to assess the performance of each method on each task instance.
Figure 3.7 shows the losses (fraction of patterns misclassified) for each of the tasks (five boxes). Each horizontal bar gives an estimate of the expected loss for a particular method on a particular task. The methods are ordered (from left to right within each box): CART-1, NBAYESC-1, KNN-CLASS1, MFBN-1, SHM-1, SHM-2, GM-1, ESHM-1, LARC-1—this is the same


Page 86
Figure 3.7
Estimates of expected fractions of misclassified patterns for nine
methods trained on five different sizes of training sets.
Table 3.1
Average time in minutes required to train and test the methods from Figure 3.7 for each of the training set sizes.
Training set size
Method
120
240
480
960
1920
CART-1
0.6
1.1
1.9
4.7
5.1
NBAYESC-1
0.0
0.0
0.0
0.0
0.0
KNN-CLASS-1
0.2
1.0
4.6
25.7
192.9
MFBN-1
19.8
64.4
130.5
216.9
344.6
SHM-1
2.3
5.4
11.8
21.8
46.0
SHM-2
3.7
7.8
22.3
41.7
77.6
GM-1
41.7
85.8
176.8
238.0
396.3
ESHM-1
19.4
44.7
95.5
186.7
358.9
LARC-1
0.2
0.5
1.2
3.0
6.5


Page 87
Figure 3.8
Estimates of expected negative log-probability of the
true digit identity given by each of eight methods trained
on five different sizes of training sets.
ordering as is given top to bottom in the lower-left hand region of the figure. Each vertical bar gives an estimate of the error (one standard deviation) for the corresponding estimate of the expected loss. Numbers in the boxes lying beneath the x-axis are p-values (in percent) for a paired t-test. Choose your favorite method from the list in the lower left-hand corner of the figure and scan from left to right. Whenever you see a number, that means that another method has performed better than your favorite method, with the given statistical significance. A low p-value indicates the difference in the misclassification rates is very significant. More precisely, a p-value is an estimate of the probability of obtaining a difference in performance that is equal to or greater than the observed difference, given that we assume the two methods actually perform equally well (the null hypothesis).
The ESHM-1 and LARC-1 methods clearly outperform all other methods for all tasks. If we scan the p-values for these two methods from left to right, we see that there is only a single method that performs significantly better than ESHM-1, and that is LARC-1 on the task with the smallest training set size (120). I found that the performance of ELARC-1 (ensemble of logistic autoregressive networks, not shown) was indistinguishable from plain LARC- 1


Page 88
with respect to classification error. In contrast, ESHM- I performs significantly better than SHM- 1. It is of particular interest that LARC- 1, which contains no latent variables, performs slightly better than the methods that contain latent variables.
GM-1 performs the best out of all approximate maximum likelihood methods, including MFBN-1, SHM-1, and SHM-2. However, GM-1 required an order of magnitude more training and validation time than SHM- 1. For this reason, an ensemble of logistic multiple cause networks was not considered for the Gibbs sampling estimation method. Table 3.1 shows the average time taken to train and test each method for each training set size on a 195 MHz MIPS R4400 processor. MFBN- I also required an order of magnitude more training and validation time than SHM-1, and so an ensemble of mean field Bayesian networks was not considered.
Figure 3.8 shows the performance results for soft decisions (the loss is the negative log-probability of the true class under each model). The methods are: NBAYESC-1, KNN-CLASS-1, MFBN-1, SHM-1, SHM-2, GM- , ESHM-1, LARC-1, ELARC- 1. In this case ELARC-1 performs significantly better than LARC-1. Also, in this case ESHM-1 performs slightly better than LARC-1 and ELARC- 1.


Page 89
4
Unsupervised Learning
The Bayesian networks described in the previous chapter are supervised, in the sense that they are trained with pattern - class label pairs. The task in unsupervised learning is for the learning machine to extract ''meaningful'' hidden causes for the sensory input. Perceptual tasks such as vision and speech recognition can be viewed as inference problems where the goal is to estimate the posterior distribution over the hidden causes (e.g., depth in stereo vision) given the sensory input (c.f. [Hinton and Sejnowski 1986; Hinton et al. 1995; Dayan et al. 1995]). The recent flurry of research in independent component analysis [Comon, Jutten and Herault 1991; Bell and Sejnowski 1995; Amari, Cichocki and Yang 1996; MacKay 1997b] exemplifies the importance of inferring the continuous-valued latent variables of input data. The latent variables found by this method are linearly related to the input, but perception requires nonlinear inferences such as decision-making and nonlinear continuous inferences, such as depth estimation. In this chapter, I show how the wake-sleep algorithm, slice sampling, and variational techniques can be used to learn networks of binary units and nonlinear continuous units in an unsupervised fashion.
4.1 Extracting structure from images using the wake-sleep algorithm
An interesting problem relevant to vision is that of extracting independent horizontal and vertical bars from an image [Foldiak 1990; Saund 1995; Zemel 1993; Dayan and Zemel 1995; Hinton et al. 1995]. Figure 4.1 shows 48 examples of the binary images I am interested in. Each image is produced by randomly choosing between horizontal and vertical orientations with equal probability. Then, each of the 16 possible bars of the chosen orientation is independently instantiated with probability 0.25. Finally, additive noise is introduced by randomly turning on with a probability of 0.25 each pixel that was previously off. The production of these images involves three levels of hierarchy: the first and lowest level represents pixel noise, the second represents bars that consist of groups of 16 pixels each, and the third represents the overall orientation of the bars in the image.


Page 90
Figure 4.1
Examples of training images whose production involved three
levels of hierarchy. First, an orientation (i.e., horizontal or vertical)
is randomly chosen with fair odds. Second, each bar of the chosen
orientation is randomly instantiated with probability 0.25. Third,
additive noise is introduced by randomly turning on with a
probability of 0.25 each pixel that was previously off.
4.1.1 Wake-sleep parameter estimation
Using the wake-sleep algorithm, I trained a stochastic binary Helmholtz machine that has 4 top-layer ("meta-cause") variables, 36 middle-layer variables, and 256 bottom-layer image variables. Each conditional distribution is modelled using logistic regression. Learning is performed through a series of iterations, where each iteration consists of one bottom-up wake phase sweep used to adjust the generative network parameters and one top-down sleep phase sweep used to adjust the recognition network parameters. Every 5000 iterations, the recognition network is used to obtain an estimate (with error bars) of the lower bound on the data log-likelihood under the generative network. To do this, 1000 recognition sweeps are performed without learning. During each recognition sweep, binary values for the hidden variables are obtained for the given training image. The log-probability of the values of all the variables under the generative network minus the log-probability of the hidden variable values under the recognition network gives an unbiased estimate of the log-likelihood bound (3.40). In this way, I obtain 1000 i.i.d. noisy unbiased estimates of


Page 91
Figure 4.2
Variation of the lower log-likelihood bound with the
number of wake-sleep learning iterations
for the stochastic Helmholtz machine.
the log-likelihood bound. The average of these values gives a less noisy unbiased estimate. Also, the variance of this estimate is estimated by dividing the sample variance by 999.
I am interested in solutions where the generative network can construct the image by adding features, but cannot remove previously instantiated features. If the network parameters are in no way constrained to favor this type of solution, perceptually unattractive solutions are found (see Section 4.1.3). So, I constrain the parameters of the logistic regression model that connects the middle layer to the bottom layer to be positive by setting to zero any negative weights every 20th learning iteration. In order to encourage a solution where each image can be succinctly described by the minimum possible number of causes in the middle layer, I initialize the middle-layer generative biases to -4.0 which favors most middle-layer variables being inactive (value of 0) on average. All other parameters are initialized to 0.0. For the first 100,000 iterations, a learning rate of 0.1 is used for the generative parameters of the model feeding into the bottom layer and for the recognition parameters of the model feeding into the middle layer; the remaining learning rates are set to 0.001. After this, learning is accelerated by setting all learning rates to 0.01.
Figure 4.2 shows the learning curve for the first 300,000 iterations of a simulation consisting of a total of 1,000,000 iterations. Aside from several minor fluctuations, the wake-sleep algorithm maximizes the log-likelihood bound in this case. Eventually, the bound converges to the optimum value (-170 bits) shown by the solid line. This value is computed by estimating the average log-


Page 92
likelihood of the data under the method that was used to produce the data, i.e., the negative entropy of the training data.
By examining the generative parameters after learning, we see that the wake-sleep algorithm has extracted the correct 3-level hierarchical structure. Figure 4.3 shows the parameters for the generative logistic regression models feeding into and out of the middle layer in the network. A black blob indicates a negative parameter and a white blob indicates a positive parameter; the area of each blob is proportional to the magnitude of the parameter (the largest value shown is 7.77 and the smallest value shown is -7.21). There are 36 blocks arranged in a 6 ´ 6 grid and each block corresponds to a middle-layer variable. The 4 blobs at the upper-left of a block show the parameters that connect each of the top-layer variables to the corresponding middle-layer variable. The single blob at the upper-right of a block shows the bias for the corresponding middle-layer variable. The 16 ´ 16 matrix that forms the bulk of a particular block shows the parameters that connect the corresponding middle-layer variable to the bottom-layer image. These matrices clearly indicate that 32 of the 36 middle-layer variables are used by the network as "feature variables" to represent the 32 possible bars. These feature variables are controlled mainly by the right-most top-layer "orientation" variable - the parameters connecting all the other top-layer variables to the feature variables are nearly zero. If the orientation variable is off, the probability of each feature variable is determined mainly by its bias. Vertical feature variables have significantly negative biases, causing them to remain off if the orientation neuron is off. Horizontal feature variables have only slightly negative biases, causing them to turn on roughly 25% of the time if the orientation variable is off. The parameters connecting the orientation variable to the vertical feature variables are significantly positive, so that when the orientation variable is on the total input to each vertical feature variable is slightly negative, causing the vertical feature variables to turn on roughly 25% of the time. The parameters connecting the orientation variable to the horizontal feature variables are significantly negative, so that when the orientation variable is on the total input to each horizontal feature variable is significantly negative, causing the horizontal feature variables to remain off. Since the parameters connecting the top-layer variable to the 4 middle-layer nonfeature variables are negative, and since the nonfeature variables have large negative biases, the nonfeature variables are usually inactive. Because the bottom-layer biases (not shown) are only slightly negative, a pixel that is not turned on by a feature variable still has a probability of 0.25 of being turned on. This accounts for the additive noise.


Page 93
Figure 4.3
Parameters for connections that feed into and out of the middle-layer
variables in the generative network. A black blob indicates a negative
parameter and a white blob indicates a positive parameter; the area of
each blob is proportional to the parameter's magnitude (the largest
value shown is 7.77 and the smallest value shown is -7.21).
4.1.2 Automatic clean-up of noisy images
Once learned, the recognition network can nonlinearly filter the noise from a given image, detect the underlying bars, and determine the orientation of these bars. To clean up each of the training images shown in Figure 4.1, I apply the learned recognition network to the image and obtain middle-layer variable


Page 94
Figure 4.4
Filtered versions of the training examples from Figure 4.1
extracted using the estimated recognition network.
values which reveal an estimate of which bars are on. The results of this procedure are shown in Figure 4.4 and clearly show that the recognition network is capable of filtering out the noise. Usually, the recognition network correctly identifies which bars were instantiated in the original image. Occasionally, a bar is not successfully detected. In two cases a bar is detected that has an orientation that is the opposite of the dominant orientation; however, usually the recognition network chooses a single orientation. Inspection of the original noisy training images for the two incorrect cases shows that aside from the single-orientation constraint, there is significant evidence that the mistakenly detected bars should be on. Further training reduces the chance of misdetection.
4.1.3 Wake-sleep estimation without positive parameter constraints
If all the parameters are initialized to 0.0, the parameters that connect the middle layer to the bottom layer are not constrained to be positive, and all the learning rates are set to 0.01, the estimated generative network does not properly extract the bar structure. Figure 4.5 shows the generative network parameters that connect the middle layer to the bottom layer, after 5,000,000 learning iterations. The black bars indicate that some middle-layer variables are capable


Page 95
Figure 4.5
Parameters for generative connections that feed out of the
middle-layer variables after estimation without special initialization
of the weights, without different learning rates between layers,
and without positive weight constraints in the generative network.
of "uninstantiating" bars that may be instantiated by other variables. Although it is imaginable that such a complex scheme is still capable of modelling the training images, the log-likelihood bound for this trained Helmholtz machine is -190 bits—significantly lower than the optimum value of -170 bits.
4.1.4 How hard is the bars problem?
Although this bar extraction problem may seem simple, it must be kept in mind that the network is not given a priori topology information—a fixed


Page 96
Figure 4.6
Training examples from Figure 4.1 after a fixed random
rearrangement of the pixels has been applied. These are
indicative of the difficulty of the bars problem in the absence
of a topological prior that favors local intensity coherence.
random rearrangement of the pixels in the training images would not change the learning performance of the network. Insofar as the network is concerned, the actual training examples look like those shown in Figure 4.6 which were produced by applying a fixed random rearrangement to the pixels in the images from Figure 4.1.
4.2 Simultaneous extraction of continuous and categorical structure
The Bayesian networks presented so far in this chapter have contained discrete variables. However, some hidden variables, such as translation or scaling in images of shapes, are best represented using continuous values. Work done on continuous-valued Bayesian networks has focused mainly on Gaussian random variables that are linked linearly such that the joint distribution over all variables is also Gaussian [Pearl 1988; Heckerman and Geiger 1995]. Lauritzen et al. [1990] have included discrete random variables within the linear Gaussian framework. They consider networks that are singly-connected, so


Page 97
Figure 4.7
(a) schematically shows the dependence of the proposed
variable on its parents. (b) to (e) illustrate four quite different modes
of behavior: (b) deterministic mode; (c) stochastic linear mode; (d)
stochastic nonlinear mode; and (e) stochastic binary mode (note the
different horizontal scale). For the sake of graphical clarity, the
density functions are normalized to have equal maxima and the
subscripts are left off the variables.
that probability propagation can be used. Most work on continuous-valued Bayesian networks requires that all the conditional distributions represented by the network can be easily derived using information elicited from experts. Hofmann and Tresp [1996] consider estimating continuous Bayesian networks that may be richly connected, but they assume that all variables are observed. As far as nonlinear continuous networks with latent variables are concerned, continuous-valued Boltzmann machines have been developed [Movellan and McClelland 1992], but these suffer from long simulation settling times and the requirement of a "negative phase" during learning. Tibshirani [1992], MacKay [1995] and Bishop et al. [1997] consider estimating mappings from a continuous latent variable space to a higher-dimensional input space, effectively using multiple-cause type networks of the form shown in Figure 3.4 on page 69. In this section, I consider a hierarchical Bayesian network with variables that can adapt to be continuous or categorical, as needed by the training data [Frey 1997a; Frey 1997b].
4.2.1 Continuous sigmoidal Bayesian networks
The proposed random variable is shown schematically in Figure 4.7a.1 The parents ai affect z i via a total input,
1 Geoffrey Hinton suggested this unit as a way to make factor analysis nonlinear.


Page 98
where we constrain q ij = 0.0 if z j Ï ai. The probability density over the presigmoid activity x i for variable z i is
where [image] is a parameter associated with variable z i (formally, part of q ). The value of z i (its postsigmoid activity) is obtained by passing the presigmoid activity through a sigmoidal cumulative Gaussian squashing function:
Networks of these variables can represent a variety of structures that are brought about by a range of significantly different modes of behavior available to each variable:
Deterministic mode: If the noise variance [image] is very small, the postsigmoid activity y i will be a practically deterministic sigmoidal function of the mean (see Figure 4.7b). This mode is useful for representing deterministic nonlinear mappings such as those found in deterministic multilayer perceptrons and mixture of expert models.
Stochastic linear mode: For a given mean, if the squashing function is approximately linear over the span of the added noise, the postsigmoid distribution will be approximately Gaussian with the mean and standard deviation linearly transformed (see Figure 4.7c). This mode is useful for representing Gaussian noise effects such as those found in mixture models, the outputs of mixture of expert models, and factor analysis models.
Stochastic nonlinear mode: If the variance of a variable in the stochastic linear mode is increased so that the squashing function is used in its nonlinear region, a variety of distributions are producible that range from skewed Gaussian to uniform to bimodal (see Figure 4.7d).
Stochastic binary mode: This is an extreme case of the stochastic nonlinear mode. If the variance [image] is very large, then nearly all of the postsigmoid probability mass will lie near the ends of the interval (0, 1) (see Figure 4.7e). E.g., for a standard deviation of 150, less than 1% of the mass lies in (0.1, 0.9). In this mode, the postsigmoid activity zi appears to be binary with probability of being "on"(i.e., z i > 0.5 or, equivalently, x i > 0):


Page 99
This sort of stochastic activation is found in binary sigmoidal belief networks [Jaakkola, Saul and Jordan 1996] and in the decision-making components of mixture of expert models and hierarchical mixture of expert models.
4.2.2 Inference using slice sampling
Assuming the variables are labeled in ancestral order, the joint distribution can be written
where N is the number of variables. and [image] are the presigmoid and postsigmoid conditional densities for variable zi. (Recall that the set of parents is represented by parameter constraints.) As usual, I define [image] to allow for a constant bias.
Even for small networks of these variables, probabilistic inference can be very difficult. Not only is the inference problem combinatorial, but it involves continuous hidden variables whose distribution when conditioned on visible variables may be multimodal with peaks that are broad in some dimensions but narrow in others. I use a Markov chain Monte Carlo procedure, which consists of sweeping a prespecified number of times through the set of hidden variables. A new value is obtained for each hidden variable using slice sampling [Neal 1997] (see Section 2.2.4), based on the distribution for the variable conditioned on all other variables. If an infinite number of slice samples are drawn for each hidden variable before passing on to the next hidden variable, this procedure is equivalent to Gibbs sampling (see Section 2.2.2). In fact, detailed balance still holds if only a fixed number of slice samples are drawn for each variable before passing on to the next variable. In most cases, drawing only one slice sample for each variable before continuing on will be most efficient.
If the parent-child influences cause there to be two very narrow peaks in the conditional distribution [image] for a hidden variable (corresponding to a variable in the stochastic binary mode), the slices will almost always consist of two very short line segments and it will be very difficult for the above procedure to switch from one mode to another. To fix this problem, slice sampling is performed in a new domain, y i = F ({x i - m i}/si). In this domain the parent-child distribution [image] is uniform on (0, 1), and


Page 100
thus [image]. It follows that I can use the following function for slice sampling:
where [image] . Since x i, y i and z i are all deterministically related, sampling from the distribution of y i will give appropriately distributed values for the other two.2
4.2.3 Parameter estimation using slice sampling
I use on-line stochastic gradient ascent to perform generalized EM (see Section 3.3.2). This consists of sweeping through the training set and for each training case following the gradient of [image] , while sampling hidden unit values as described above. The parameters are changed as follows:
where [image] is the learning rate.
I designed two experiments meant to elicit the four modes of operation described above. Both experiments were based on a simple network with one hidden layer h containing two variables and one visible layer v containing two variables. The training data was obtained by carefully selecting model parameters so as to induce various modes of operation and then generating 10,000 two-dimensional examples. Before training, the log-variances were initialized to 10.0, and the other parameters were initialized to uniformly random values between -0.1 and 0.1. Training consisted of 100 epochs using a learning rate of 0.001 and 20 sweeps of slice sampling to complete each training case. Each task required roughly five minutes on a 195 MHz MIPS R4400 processor.
The distribution of the training cases in visible unit space [image] for the first experiment is shown by the contours in Figure 4.8a. After training the network, I ran the inference algorithm for each of ten representative training cases. The mean postsigmoid activities of the two hidden units are shown beside the cases in Figure 4.8a; clearly, the network has identified four classes
2 Both F( ) and F-1( ) do not have closed-form expressions, so I use the C-library erf( ) function to implement F( ) and table lookup with quadratic interpolation to implement F-1( ).


Page 101
Figure 4.8
For each experiment (a) and (b), contours show the distribution of the
2-dimensional training cases. The inferred mean postsigmoid activity of the
two hidden units after learning are shown in braces for several training
cases, marked by ´.
that it labels (0, 0), (0, 1), (1, 0) and (1, 1). Based on a 30 ´ 30 histogram, the relative entropy between the training set and data generated from the trained network is 0.02 bits. Figure 4.8b shows a similar picture for the second experiment, using different training data. In this case, the network has identified two categories that it labels using the first postsigmoid activity. The second postsigmoid activity indicates how far along the respective ''ridge'' the data point lies. The relative entropy in this case is 0.04 bits.
The above experiments illustrate how the same network can be used to model two quite different types of data. In contrast, a Gaussian mixture model would require many more components for the second task as compared to the first. Although the methods due to Tibshirani [1992] and Bishop et al. [1997] would nicely model each submanifold in the second task, they would not properly distinguish between categories of data in either task. MacKay's method may be capable of extracting both the submanifolds and the categories, but I am not aware of any results on such a dual problem.
It is not difficult to conceive of models for which naive Markov chain Monte Carlo procedures will become fruitlessly slow. In particular, if two variables are highly correlated, the procedure of sampling one variable at a time will converge extremely slowly. Also, the Markov chain method may be prohibitive for larger networks. One approach to avoiding these problems is to use the Helmholtz machine or variational methods. In the next section, I present a general variational method for inference and learning in networks of nonlinear stochastic units.


Page 102
4.3 Nonlinear Gaussian Bayesian networks (NLGBNs)
In this section, I present a unifying framework for stochastic neural networks with nonlinear hidden variables. Nonlinear units are obtained by passing the outputs of linear Gaussian units through various nonlinearities. The continuous sigmoidal networks described above are a special case of these general nonlinear networks. A general variational method can be used to maximize a lower bound on the likelihood of a training set in these networks and the bound can be used for pattern classification (see Frey and Hinton [ 1998] for details).
4.3.1 The model
Linear generative models, binary sigmoidal Bayesian networks, continuous sigmoidal belief networks, and piecewise linear Bayesian networks (such as the "rectified Gaussian Bayesian networks" in Hinton and Ghahramani [ 1997]) can all be viewed as networks of Gaussian units that apply various nonlinearities to their Gaussian states. For example, a rectified unit adds Gaussian noise to the top-down input from its parents, and then outputs this value if it is positive and outputs 0 if the value is negative.
The probability density function over the pre-nonlinearity variables x = (x1, . . . x N ) in such a nonlinear Gaussian belief network (NLGBN) is
where A i is the set of indices for the parents of unit i and f(×) is the normal density function:
is the variance of the Gaussian noise for unit i and ¦j(×) is the nonlinear function for unit j. For example, some units may use a step function (making them binary sigmoidal units with a cumulative Gaussian activation function), whereas other units may use the rectification function (making them real-valued units that encourage sparse representations). Define ¦0(x0) = 1 so that w i 0 represents a constant bias for unit i in (4.8).
In the following section, I present a generalization of the variational method developed by Jaakkola et al. [1996] for networks of binary units and show that it can be successfully applied to performing approximate inference and learning in networks with other nonlinearities. The variational method can still be applied when different types of nonlinearity are used in the same network.


Page 103
4.3.2 Variational inference and learning
One way to infer the distribution over the hidden variables in an intractable Bayesian network is to postulate a simple parametric variational distribution q(·) over the hidden variables. (The variational distribution q(·) is separate from the generative distribution p(·).) A numerical optimization method is then used to adjust the variational parameters to bring q(·) as "close" to the true posterior as possible. See Section 2.3 for a description of this inference technique. In fact, the same variational distribution can be used to perform learning via the generalized EM method described in Section 3.3.2.
Let V be the set of indices of the sensory input variables for the current input pattern, and let H be the set of indices of the hidden variables for the current input pattern, so that V È H = {1, . . . , N}. The variational bound from (3.22) is
where (·) indicates an expectation over the hidden variables with respect to q(×).
A simple parametric variational distribution is a product of Gaussian distributions:
where m i and si, i Î H are the variational parameters. By adjusting these parameters, we can obtain an axis-aligned Gaussian approximation to the true posterior distribution over the hidden variables.
In the case of NLGBNs and the variational distribution consisting of a product of Gaussians, (4.10) simplifies to
In order to make this formula concise, I have introduced dummy variational parameters for the observed variables: if x i is observed to have the value [image] we fix m i = [image] and s i= 0. In (4.12), M j(m,s)is the mean output of unit j


Page 104
when the input is Gaussian noise with mean m and variance s2:
V j(m, s) is the variance at the output of unit j when the input is Gaussian noise with mean m and variance s2:
I assume that these can be easily computed, closely approximated, or in the case of V j(×,×), bounded from above (the latter will give a new lower bound on B q||p. See Frey and Hinton [1998] for these functions in the case of linear units, binary units, rectified units, and sigmoidal units.
For unit i in (4.12), the term in curly braces measures the mean squared error under q(×) between m i and the input to unit i as given by its parents: . It is down-weighted by the model noise variance [image], since the importance of the prediction error decreases with noise variance.
Variational inference consists of first fixing m i = and s i Î V in (4.12) and then maximizing B q||p with respect to m i and log. , i Î H. (The optimization for the variances is performed in the log-domain, since log [image] is allowed to go negative.) I use the conjugate gradient method to perform this optimization, although other techniques could be used (e.g., steepest descent or possibly a covariant method [Amari 1985]). The derivatives of B q||p with respect to m i and log [image], i Î H are given in Frey and Hinton [1998]. After optimization, the means and variances of the variational distribution represent the inference statistics.
For learning, the log-likelihood of the entire training set is bounded from below by a sum of bounds for the individual training patterns. The variational expectation maximization algorithm consists of the following:
· E-Step: Perform variational inference by maximizing the total bound with respect to the multiple sets of variational parameters corresponding to the different input patterns.
· M-Step: Maximize the total bound with respect to the model parameters (w..'s and s.'s).
Notice that by maintaining sufficient statistics while scanning through the training set in the E-Step, it is not necessary to store the multiple sets of variational parameters.


Page 105
It turns out that the M-Step can be performed very efficiently. Since the values of the model variances do not affect the values of the weights that maximize B q||p in (4.12), we can first maximize B q||p with respect to the weights. As pointed out by Jaakkola et al. [1996] for their binary Gaussian belief networks, B q||p is quadratic in the weights, so we can use singular value decomposition to solve for the weights exactly. Next, the optimal model variances are computed directly.
Software that implements this variational learning algorithm for linear units, binary units, rectified units, and continuous sigmoidal units is available by following pointers at http://mitpress.mit.edu/book-home.tcl?isbn=026206202X.
4.3.3 Results on the continuous stereo disparity problem
In this section, I consider an unsupervised visual feature extraction task, and compare the representations learned by the variational method applied to two types of NLGBN and the representations learned by Gibbs sampling applied to a piece-wise linear NLGBN. If the hidden units all use a piece-wise linear activation function, then Gibbs sampling can be efficiently used for learning, as described in Hinton and Ghahramani [1997]. Zoubin Ghahramani provided the MATLAB software for Gibbs sampling. One of the NLGBNs used for variational learning contains only binary hidden units of the type described by Jaakkola, Saul and Jordan [1996]. So, in this section we see how the variational technique compares to another learning method for continuous hidden units, as well as how the generalization of the variational method from binary to continuous units compares to variational learning in binary networks.
An interesting vision problem where the latent causes are nonlinearly related to the input is the estimation of depth from a stereo pair of sensor images. In the simplified version of this problem presented in Becker and Hinton [ 1992], the goal is to learn that the visual input consists of randomly placed dots on a surface at one of two depths. These experiments involved 1-dimensional "eyes" that had 12 real-valued sensors each. The eyes viewed a 1-dimensional surface on which 4 dots with intensity drawn uniformly from [0, 5) were randomly positioned. The surface was placed at one of two different depths with equal probability, so that the rays entering the two eyes were relatively shifted leftward or rightward. The 12 sensors in each eye used Gaussian activation filters and added unit-variance Gaussian noise. Twelve examples of the sensory images obtained in this manner are shown in Figure 4.9a, where the 1-dimensional sensory image for the right eye is shown above the image for the left eye.
The stereo disparity problem is much more difficult than the bars problem discussed in Section 4.1, since there is more overlap between the underlying


Page 106
Figure 4.9
Learning in NLGBNs using the variational method and Gibbs sampling, for
noisy stereo disparity patterns. (a) shows some training examples. (b) shows
the middle-layer to visible-layer weights learned by the variational method
with rectified units. (c) shows the weights learned by the variational method
with binary units. (d) shows the weights learned by Gibbs sampling with
rectified units.
features. To see this, imagine a "multi-eyed" disparity problem in which there are as many eyes as there are 1-dimensional sensors. We expect the depth inference to be easier in this case, since there is more evidence for each of the two possible directions of shift. Imagine stacking the sensory images on top of each other, so that each resulting square image will contain blurred diagonal bars that are oriented either up and to the right or up and to the left. Extracting disparity from this data is roughly equivalent to extracting bar orientation in the data from the previous section.


Page 107
100 iterations of variational EM were used to train a three-layer NLGBN with 1 binary top-layer unit, 20 rectified middle-layer units, and 24 linear visible bottom-layer units. The resulting weights projecting from each of the 20 middle-layer units to the two sets of 12 sensors are shown in Figure 4.9b. The algorithm has extracted features that are spatially local and represent each of the two possible depths. Figure 4.9c shows the weights learned by variational EM in a network where all of the hidden units are binary. The lower bound on the likelihood for the trained binary-rectified-linear NLGBN is -43.7 nats, whereas the bound for the trained NLGBN with binary hidden units is -46.1 nats, indicating that the continuous-valued hidden variables give superior likelihood bounds. The Gibbs sampling method with rectified hidden units produced the weight patterns shown in Figure 4.9d. (It is not straightforward to estimate the data likelihood using Gibbs sampling.)
4.3.4 Pattern classification using the variational bound
In contrast to the Monte Carlo method, the variational method provides an easily computed approximation to log p({x i}i Î v) in the form of the variational bound. If one NLGBN is trained on each class of real-valued data, a test pattern can be classified by selecting the class whose corresponding NLGBN gives the highest bound. This method is used in [Frey and Hinton 1998] to show that multi-layer NLGBNs with rectified units are competitive at classifying grey-level images of handwritten digits.


Page 109
5
Data Compression
The goal of data compression is to exploit the redundancy in input patterns so as to represent individual patterns concisely on average. In this book, I focus on lossless data compression, in which the original pattern can be completely recovered from the concise representation. A source code maps each input pattern v to a codeword, such that for each valid codeword there is a unique input pattern. I will consider sources where the patterns are i.i.d. (independent and identically drawn) from a distribution Pr(v).
Shannon's noiseless source coding theorem [Shannon 1948] states that the average codeword length cannot be less than the entropy of the source:
where [image] is the length in bits of the codeword for input pattern v, and [image] is the entropy of the source:
Traditional approaches to data compression have focused on producing source codes whose codeword lengths are nearly optimal, where the optimal length of the codeword for v is log2 Pr(v).
Arithmetic coding [Rissanen and Langdon 1976; Witten, Neal and Cleary 1987] is a practical algorithm for producing near-optimal codewords when the source distribution Pr(v) is known. If v is binary-valued, Pr(v) can be easily estimated and arithmetic coding can be used to produce near-optimal "fractional bit" codewords. If v is high-dimensional and the distribution is quite complex (e.g., images of faces), it may be desirable to estimate a more sophisticated flexible probability model P(v). Unfortunately, even if such a model can be estimated so that P(v) » Pr(v), there may not be a practical way to encode v using the model. For example, an arithmetic encoder requires a table of the probabilities for all possible inputs. For a 16 ´ 16 binary image, this table would have 2256 entries! Not only do we need a model that provides a probability P(v), but we also need a model that somehow decomposes P(v) in a way that allows the encoder to encode the variables one at a time (or in small groups).
Graphical models provide a structured description of P(v), and so they seem like a good place to look for the source models described above. How-


Page 110
ever, it turns out that Markov random fields do not provide the right type of structure. For example, the Boltzmann machine [Hinton and Sejnowski 1986] (a Markov random field that learns) is poorly suited to data compression, because it does not decompose P(v) in a way that is suitable for efficient piecewise compression. (A method such as Markov chain Monte Carlo must be used to estimate the partition function, which normalizes the probabilities.) On the other hand, Bayesian networks do provide an ideal structure for data compression.
In Section 5.1, I show how Bayesian network source models that do not have latent variables can be used very efficiently to compress data. Then, in Section 5.2, I go on to discuss source models that have many latent variables. Values can be chosen for the latent variables and the entire configuration can be encoded. In this way, a "multi-valued source code" with many codewords for each input pattern is obtained. In many cases, these codewords cannot be mixed together in a tractable way. To remedy this problem, I show how extra information can ride "piggyback" on the choice of codeword and derive the communication rate for this "bits-back" procedure. In Section 5.3, we see that a broad class of approximations to maximum likelihood parameter estimation actually minimizes this communication rate. In Section 5.4, I outline the ''bits-back coding" algorithm, which is a practical implementation of the "bits-back'' idea. It turns out that by using an arithmetic decoder in the bits-back encoder and an arithmetic encoder in the bits-back decoder, we can achieve a practical communication rate that is nearly optimal. Finally, in Section 5.5, I present compression results for Helmholtz machine source models that are adapted using the wake-sleep algorithm.
5.1 Fast compression with Bayesian networks
Suppose we have at hand a Bayesian network for the binary variables in v, such that P(v) » Pr(v). As discussed in Section 1.2.3, the joint distribution for a Bayesian network can be written
where N is the number of variables in v, and ak are the parents of v k. This decomposition is very well-suited to arithmetic coding.
In order to encode v, we pick an ancestral order for the variables. I will assume without loss of generality that v1, . . . v N is such an ordering. Compression begins with v1, whose observed value is fed into the arithmetic en-


Page 111
coder, along with its distribution{P(v1 = 0), P(v1 = 1)}, which is part of the network specification. Next, we compute {P(v2 = 0|a2), P(v2 = 1|a2)} using the conditional probability given in the network specification as well as the values of v2's parents (i.e., either {v1} or Æ) which have already been encoded. We feed the observed value of v 2 into the arithmetic encoder, along with its distribution. Encoding continues in this fashion until all network variables have been encoded.
For this procedure to work as described, the Bayesian network must be fully visible. That is, all network variables are part of the input pattern. Suppose there are some latent variables h that are not part of the input. Then, the Bayesian network models P(v, h). These variables may be important for representing higher-order structure in the input v, as discussed extensively in Chapter 3. Now, the decomposition in (5.3) cannot be used.
If there aren't many latent variables, we can use a procedure that is similar to the one described above. We pick an ancestral order and proceed as described above, encoding only the observed variables. Whenever we encounter an observed variable that is not dependency-separated from an unobserved variable by the variables that have been encoded so far, the unobserved variable must be integrated out, by summing over its values. The complexity of this encoding procedure is usually exponential in the number of unobserved variables. Sometimes, the graphical structure of the network permits this procedure to be done in a very efficient way. For example, the latent variables in a hidden Markov model [Rabiner 1989] with a fixed state space size can be integrated out in a way so that the encoding complexity is linear in the number of latent variables (number of time steps).
5.2 Communicating extra information through the codeword choice
In general, when the latent variables in a Bayesian network cannot summed away to compute P(v) in a tractable way, we are left with the option of picking values for them. Then, the entire set of variables {v, h} can be encoded using the procedure described above, as if all values were observed. This can be viewed as a multi-valued source code, where there are many codewords for each input v. The codeword depends on which values are chosen for the latent variables. Often, as with the hidden Markov model, these codewords can be mixed in an efficient way. However, there is an interesting class of multi-valued source codes (e.g., Bayesian networks with latent variables) for which the multiple codewords cannot be mixed in a tractable manner.


Page 112
Figure 5.1
A scheme in which auxiliary data is communicated along with the symbol in
order to achieve optimal compression when the source code is multi-valued.
The approach I take to solve this problem [Frey and Hinton 1996; Frey and Hinton 1997] is motivated by the "bits-back" argument of Hinton and Zemel [1994], which they used to develop a Lyapunov function for machine learning. It turns out that Wallace [1990] devised a similar argument to construct minimum-length integer-length messages for use in minimum-message-length inference. By selecting codewords through the use of extra auxiliary data, the auxiliary data can ride "piggyback" on the codewords for the symbols that we are encoding. Compared to the optimal single-valued source code obtained by mixing together the codewords for an input pattern, the bits communicated in the auxiliary data will make up for the lengths of the suboptimal codewords that are sent. In particular, the communication rate will be less than the rate that would be obtained by always picking the shortest codeword. A block diagram for this communication scheme is shown in Figure 5. 1. A simple example will help illustrate this procedure.
5.2.1 Example: A simple mixture model
Consider a source that outputs real numbers that are distributed according to a mixture of two Gaussians. These numbers are rounded to some precision to form a set of symbols. The component distributions and the output distribution are shown in Figure 5.2a, where the rounding effect is left out for the sake of graphical simplicity.
The most natural source code to use in this case is one that requires one bit to specify from which Gaussian a given symbol was produced plus however many bits are needed to code the symbol using that Gaussian. However, the identity of the Gaussian that produced a given symbol is often ambiguous. In particular, a number near v0 could well have come from either Gaussian. In these cases the source model maps each symbol to two codewords—one for each Gaussian—producing a multi-valued source code. If we were to always pick the shorter of the two codewords, we would effectively be assuming the symbols were distributed as in Figure 5.2b. However, this distribution is obviously incorrect—it is not even normalized—and will lead to suboptimal compression.


Page 113
Figure 5.2
The most natural source model may produce multiple
codewords for a given symbol. (a) shows a source with a single
binary hidden variable which identifies from which Gaussian,
G1 or G2, the symbol value v is sampled. Values of v near v0 are
likely to have come from either Gaussian. (b) shows the resulting
coding density effectively used if we were to always pick the
shorter codeword. This density wastes coding space because it
is wrongly shaped and has an area significantly less than unity.
The obvious way around this problem is to use a single-valued code that is based on a summation of the mixture component probabilities. That is, we assign a new codeword to each symbol based on its total probability mass, obtained by summing the contributions from each Gaussian. Although this procedure is obviously computationally feasible for this example, there are more complicated models where it is not (see Section 5.5). In fact, the same rate can be achieved by using the original multi-valued source code and communicating extra information through the choice of codeword. This may seem surprising, since for a given symbol both codewords in the multi-valued source code are longer than the codeword in the single-valued source code.
Consider a sender that wishes to encode a rounded value v' that requires 2 bits if encoded using G1 and requires 3 bits if encoded using G2 (i.e., v' is twice as likely under G1 as it is under G2). Including the single bit required to specify which Gaussian is being used, an optimal source code (where the Gaussian identity is explicit) will thus have codewords with lengths [image] = 3 bits and [image] = 4 bits. If the sender always picks the shorter codeword, the average codeword length is 3 bits.
Suppose instead that whenever the sender must communicate the particular symbol v', the sender chooses between the two codewords with equal probability. (In general, the ratio of choices will depend on v'.) It appears the average codeword length in this case is ( [image] + [image]) = 3.5 bits, which is higher than that obtained by always choosing the shortest codeword. However, this cost is effectively lowered because the receiver can recover information from the


Page 114
choice of codeword in the following manner. Say the sender has high-entropy auxiliary data available in the form of a queued bit stream with 0 and 1 having equal frequency. When encoding v', the sender uses the next bit in the auxiliary data queue to choose between G1 and G2. The sender then produces a codeword that will have an average length of 3.5 bits (it is important to note that this codeword specifies which of G1 and G2 is being used).
When decoding, the receiver reads off the bit that says which Gaussian was used and then determines the rounded value v' from the codeword. Given the decoded value, the receiver can run the same encoding algorithm that the sender used, and determine that a choice of equal probability was made between G1 and G2. Since the receiver also knows which Gaussian was selected, the receiver can recover the queued auxiliary data bit that was used to make the choice. In this way, on average 1 bit of the auxiliary data is communicated at no extra cost. I refer to these recovered bits as bits-back.
If the auxiliary data is useful, the average effective codeword length is reduced by 1 bit due to the savings, giving an effective average length of 2.5 bits—less than the 3 bits required by the shortest codeword. I refer to this method of source coding as bits-back coding. It is important to note that the ratio of choices between G1 and G2 depends on the symbol being encoded. For example, if the rounded value is far to the right of v0 in Figure 5.2a, then picking the codewords equally often would be very inefficient, since the codeword under G1 would be extremely long, making the benefit of the single recovered bit negligible. In this case the sender should pick G1 much less often and as a result the sender will read off only "part" of a bit from the auxiliary data queue to determine which codeword to use. As we will see below, choosing between the two codewords with equal probability is not optimal in the above example.
5.2.2 The optimal bits-back coding rate
The rate for bits-back coding can be determined by defining a distribution that is used to select codewords for a given input symbol (pattern), v:
where h is a binary vector representing the index of the selected codeword for input v. (It is represented as a vector, since it too must be encoded.) Letting [image] (v,h) be the length of the hth codeword1 for a specific pattern v, the expected length of the two-part codeword for v is
1 The codewords may have fractional lengths produced, say, by arithmetic coding


Page 115
The expected bits-back for v is the information content (entropy) of the distribution used to select codewords:
The difference between (5.5) and (5.6) gives the communication cost for v:
The overall rate [image] for bits-back coding is given by averaging this cost over the source distribution, Pr (v):
It is easily proven from (5.8) that for each v the codeword selection distribution which minimizes the bits-back coding rate is the Boltzmann distribution:
I denote by * those quantities determined from the Boltzmann distribution. Note that this distribution depends on the input symbol, v. The optimal rate for a given multi-valued source code is achieved if for each input symbol a codeword is selected using the corresponding Boltzmann distribution. By substituting (5.9) into (5.8), we find that the optimal bits-back coding rate is
This rate is the same as the rate for a single-valued source code that has codeword lengths which properly reflect the total codeword space associated with each symbol in the multi-valued source code.
In the mixture of Gaussians example, where for symbol v' we had [image] = 3 bits and [image] = 4 bits,


Page 116
This is the minimum [image] (v') for the given example. A slightly higher than optimal [image] (v') of 2.5 bits was obtained above using [image] (G1|v' ) = [image] (G2|v' ) = 0.5.
5.2.3 Suboptimal bits-back coding
For complex source models, the summation in the denominator of (5.9) is usually intractable; in these cases, it is not possible to obtain the exact Boltzmann distribution. When it is possible to obtain the exact Boltzmann distribution, the denominator in (5.9) can often be directly used to create a single-valued source code. The advantage of bits-back coding is that when the multi-valued source code is unmixable, an approximation to the Boltzmann distribution can be used. There are a variety of practical algorithms for obtaining such an approximation, including Markov chain Monte Carlo methods [Geman and Geman 1984; Hinton and Sejnowski 1986; Ripley 1987; Potamianos and Goutsias 1993], variational methods [Chandler 1987; Peterson and Anderson 1987; Zhang 1993; Saul, Jaakkola and Jordan 1996], and inverse model methods [Hinton et al. 1995; Dayan et al. 1995] (see Section 5.5). The rate for an arbitrary codeword selection distribution Q(h|v) can be compared to the optimal rate given by the Boltzmann distribution:
This is the information divergence (a.k.a. relative entropy) between the codeword selection distribution and the Boltzmann distribution, averaged over the source distribution. It is always non-negative and yields the increase in coding rate caused by the approximation to the Boltzmann codeword selection distribution.
A suboptimal selection distribution of particular interest is Qshort(h|v), which always picks the shortest codeword, hshort(v). (This is analogous to the two-


Page 117
part codes discussed by Rissanen [1989].) In this case, the rate increase compared to the optimal rate is
Bits-back coding makes gains over shortest codeword selection by approximately taking into account the entire codeword space associated with an input symbol, as opposed to just the codeword space associated with the shortest codeword. If several of the shortest codewords have roughly equal lengths, or if there are a large number of codewords with lengths somewhat larger than the shortest, then Q*(hshort(v)|v) is significantly less than unity indicating that picking the shortest codeword is far from optimal. Relative to Rissanen's work, bits-back coding provides a tractable way to approximate the stochastic complexity [Rissanen 1989] and furthermore communicate at this rate.
5.3 Relationship to maximum likelihood estimation
The whole idea of a multi-valued source code may seem absurd. Why waste codeword space by associating multiple codewords with each symbol? An answer to this question must be closely related to the structure of the source model. In addition to the input pattern being encoded, it is often useful and natural to consider extra latent variables whose purpose is to capture high-order structure. For example, when modelling grey-scale images, it may help to create a latent variable that measures overall image contrast. The codeword for a particular image will include a binary representation of this contrast value. However, there may be several quite different contrast values that are equally plausible, leading to several different codewords.
A generative model of the type described above typically provides a parameterized distribution P(h|q H) that can be used for encoding the set of latent variables h, as well as a distribution P(v|h,q V) to be used for encoding the input symbol v for a given setting of the latent variables. Such a codeword will have an optimal length (e.g., obtained using arithmetic coding) given by


Page 118
Note that the generative structure implies that P(v|h,q V) is easy to compute. (Rissanen [1989] refers to this type of code as a two-part code.)
The set of parameters q = {q H, q V} must be fixed by hand, estimated using a stored data set, or adapted on-line. Estimating these parameters is a difficult task when there are latent variables. The popular technique of maximum likelihood estimation minimizes the following cost:
Combining (5.14) and (5.15), we find that maximum likelihood estimation minimizes
which is equal to the optimal bits-back coding rate (5.10). In contrast, maximum likelihood estimation does not minimize the rate for an encoder that always picks the shortest codeword.
Often, maximum likelihood estimation is not tractable when the generative model is overly complex. In these cases, it is possible to use various approximations to maximum likelihood estimation. A common approach [Peterson and Anderson 1987; Neal and Hinton 1993; Zhang 1993; Hinton et al. 1995; Dayan et al. 1995; Saul, Jaakkola and Jordan 1996] is to minimize an upper bound on C, thus guaranteeing that the cost is lower than a certain value. (This is described in Section 3.3, where it is called maximum likelihood-bound estimation.) The logarithmic term in (5.16) is first bounded by introducing an extra distribution Q(h|v) and using Jensen's inequality:
Inserting this bound into (5.16), we get an upper bound on C:


Page 119
Figure 5.3
Feeding a random number into an arithmetic decoder with appropriate
probabilities (shown by the partition heights within a column) selects code-
words (shaded partitions), while at the same time conserving information.
which is equal to the suboptimal bits-back coding rate (5.8). So, these methods—including the algorithms presented in Section 3.4—minimize the suboptimal bits-back coding rate. As with exact maximum likelihood estimation, these methods do not minimize the rate for an encoder that always picks the shortest codeword.
5.4 The "bits-back" coding algorithm
To implement the communication scheme shown in Figure 5.1, we need a general method of recovering the auxiliary data bits from the codeword choices. In the mixture of Gaussians example, we considered a specific input symbol for which there were two codewords. These codewords were selected equally often so that a single bit could be used for bits-back. If the codeword selection distribution is dyadic2, Huffman decoding [Huffman 1952] can be used to pick codewords. In this section, I consider the case of an arbitrary codeword selection distribution. Software that implements this bits-back coding algorithm can be found by following the pointers given at http://mitpress.mit.edu/book-home.tc1?isbn=026206202X.
In the case of an arbitrary codeword selection distribution, it is not obvious how random codeword choices can be made without losing auxiliary data information. To address this problem, consider the operation of an arithmetic
2 i.e., each probability is an integral power of 2.


Page 120
Figure 5.4
The block diagrams for (a) the bits-back encoder and (b) the bits-back decoder.
decoder [Rissanen and Langdon 1976; Witten, Neal and Cleary 1987]. It receives a finite-precision number on [0,1) and extracts from it a series of decisions according to a table of probabilities. If a collection of uniformly distributed finite-precision numbers on [0, 1) is decoded in parallel, we will obtain a collection of decisions whose distribution exactly matches the table of probabilities. Figure 5.3 shows how an arithmetic decoder can be used to conserve the information in the auxiliary data when making random codeword choices. The probabilities associated with the decisions form the table of the arithmetic decoder, while the auxiliary data defines a random number to be decoded. Each column of the figure corresponds to a single codeword choice and is partitioned into several possible outcomes with the height of each partition proportional to the probability of the corresponding outcome. It is easy to see that if the random number defined by the auxiliary data is uniform, codeword choices will be made according to the codeword selection distribution (as shown for a particular case by the shaded partitions). It is also easy to see


Page 121
Figure 5.5
The need for extra auxiliary data is eliminated by
feeding the codeword bits back into the bits-back
encoder as auxiliary data.
that by applying an arithmetic encoder to the sequence of decisions, we can regenerate the random number.
Figures 5.4a and 5.4b show block diagrams for the bits-back encoder and decoder respectively. When the bits-back encoder acquires a symbol v, it uses the codeword selection distribution Q(h|v) with an arithmetic decoder to choose codeword h, while consuming some auxiliary data bits. The multi-valued source code is then used to produce a codeword of length [image] (v| h). (For the experiments described in the next section, the multi-valued source code is implemented using an arithmetic encoder in conjunction with a source model distribution P(h|q )P(v| h,q ) as described in Section 5.3). When the bits-back decoder receives the codeword, it first decodes v and h using the multi-valued source code. It then uses the codeword selection distribution Q(h|v) with an arithmetic encoder to recover auxiliary data bits back from the codeword choice h.
Insofar as algorithm complexity goes, for an arbitrary codeword selection distribution, the codeword selection procedure described above requires an arithmetic encoder/decoder pair. If codewords are produced using arithmetic coding, the incremental cost of the codeword selection procedure is not overwhelming. In a hardware implementation, the codeword selection procedure can run in parallel with codeword production.
5.4.1 The bits-back coding algorithm with feedback
In practice, when encoding a block of symbols, extra auxiliary data is often not readily available. One solution to this problem is to use the binary form of a portion of the block of symbols for auxiliary data. However, so that the bits-back are efficiently utilized, this portion of symbols should first be source coded. Figure 5.5 shows a scheme for using the same multi-valued source code for doing just this, when the codewords have integer lengths. In order to encode a block of symbols, some initial primer bits (e.g., a few unencoded


Page 122
source symbols) are first placed in the queue. When the next symbol is bits-back encoded, some of the bits in the storage queue are used for auxiliary data. The resulting codeword is fed back into the storage queue so that it can (possibly) be used as auxiliary data later on. Once the entire block of symbols is encoded, the bits-back decoder proceeds to remove the codewords from the storage queue in reverse order. Since the decoder has no way of knowing a priori how long each codeword is, it is essential that the encoder reverse the bits within each codeword before feeding the codeword into the storage queue. The source symbols are decoded in reverse order compared to the order in which they were encoded. As decoding proceeds, the recovered bits-back are fed into the opposite end of the storage queue and will later be used as codeword bits or primer.
This method is inherently block-oriented, since each block must be decoded in the opposite order in which it was encoded. As a consequence, a block delay is introduced, which is often undesirable. Shorter block lengths will lead to extra overhead due to the primer and also due to framing information (such as a codeword used to indicate the end of the block). However, if the block delay is tolerable, this scheme nicely eliminates the need for extra auxiliary data.
When the multi-valued source code is implemented using arithmetic coding, the above feedback procedure cannot be used as defined. An arithmetic encoder produces a sequence of codeword bits and in general there is no way to break apart this sequence into pieces of integer length such that each piece corresponds to one symbol. This problem is easily solved by dividing the block of symbols into sub-blocks. The arithmetic encoder used to produce codeword bits is halted after each sub-block of symbols is processed. The resulting series of codeword bits is reversed and fed into the storage queue as described above. Practical arithmetic encoders usually waste only a few bits (2 in the implementation described in Witten, Neal and Cleary [19871) when encoding is terminated. The sub-block size should be chosen so as to minimize the effect of this wastage. For example, if the optimal bits-back coding rate is 1 bit/symbol, then choosing a sub-block size of 1000 symbols/sub-block will lead to a rate increase of only 0.2%. On the other hand, if the optimal bits-back coding rate is 1000 bits/symbol, arithmetic encoding can be terminated after each symbol (i.e., the sub-block size is 1 symbol/sub-block) and the rate will increase by only 0.2%.
5.4.2 Queue drought in feedback encoders
At first sight, it may appear that queue drought is a serious problem. This can occur if the arithmetic decoder in the bits-back encoder uses up all of the bits in the storage queue and still can't make a codeword choice. In fact, this is


Page 123
usually not a problem because practical arithmetic decoders/encoders [Witten, Neal and Cleary 1987] use a coding value with a restricted size (32 bits in my implementation). Consequently, in my implementation no more than 32 auxiliary data bits will ever be drawn from the storage queue when making a codeword choice. In degenerate cases where the codeword selection distribution places very little mass on one or more short codewords, it is possible for a queue drought to occur when a sequence of very short codewords are chosen that consistently draw a large number of bits each from the storage queue. However, even in such degenerate cases, the sequence of events that leads to a queue drought is highly atypical. I have found that in practice queue drought is not a problem, as long as a reasonable amount of primer (say 20 patterns) is used.
5.5 Experimental results
In this section, I present two sets of results for bits-back data compression. The source models are Helmholtz machines trained using the wake-sleep algorithm (see Section 3.4). The first data set consists of simple patterns of horizontal and vertical bars. The second data set consists of binary images of handwritten digits.
5.5.1 Bits-back coding with a multiple-cause model
In this section, I describe how bits-back coding can be applied to a binary Bayesian network source model, that has one layer of hidden binary variables. Then, I present compression results when the model is fit to images of horizontal and vertical bars using the wake-sleep algorithm described in Section 3.4.4. I compare the compression efficiency of the one-to-many bits-back source coding algorithm with the one-to-one source code obtained using approximate shortest codeword selection, and also with the UNIX gzip utility. The multi-valued source code has over 68 billion codewords for each input symbol, and there is no tractable way to mix them, as there is with a hidden Markov model. For a given symbol, most of these codewords are extremely long and therefore play a negligible role in the source code. However, it turns out that the rate for an algorithm that uses a tractable approximation to shortest codeword selection is significantly higher than the suboptimal bits-back coding rate. This indicates that multiple codewords should in some way be accounted for.
It turns out that there isn't an efficient way to convert the multi-valued source code for the sigmoidal Bayesian network into a single-valued source code that achieves a rate that is comparable to the bits-back coding rate. To perform such a conversion, we must compute most of the probability mass correspond-


Page 124
ing to the codewords for a given data vector. Because of the combinatorial way in which the latent variables h interact to produce P(v|h), the marginal probability mass P(v) cannot be computed in a tractable manner. v could be encoded bit by bit using Gibbs sampling to collect statistics. However, this procedure would require the computationally taxing simulation of a Markov chain for each element in v.
In order to use bits-back coding, we need a codeword selection distribution that is close to P(h|v,q ). The Helmholtz machine with the wake-sleep learning algorithm provides an estimate of the optimal codeword selection distribution. The learning algorithm jointly estimates the generative network P(v,h|q ) and a recognition network Q(h|v,f ) » P(h|v,q ). So, an input pattern can be encoded as follows. The sender first uses an ancestral order for the recognition network to compute the probability for the first latent variable (in the ancestral order). This probability and some auxiliary data are then fed into an arithmetic decoder which outputs a value for the first latent variable. Given the input pattern and the value for the first latent variable, the sender then computes the probability for the second latent variable, and so on. Once h has been chosen in this manner, the entire configuration for {v,h} is arithmetically encoded using the method described in Section 5.1.
The receiver decodes the entire configuration for {v,h} and then computes the probability for the first latent variable using the same ancestral order that was used by the sender. This probability and the value for the variable are fed into an arithmetic encoder. While this process is repeated for the remaining latent variables, the arithmetic encoder will output auxiliary data bits.
So that we can compare the performance of bits-back coding with the actual entropy rate of the source, I used a synthetic source to produce 6 ´ 6 binary images. The images are iid., and each image is produced by turning on each of the 12 possible horizontal and vertical bars with probability 0.2. (Both types of bars may appear in the same image.) The entropy rate of this source is 8.6 bits/image.
The multiple-cause network had a single hidden layer containing 36 binary variables, in addition to the visible layer containing 36 binary variables. In order to avoid the need for extra auxiliary data, bits-back coding with feedback was used (see Section 5.4.1). The images were grouped into sub-blocks of size 20 images/sub-block and a block size of 200 sub-blocks/block was used. Before each block was encoded, the first sub-block of binary images was used to prime the storage queue. After each block of images was communicated, both the encoding model and the decoding model were adapted using the wake-sleep algorithm with a gradient descent step size of 0.01. The parameters for both the generative network and the recognition network were initialized to 0.0 before any images were processed. I also approximated short-


Page 125
Figure 5.6
Experimental results for a Helmholtz machine with one hidden
layer of binary units applied to binary synthetic images.
est codeword selection by picking for each image v the configuration h that maximized Q(h|v,f ). (The quality of this approximation is discussed below.) Choosing the configuration h that maximizes Q(h|v,f ) can be done efficiently by considering one latent variable at a time. Figure 5.6 shows the number of codeword bits communicated as a function of the number of blocks encoded for both of these methods. The curves for the uncoded binary image data and the Shannon limit given by the entropy rate of the source are also given. The curve for the UNIX gzip utility with the "-best" option is shown for comparison. (Although the UNIX gzip utility is not really meant for image compression, I include it as a reference point for the reader.) It is evident that if we were to compare the Helmholtz machine with gzip, we would arrive at different conclusions depending on whether we used approximate shortest codeword selection or bits-back coding. The bits-back coding curve is clearly superior to the curve for approximate shortest codeword selection.
Table 5.1 gives a comparison of the rates obtained for the next block after 400 blocks of images were processed. The rate for approximate shortest codeword selection is significantly higher than the rate for bits-back coding, indicating that a significant practical savings can be made by using the new algorithm as opposed to shortest codeword selection. However, the communication rate for a logistic autoregressive network that was trained on-line (using


Page 126
Table 5.1
Rate comparisons of software-implemented codes for synthetic images.
Rate (bits/image)
Uncoded binary images
36.0
UNIX gzip utility (with ''-best'' option)
22.7
Approximate shortest codeword selection
(Helmholtz machine)
21.0
Bits-back coding (Helmholtz machine)
12.3
Logistic autoregressive network
10.7
Shannon limit
8.6

a learning rate of 0.01) is also given in Table 5.1, and is significantly lower than the rate for the Helmholtz machine. It appears that although bits-back coding opens the door to new multi-valued source codes, the ones studied in this section are not yet competitive with simpler compression methods.
How close does the approximation to shortest codeword selection come to actually picking the shortest codeword for each data vector? Since there are over 68 billion codewords for each image in the above example, we cannot make a direct comparison by actually searching for the shortest codeword. However, consider the same type of multiple-cause network, except with 9 hidden variables and 9 visible variables, applied to similar synthetic data, except with an image size of 3 × 3 and a bar probability of 0.1. This network is small enough that an exhaustive codeword search is possible. After processing 400 blocks of 1000 images each, I found that the approximation to shortest codeword selection gave a rate of 5.92 bits/image and exact shortest codeword selection gave a rate of 5.87 bits/image. These two rates are indistinguishable in the first decimal place. I expect that the results for the approximation used for the larger network are also close to the results that would have been obtained if an exhaustive search had been performed.
5.5.2 Compressing handwritten digits
The first step in efficient document compression is to segment the document into text components and nontext components. The text components may then be further segmented into isolated characters and then compressed using a model of the characters. Here, I consider modelling handwritten digits using a hierarchical Helmholtz machine and then compressing them using bits-back coding. Figure 5.7 shows 50 examples of the binary images that were used in this experiment. The binary Bayesian network source model had three hidden


Page 127
Figure 5.7
Examples of 8 ´ 8 binary images of handwritten digits.
Table 5.2
Rate comparisons for software-implemented source codes on the binary digit data.
Rate (bits/image)
Original binary file
64
Approximate shortest codeword selection
(Helmholtz machine)
60
gzip -best
39
Bits-back coding (Helmholtz machine)
33

layers of binary variables and one bottom layer of 64 visible variables. From top to bottom, the three layers of causes had 16 variables, 20 variables, and 24 variables, giving a total of 60 latent variables (260 codewords for each input pattern). Both the top-down and the bottom-up networks were fully-connected from layer to layer, but had no connections within each layer. The Helmholtz machine was fit to a training set of 100,000 images, using the wake-sleep algorithm.
A comparison of the average rates obtained on the training set using approximate shortest codeword selection and bits-back coding with the estimated binary Bayesian network, as well as the rates obtained by the UNIX gzip utility with the -best option, are given in Table 5.2. The rate for shortest codeword selection is again significantly higher than the rate for bits-back coding, indicating that a significant practical savings can be made by using bits-back coding.


Page 128
5.6 Integrating over model parameters using bits-back coding
As noted in Section 5.3, bits-back coding is closely related to statistical inference. In fact, the optimal bits-back coding rate is equivalent to Rissanen's stochastic complexity [Rissanen 1989] if we interpret the choice of codeword as a model parameter. Also, if the codewords are constructed by choosing a prior over codeword identities (P(h|q H) in (5.14)), bits-back coding effectively integrates over a discrete set of models.
Both of these relationships lead to an interesting application for bits-back coding. Suppose we are interested in encoding blocks of source symbols and that the source changes from block to block, but not within any single block. Given a model with a continuous parameter vector q , there is a single block codeword with length -log2 [image] for each block of source symbols, [image] . According to the principles of Bayesian analysis, we ought to encode [image] by integrating over the entire continuum of models, giving a codeword of length - [image] . In practice, this integral is usually impossible to compute and an approximation must be used. One approximation is to use the maximum a posteriori (MAP) model (i.e., [image]), for which the parameters [image] are communicated using some (hard to determine) precision.
In fact, bits-back coding can be used to communicate each block of symbols using the entire continuum of models, as long as a good approximation to the posterior distribution, [image] , is available. This distribution is used as the model selection distribution (in place of the codeword selection distribution) and the model parameters are communicated to an arbitrary precision. Whereas with the MAP approach, greater precision eventually leads to an increase in coding rate, with the bits-back coding approach, greater precision usually leads to a decrease in coding rate. Intuitively, this can be seen as an interaction of two processes. First, the extra codeword length caused by greater precision is partly recovered as bits-back. Second, greater precision usually leads to a more accurate approximation to the posterior distribution, and therefore shorter codewords on average. The latter process dominates except in the unusual case when the quantized version of [image] has a lower entropy relative to [image] than the unquantized version. I am currently exploring the use of bits-back coding for integrating over continuous parameter spaces.


Page 129
6
Channel Coding
Our increasingly wired world demands efficient methods for communicating discrete messages over physical channels that introduce errors. Examples of real-world channels include twisted-pair telephone wires, shielded cable-TV wire, fibre-optic cable, deep-space radio, terrestrial radio, and indoor radio. Each of these channels is subject to information-theoretic limitations, physical degradation, and governmental regulation. The prime information-theoretic limitation is Shannon's limit, which gives the maximum average number of information bits that can be communicated per second over a specific channel for a given set of transmitter constraints (e.g., transmission power). Examples of physical degradation include attenuation, thermal noise, self-interference (inter-symbol-interference), multiple-user interference, multiple-path radio reflections, and power limitations in practical circuits. Examples of governmental regulations include transmission power limits, bandwidth usage, and information packet sizes. Together, all of these restrictions and many more define the practical channel coding problem of how to communicate discrete messages reliably.
Despite the multi-faceted nature of the practical channel coding problem, the bottom line is nonetheless quite straightforward. (See MacKay [1998] for an excellent introduction to information theory and its connections with probabilistic inference.) In order to communicate, the transmitter sends a finite-duration real-valued signalling waveform. This waveform is determined by a binary information sequence, which we usually assume is uniformly distributed over all possible information sequences. The duration of this waveform may correspond to a relatively short block of information or an infinite-length limiting-case block of information. Once the transmitter has produced a signalling waveform, it is transformed stochastically by the channel and a received waveform or channel output waveform is obtained at the output of the channel. The receiver then uses the received waveform to make a guess at the information sequence.
Physical channels are usually band-limited, meaning that for practical purposes the channel output waveform will not have any frequency components above some limit W Hz. Many channels are also linear (or we assume they are), so that the frequency components of the signalling waveform that are above W Hz will not influence the channel output. Because of this, we need only consider signalling waveforms that are also band-limited to W Hz. Using Nyquist sampling at a rate of 1/Dt = 2W samples/second, a signalling


Page 130
waveform defined on [0,N Dt] can be represented exactly by the discrete-time sequence [image] . The transmission of each sample a i is called a channel usage. Similarly, the channel output waveform can be represented exactly by the discrete-time sequence
Since the information sequence is effectively random, for multiple trials different signalling sequences will be produced according to some (usually discrete) distribution p(a). The channel output sequence is probabilistically related to this sequence by a channel model p(y|a).
For a fixed level of additive noise, the transmitter can communicate in an error-free fashion simply by using a very powerful signalling waveform. However, this is an uninteresting and practically expensive solution to the channel coding problem. In practice, a limit is placed on the average transmission power:
It turns out that the information rate (in bits/channel usage) that can be communicated with arbitrarily low probability of bit error, is bounded from above by the capacity C of the channel:
where the power constraint in (6.1) is enforced during the maximization. This optimal information rate was introduced by Shannon [1948], and is just the mutual information between the channel input sequence and the channel output sequence. (As a practical note, to lower the bit error rate or to use an information rate that is closer to C, we must generally use longer signalling waveforms.)
The channel coding design game essentially consists of devising encoders (ways to map information sequences to signalling sequences) and decoders (ways to guess at what the information sequence is for a given received sequence). In this book, I am mainly interested in conveying to the reader the insight and breadth of application offered by describing channel coding problems using Bayesian networks and using the probabilistic inference algorithms presented in Chapter 2 to perform decoding. For this reason, I begin this chapter by distilling out the essence of the channel coding problem and presenting a simple prototypical problem that will be the focus for the remainder of the chapter. In the prototypical problem, the transmitter sends a discrete-time binary sequence of +1's and -1's (this is called binary signalling), and each of


Page 131
these values is corrupted by additive Gaussian noise. So, the encoder maps each information sequence to a binary signalling sequence, and given a received noisy sequence, the decoder makes a guess at the binary information sequence. It turns out that the solution to this problem has far-reaching consequences in multi-level (nonbinary) coding [Imai and Hirakawa 1977], mainly due to recent proofs by Wachsmann and Huber [1995] and Forney [1997].
In Section 6.2, I show how Bayesian networks and probability propagation can be used to describe and decode Hamming codes, convolutional codes, turbocodes, serially-concatenated convolutional codes, and low-density parity-check codes. In Section 6.4, I introduce "trellis-constrained codes", which are a trellis-based generalization of all of the above codes. In Section 6.7, I present a method for speeding up iterative decoders that are implemented on serial machines.
6.1 Review: Simplifying the playing field
The real-valued signalling sequences described above are the price to pay for an efficient description of digital communication in the real world, where signal amplitudes are usually real-valued. The channel coding problem would be much simpler to pose and implement if (1) signal levels were discrete, (2) the channel model was simple, and (3) the mapping from information sequences to channel inputs was assumed to be of a relatively simple form. While this approach can simplify the problem, it can also lead to a communication rate that is far below the general capacity given in (6.2). In this section, I simplify the coding problem in the ways described above, while attempting to argue that if done properly, the simplification will lead to a communication rate that is practically very close to capacity.
6.1.1 Additive white Gaussian noise (AWGN)
A channel model that is simple and works well in practice is the AWGN channel. Additive white Gaussian noise with single-sided spectral density N0 is added to the signalling waveform to obtain the channel output waveform. Assuming the channel is bandlimited to W Hz (as described above), the decoder can apply a low-pass filter with bandwidth W Hz and sample the noisy waveform at the Nyquist rate to get a discrete-time sequence [image] . It turns out that an AWGN channel simply adds independent Gaussian noise to each input value a i where the variance of the noise is related to N0 by s2 = N0/2:


Page 132
If the decoder applies a low-pass filter with a higher bandwidth, then frequency components of the AWGN that are above W Hz will increase the effective noise on the sequence [image].
The AWGN channel leads to an appealing formulation of maximum likelihood (ML) signal detection. The log-probability density of the received sequence given the signalling sequence is
So, ML signal detection for the AWGN channel consists of finding the allowed signalling sequence a that is closest to y in Euclidean distance.
6.1.2 Capacity of an AWGN channel
For the AWGN channel, each channel output y i depends only on a i, and not any a j, j ¹ i. Consequently, the signalling distribution p(a) that will give the highest mutual information is of product form:
(This distribution allows us to stuff as much information into each a i as possible.) In this case, the capacity in (6.2) simplifies to
bits per channel usage. Note that for a product-form signalling distribution, the power limit in (6.1) becomes [image] .
It turns out that the maximum in (6.6) is obtained by a Gaussian signalling distribution with variance P (Cover and Thomas [1991]), and the capacity is


Page 133
For example, if P = 3s 2, then C = 1 bit/channel usage. For reasonable power levels, it is not possible to deterministically map C bits of information to a value a i that will have a Gaussian distribution (or even close to Gaussian). For example, try mapping 1 bit of information to a variable whose distribution is close to Gaussian!
The optimality of a Gaussian signalling distribution leads to a new type of coding concept called shaping. A signalling technique has good shape if the marginal signalling distributions are nearly Gaussian. If the signalling shape is poor, then the capacity given in (6.7) cannot be achieved no matter how good a code is used. For example, if binary signalling is used [image], then the channel capacity cannot be achieved, as shown in Section 6.1.6 and Figure 6.2.
The interplay between shaping and coding is very important. As another example, here is a method that has an excellent signalling shape, but uses a poor code. We first construct a table that maps each information vector u to a real value Cu in a way so that a uniform distribution over information vectors induces a nearly Gaussian distribution over Cu. For a given information vector u, the transmitter simply sends a constant waveform, a i = cu, i = 0, . . . , N - 1. Using this method, each marginal distribution p(a i) can be made to be as close to Gaussian as desired, by increasing N and refining the map from u to cu. However, because the waveform is constant there is no way to introduce a good code. A fruitful structure that leads to a nice mix between coding and shaping is the signal constellation.
6.1.3 Signal constellations
Since the information sequence is discrete and the signalling sequence is determined from the information sequence, the allowable set of signalling sequences is also discrete. How should we specify the set of allowed signalling sequences? One way is to require that the signalling variable at each time step be a member of a fixed signal set. Figure 6.la shows the signalling points for two signalling variables a0 and a1, where each variable can take on one of eight values. Even if a good code is used with these signalling points, the marginal signalling distributions are quite far from Gaussian and so the rate will be below capacity. Instead, consider breaking the signalling sequence into a series of groups (i.e., subspaces) containing n values each. A discrete set of values (called a constellation) is then judiciously chosen within each n-dimensional subspace in a way that leads to marginal signalling distributions that are close to Gaussian.


Page 134
Figure 6.1
Signal constellations can be used to increase the Euclidean
distance between signalling points (indicated by crosses).
(a) A naive constellation for an n = 2 signalling set with 64 points.
(b) The same 64 points can be rearranged in order to reduce the
transmission power. (c) The constellation from (b) scaled
up so that its transmission power is the same as the power
for (a)—notice that the nearest-neighbor distance has increased.
Another way to understand signal constellations is through a sphere-packing argument. Consider the 2-dimensional constellation shown in Figure 6.1 a that corresponds to the naive approach described above. For a fixed noise level in an AWGN channel, detection error falls off with distance between nearest-neighbor signals. Imagine centering a 2-dimensional sphere on each signal point as shown. Now, by trying to pack the spheres as tightly as possible, we obtain the constellation shown in Figure 6.1b. The nearest-neighbor distance has not changed, but the transmission power has decreased (since the sum of squared distances to signalling points is lower). In order to use the same power as the naive approach uses, we can now increase the Euclidean distance between nearest-neighbors as shown in Figure 6. 1c. This will increase the noise-tolerance of the system, and so increase the communication rate relative to the naive approach. For higher-dimensional constellations, this sphere-packing gain becomes more valuable. (This simple example ignores the increase in the number of nearest neighbors from 4 to 6. See Lee and Messerschmitt [1994] for more details.)
6.1.4 Linear binary codes can get us to capacity
Although the design of optimal high-dimensional constellations is straightforward in theory, it is very difficult to implement practical encoders and decoders that use these constellations. Consequently, we must approximate optimal constellations by practical ones. Ways of doing this include trellis codes (a.k.a. coset codes) [Ungerboeck 1982; Calderbank and Sloane 1987; Forney 1988], which your telephone modem probably uses. Alternatively, Wachsmann and Huber [1995] and Forney [1997] have shown that by using a technique called


Page 135
multilevel coding [Imai and Hirakawa 1977], we can achieve the capacity in (6.2) by combining several relatively simple linear binary codes. That is, optimal constellations can be well approximated if we can design appropriate linear binary codes. I refer to these new proofs to justify my focus on linear binary codes in this book.
A binary code maps each binary information vector u of length K to a binary codeword vector x of length N. The rate R of such a binary code is defined as
I will sometimes highlight the mapping by writing the codeword for u as x(u). A binary code is linear if for any u1 and u2, x(u1 Å u2) = x(u1) Å x(u2), where, "Å" indicates component-wise modulo 2 addition (0 Å 0 = 0, 0 Å 1 = 1, 1Å 0 = 1, and 1 Å 1 = 0). Note that this form of linearity is highly nonlinear in the sense of continuous algebra (where 1 + 1 = 2). In general, linear codes are easier to analyze than nonlinear ones.
Each bit in the codeword can be transmitted using binary signalling, also called binary antipodal signalling. (If the binary signal is modulated by a carrier so that it is a passband signal, it is sometimes called binary phase-shift keying (BPSK).) For x i = 1 we transmit [image] and for x i = 0 we transmit [image]. In this way, the average transmitted power is P. For an AWGN channel, we can write the probability density of channel output yi directly in terms of xi (bypassing a i):
A simple linear binary code is the repetition code. Each information bit is transmitted m times, so that R = K/mK = 1/m. Using (6.9), the probability density of channel outputs y0, . . . , y m-1 given x0 is
where the constant of proportionality does not depend on x0. By basing the decoding decision on [image] , the receiver effectively reduces the noise


Page 136
variance by a factor of 1/m. It turns out that this is a very poor code, because the suppression of noise comes at too high a cost in terms of decreasing the code rate.
6.1.5 Bit error rate (BER) and signal-to-noise ratio (E b/N 0)
For many engineering applications, the distortion value of interest is the probability P b that an information bit will be guessed incorrectly by the decoder. When analytic methods are not available for computing P b, we must resort to simulation. Often the simulation results are summarized as a point estimate called the bit error rate (BER). The BER is usually simply the observed fraction of information bit errors. When it is not possible to simulate the transmission of enough words to accurately pin down the probability of bit error, techniques such as the one described in Frey and MacKay [1998b] can be used to produce a confidence interval.
To compare the BERs of different coding schemes, we need a relatively robust measure of the noise level that each system is being exposed to. Simply stating the noise variance for an AWGN channel is not sufficient, since one system may be transmitting at a much higher power than another. Also, as shown above, performance can be improved in a trivial fashion simply by repeating signals. A reasonably robust measure of the noise level is
where P is the transmitter power, N0 is the single-sided spectral density of the AWGN, s 2 is the AWGN variance, and R is the rate of the code. E b/N0 is the ratio between the power that is transmitted per information bit, and the AWGN power. It is usually given in units of decibels (dB),
Notice that although dividing by R in (6.11) does cancel the effect of the improvement obtained trivially by repeating signals, it does not take into account the increased bandwidth needed for lower rate codes. In fact, in the next section we see that the minimum E b/N0 needed for error-free communication depends on the rate. So, when comparing one coding system to another that uses a lower rate, we must keep in mind that there is usually some way to modify the former system so as to lower its rate and at the same time lower the E b/N0 it needs to achieve error-free communication at that rate.


Page 137
6.1.6 Capacity of an AWGN channel with binary signalling
Engineering bandwidth restrictions aside, what are the communication limits for an AWGN channel when we use binary + 1/-1 signalling (binary antipodal signalling with P = 1)? (Without loss of generality, we will assume that P = 1.) The answer to this question depends on whether we are willing to tolerate a certain non-vanishing BER. Before considering the non-vanishing BER case in the next section, I will address the simpler case of a vanishing BER. More specifically, what is the minimum E b/N0 needed to communicate error-free using a rate R code on an AWGN channel with binary signalling?
The mutual information between the channel input a i and the channel output y i at time step i gives the number of bits that can be communicated per channel usage on average. For an AWGN channel with binary signalling, the mutual information as a function of the noise variance is
The first term is the entropy of y i given a i, which is just the entropy of a Gaussian distribution, 0.5 log2(2ps 2e). Since p(y i) is a mixture of two Gaussians, the second term is
which can be approximated quite well using a Monte Carlo method. In this fashion, it is possible to obtain a good estimate of M(s 2).
To communicate error-free, the rate of the code must be less than the mutual information between the channel input and the channel output: R < M(s 2) [Shannon 1948]. Inserting s 2 = 1/(2RE b/N0) (see (6.11)) into this inequality, we get [image] . After rearrangement, we have


Page 138
Figure 6.2
(a) The minimum E b/N0 needed for error-free communication with a rate
R code, over an AWGN channel using binary signalling and optimal
(Gaussian) signalling. (b) The minimum achievable BER as a function
of E b/N0 for several different code rates using binary signalling.
This bound (based on an interpolated inverse of a Monte Carlo estimate of M(s 2)) is shown in 6.2a, along with the minimum E b/N0 required by optimal (Gaussian) signalling (see Section 6.1.2).
For example, an R = 1/2 code requires E b/N0 > 0.2 dB. To communicate error-free without coding (R = 1), an infinite E b/N0 is needed.
A standard result from information theory is that regardless of rate, an E b/N0 of at least loge2 = -1.5917 dB) is required for error-free communication [Cover and Thomas 1991]. This limit is apparent from the convergence of the curve as R ® 0.
6.1.7 Achievable BER for an AWGN channel with binary signalling
If we are willing to tolerate a certain non-vanishing BER while using a rate R code, it turns out we can use a lower E b/N0. than described in the previous section. One way to pose the problem for this scenario is: For an optimal code with rate R and a specified BER, what is the minimum required E b/N0? We can think of this as a two-stage problem. First, we find a shorter representation for the information vector. This representation will obviously be lossy, since a uniformly random vector of information bits cannot be represented losslessly on average by a shorter binary vector. Second, we use a new optimal code


Page 139
to communicate this shorter representation error-free over the channel with the largest tolerable noise variance. Since the representation is shorter than the information vector, the new code rate R' will be lower than the old one: R' < R. So, the tolerable noise variance for error-free communication of the lossy representation will be higher than the tolerable noise variance for error-free communication of the information vector.
We would like to use a representation that is as short as possible, so that R' will be as low as possible and the tolerable noise variance will be as large as possible. However, shorter representations are also more lossy and will lead to higher BERs. What is the minimum ratio between the length of the representation and the length of the information vector, such that the error rate does not rise above the specified BER? It turns out that the minimum ratio is just the mutual information between a uniformly random bit and its noisy duplicate, where the probability that the value of the duplicate is flipped is BER. (This can be viewed as a result of rate-distortion theory applied to a Bernoulli source [Cover and Thomas 1991].) This mutual information is
The new code rate is
For a specified R and BER, we can compute R', determine the maximum tolerable noise variance s2 = M-1(R'), and compute the minimum E b/N0 = 1/(2s 2R) (to compute E b/N0, we use the original R, not R'). Figure 6.2b shows the achievable BER as a function of E b/N0 for several different rates. For each rate, the value for E b/N0 at which the BER converges to zero is the same as the value shown in Figure 6.2a. These achievable BER curves are used as guides for ascertaining the performances of codes and decoders later in this chapter.
6.2 Graphical models for error correction: Turbocodes, low-density parity-check codes and more
A critical component of a channel coding system is the decoder. Even if the code gives excellent performance when optimal decoding is used, if there is no way to implement a practical decoder that gives similar performance, it is not clear that the code is of any use. Channel decoders can be broken into two classes: algebraic and probabilistic. Algebraic decoders for binary codes usually quantize the channel output to two or three levels. The received vector


Page 140
y is interpreted as a copy of the binary codeword vector x, with some of the bits flipped. Alternatively, received values that are highly ambiguous (e.g., the value 0.1 when +1/- 1 binary signalling is used) are considered as erasures—i.e., the corresponding bit in y is assumed to be unknown. In both cases, decoding is a matter of using linear algebra (in a finite field) to find the binary codeword vector x that is closest to y in Hamming distance (dimensions that are erased are ignored). There are many techniques for algebraic decoding [Lin and Costello 1983; Blahut 1990; Wicker 1995] and algebraic decoders usually take advantage of special structure that is built into the code to make decoding easier. However, it is obvious that by using such a coarsely quantized form of the channel output, these decoders are suboptimal. For example, the value 0.1 from above does provide some evidence that a signal value of +1 was sent.
Probabilistic decoders are designed to make as much use as is practically possible of the real-valued unquantized channel output. The goal of probabilistic decoding is either maximum likelihood (ML) information sequence detection, or maximum a posteriori (MAP) information bit detection:
Obviously, ML sequence detection minimizes the word error rate (we usually assume that all words are equally likely a priori), while MAP bit detection minimizes the BER. So, by definition, optimal probabilistic decoders are superior to optimal algebraic decoders. However, can we implement useful probabilistic decoders? The success of algebraic decoders is due to the way they take advantage of the algebraic structure of a code. Is there an analogous structure that probabilistic decoders can use? In this section, I show how graphical models can be used to describe probabilistic structure for channel codes and how the inference algorithms that make use of this structure can be used for probabilistic decoding.
6.2.1 Hamming codes
Hamming codes are an extension of the notion of adding a single parity-check bit to a vector of information bits. Instead of adding a single bit, multiple bits are added and each of these parity-check bits depends on a different subset of the information bits. Hamming developed these codes with a special algebraic structure in mind. Consequently, they are really meant for binary channels where the noise consists of randomly flipping bits. However, Hamming codes


Page 141
are short and easy to describe, so they make a nice toy example for the purpose of illustrating probabilistic decoding.
An (N, K) Hamming code takes a binary information vector of length K and produces a binary codeword of length N. For an integer m ³ 2, N and K must satisfy N = 2m - 1 and K = 2m - m - 1. The Bayesian network for a K = 4, N = 7 rate 4/7 Hamming code is shown in Figure 6.3a. The algebraic structure of this code can be cast in the form of the conditional probabilities that specify the Bayesian network. Assuming the information bits are uniformly random, we have P(u k = 1) = P(u k = 0) = 0.5, k = 0, 1, 2, 3. Codeword bits 0 to 3 are direct copies of the information bits: P(x k|u k) = d (x k|u k), k = 0, 1, 2, 3. Codeword bits 4 to 6 are parity-check bits:
Assuming binary antipodal signalling with power P over an AWGN channel, the conditional channel probabilities p(y i|x i), i = 0, 1, 2, 3, 4, 5, 6 are given by (6.9), where s2 is related to E b/N0 by (6.11).
This code is small enough that we can compute the MAP bit values in (6.18) exactly using Bayes' rule. The BER-E b/N0 curve for MAP bit decoding and the achievable BER (see Section 6.1.7) at rate 4/7 are shown in Figure 6.3b. Although there is an 8 dB gap between these curves at a BER of 10-6, the MAP decoder gives a significant improvement of 2dB over uncoded transmission.
By making hard decisions for the channel outputs (calling a value below 0 a ''0'' and calling a value above 0 a "1"), an algebraic decoder can be used. This decoder applies a parity-check matrix to the received binary word in order to try to locate any errors. (See Lin and Costello [1983] for details.) In this fashion, it can correct up to one bit error per codeword. The curve for algebraic decoding is also shown in Figure 6.3b. Algebraic decoding gives an improvement of only 0.5 dB over uncoded transmission at a BER of 10-6. Although this may seem surprising, keep in mind that the receiver for the uncoded transmission is allowed to average the channel output to reduce the effective noise (see Section 6.1.4) 7/4 times longer than the receiver for the algebraic decoder.
One way to approximate the probabilities P(u k|y) used for MAP bit decoding is to apply the probability propagation inference algorithm (Section 2.1 ) to the Bayesian network shown in Figure 6.3a [Frey and MacKay 1998a]. Probability propagation is only approximate in this case because the network is multiply-connected or "loopy" (e.g., u0-x4-u1-x5-u0. Once a channel output vector y is observed, propagation begins by sending a message from y k


Page 142
Figure 6.3
(a) The Bayesian network for a K = 4, N = 7 Hamming code. (b) BER
performance for the maximum a priori decoder, an iterative probability
propagation decoder, the optimal algebraic decoder, and uncoded
transmission.
to xk for k = 0,1,2,3,4,5,6. Then, a message is sent from x k to u k for k = 0, 1, 2, 3. An iteration now begins by sending messages from the information variables u0, u1, u2, u3 to the parity-check variables x4, x 5, x6 in parallel. The iteration finishes by sending messages from the parity-check variables back to the information variables in parallel. Each time an iteration is completed, new estimates of P(u k|y) for k = 0, 1, 2, 3 are obtained. The curve for probability propagation decoding using 5 iterations is shown in Figure 6.3b. It is quite close to the MAP decoder, and significantly superior to the algebraic decoder.
For this simple Hamming code, the complexities of the probability propagation decoder and the MAP decoder are comparable. However, the similarity in performance between these two decoders raises the question: "Can probability propagation decoders or other approximate decoder give performances comparable to MAP decoding in cases where MAP decoding is computationally intractable?" Before exploring a variety of systems where probability propagation in multiply-connected networks gives surprisingly good results, I will review convolutional codes, whose graphical models are essentially singly-connected chains. For these networks, the probability propagation algorithm is exact and it reduces to the well-known forward-backward algorithm [Baum and Petrie 1966] (a.k.a. BCJR algorithm [Bahl et al. 1974]). (For more on this connection, see Smyth, Heckerman and Jordan [1997].)


Page 143
Figure 6.4
The linear feedback shift register (LFSR) configurations
for rate 1/2 convolutional codes with maximum dmin.
(a) A systematic nonrecursive convolutional code (dmin = 7).
(b) A nonsystematic nonrecursive convolutional code (dmin = 10).
(c) A systematic recursive convolutional code (dmin = 10).
6.2.2 Convolutional codes
Convolutional codes are produced by driving a finite state machine with information bits. The outputs of the finite state machine (which may include copies of the inputs) are then used as codeword bits. A code for which each information bit appears as a codeword bit is called systematic. Typically, linear convolutional codes are used, and any code in this class can be represented by a linear feedback shift register (LFSR). An example of a systematic code of this type with a memory of 7 bits is shown in Figure 6.4a. Each box represents a 1-bit memory element and D is a delay operator: D nuk = u k-n. In this example, there is no feedback from the shift register to its input; a convolutional code of this type is called nonrecursive. An output is produced by


Page 144
adding (modulo 2) values "tapped" from the memory chain. The output taps for this rate 1/2 systematic nonrecursive convolutional code were chosen to maximize the minimum distance dmin between codewords [Lin and Costello 1983]. For this code, dmin = 7, meaning that the codeword vectors for any two information vectors will differ in at least 7 places. Using the delay operator, this code can be described by the following two equations:
where G(D) is called the generator polynomial. This polynomial is often expressed in octal form by letting the coefficient of D O be the least significant bit and the coefficient of D 7 be the most significant bit. In this case the octal representation is 3538.
Since dmin plays the central role in determining the error-correcting capabilities of a code at high signal-to-noise ratio E b/N0, we would like to use codes that have large dmin. One way to obtain a greater dmin for convolutional codes is to use a larger memory. However, it turns out that decoding complexity increases exponentially with the size of the memory. In fact, it is possible to increase the minimum distance of any systematic nonrecursive convolutional code without using more memory. Figure 6.4b shows a rate 1/2 nonsystematic nonrecursive convolutional code that has dmin = 10. (The two sets of output taps that maximize dmin were found using a method described in Lin and Costello [1983].) This code can be described as follows:
For a nonsystematic convolutional code, there are two generator polynomials corresponding to the two sets of output taps. For this code, the octal representation is (2478, 3718).
Although the performance of the nonsystematic code described above is better than the systematic one at high E b/N0, it is the other way around for values of E b/N0 near the Shannon limit. Berrou and Glavieux [1996] have argued that a nice compromise between these codes is a systematic recursive convolutional code. The code in Figure 6.4b can be converted to a systematic code by taking one set of the output taps (either one will do) and using them as feedback to the input of the shift register, making a LFSR. If we do this with the upper set of taps, we obtain the rate 1/2 systematic recursive convolutional code shown in Figure 6.4c. This code can be described by the following two equations:


Page 145
The second equation is to be interpreted as
which can be derived from the figure.1 The former expression allows us to retain the G(D) notation, which in this case is 2478/3718. From the point of view of linear algebra, we have obtained this new code simply by dividing G1 (D) and G 2(D) from above by G1 (D). It can be shown that this operation does not change the algebraic structure of the code. For example, the new code has dmin = 10 as before. However, as we saw in the previous section, there is more to channel coding than algebraic structure. It turns out that this systematic recursive code performs better than the above nonsystematic nonrecursive code at low E b/N0.
6.2.3 Decoding convolutional codes by probability propagation
Bayesian networks for nonsystematic and systematic convolutional codes are shown in Figures 6.5a and 6.5d. In the former case, both codeword bits at stage k depend on the encoder state as well as the information bit, whereas in the latter case, one codeword bit is simply a direct copy of the information bit. Notice that because of the dependency of at least one codeword bit at stage k on the encoder state and the information bit, these networks are not singly-connected. However, they can be converted to singly-connected networks in the following way. By duplicating the information bits, we obtain the networks shown in Figures 6.5b and 6.5e (see Section 2.1.4). By grouping each state variable with one of these duplicates as shown by dashed loops, we obtain the singly-connected networks shown in Figures 6.5c and 6.5f (see Section 2.1.4).
In the new networks, each state variable actually contains a copy of the current information bit. We can interpret each state variable as a binary number whose least significant bit (LSB) is a copy of the current information bit and whose most significant bit (MSB) is the oldest value in the LFSR (i.e., the value in the memory element that appears on the far right in the LFSRs shown
1 In fact, this representation is algebraically consistent. We can, for example, multiply the numerator and the denominator in (6.22) by a polynomial in D without changing the set of output sequences that the LFSR can produce. See Wicker [1995] for a textbook treatment.


Page 146
Figure 6.5
The multiply-connected Bayesian network (a) for a nonsystematic
convolutional code can be converted to a singly-connected network by
duplicating the information variables (b) and then grouping together
information variables and state variables (c). The multiply-connected
network for a systematic convolutional code can be converted to a
singly connected one (d)—(f).
in Figure 6.4). Let s k/2 be the binary number obtained by cutting off the LSB of s k, and let s k % 2 be the value of the LSB of s k. Let ¦(s k-1) be the binary number obtained by cutting off the MSB of s k-1 and replacing the LSB of sk = 1 with the value of the LFSR feedback bit obtained by adding (modulo 2) the values of the bits in s k-1/2 that correspond to the LFSR feedback taps. So, ¦( s k-1) is the value of the new state at stage k, excluding information bit u k. Finally, let g(s k) be the bit obtained by adding (modulo 2) the values of the bits in sk that correspond to the LFSR output taps. If there are two sets of taps, then there will be two output functions g1 (s k) and g2 (s k).
Now, we can specify the conditional probabilities for the convolutional code Bayesian networks. For the sake of brevity, I will consider only the systematic code shown in Figure 6.5f. Assuming the information bits are uniformly random, we have P(u k = 1) = P(u k = 0) = 0.5, k = 0, . . . , K - 1. The state transition probabilities are


Page 147
Figure 6.6
The computation of P(uk|y), P(sk|y), and P(xi|y) by probability propagation
using the forward-backward message-passing schedule, which minimizes the
total number of messages passed. Arrows represent messages in transit,
whereas solid dots represent messages waiting to be sent.
where we assume s-1 = 0 to initialize the chain. The codeword bit probabilities are
Assuming binary antipodal signalling with power P over an AWGN channel, the conditional channel probabilities p(y i|x i), i = 0, . . . , 2K - 1 are given by (6.9), where s 2 is related to E b/N0 by (6.11).
Probability propagation in the singly-connected Bayesian networks for convolutional codes can be used to compute the a posteriori bit probabilities P(u k|y) exactly. The MAP values [image] can be obtained by applying a threshold of 0.5 to these probabilities. Although the probability messages can be passed in any order, the forward-backward message-passing schedule gives the lowest number of total messages passed, and so it is most appropriate for


Page 148
Figure 6.7
Performances of the following 7-bit memory
LFSR convolutional codes with maximum dmin:
systematic nonrecursive convolutional code
(SNRCC), systematic recursive convolutional
code (SRCC), nonsystematic nonrecursive
convolutional code (NSNRCC).
decoding on a serial machine. Figure 6.6 shows how messages are passed according to this schedule in the Bayesian network for a simple systematic convolutional code. First, probability messages are propagated from the observed channel output variables (crossed vertices) to the "backbone" of the chain (the state variables). Then, the messages are buffered as shown. (See Section 2.1.3 for an explanation of buffered messages in probability propagation.) Pictorially, when a message arrives at a vertex on an edge, but is buffered and not propagated on to the other neighbors, I draw a small dot adjacent to each of the other edges. Each of these dots can be turned into an arrow (indicating a message is being passed) at any time. Next messages are passed forward along the chain, and then backward along the chain. Finally, messages are propagated to the information bits and to the codeword bits. (It is not necessary to propagate probabilities to the observed variables, since P(y i|y) is trivial to compute.) Notice that this algorithm computes P(u k|y), P(s k|y), and P(x i|y). If all we need are the information bit probabilities P(u k|y), then it is not necessary to propagate the last set of messages shown in the figure.
Figure 6.7 shows the performances of the three convolutional codes described above. The systematic nonrecursive convolutional code has a BER that is significantly higher than the BERs for the other two codes at reasonably high E b/N0. The nonsystematic nonrecursive convolutional code and the systematic recursive convolutional codes have similar BERs, except for at low


Page 149
E b/N0, where the systematic code has a significantly lower BER. The software used to obtain these results is available by following pointers at http://mitpress.mit.edu/book-home.tcl?isbn=026206202X.
6.2.4 Turbocodes: parallel concatenated convolutional codes
Although the convolutional codes and decoder described above give roughly a 5.7 dB improvement over uncoded transmission at a BER of 10-5, they are still roughly 3.7 dB from Shannon's limit at this BER. Up until the last few years, a serially-concatenated Reed-Solomon convolutional code [Lin and Costello 1983] was considered to be the state of the art. At a BER of 10-5, this system is roughly 2.3 dB from Shannon's limit. However, in 1993, Berrou, Glavieux, and Thitimajshima introduced the turbocode and the practical iterative turbodecoding algorithm. Their system was roughly 0.5 dB from Shannon's limit at a BER of 10-5. Also, these binary codes have been successfully combined with multi-level coding to obtain bandwidth-efficient coding within 0.7 dB of Shannon's limit [Wachsmann and Huber 1995].
The original presentation of turbocodes lacked a principled framework. For example, it was not at all clear how decoding should proceed when there were three or more constituent convolutional codes instead of two [Divsalar and Pollara 1995]. However, it turns out that the turbocode can be concisely described as a multiply-connected graphical model, and that the turbodecoding algorithm is just probability propagation in this model [Frey and Kschischang 1997; Frey and MacKay 1998a; Kschischang and Frey 1998; McEliece, MacKay and Cheng 1998]. This general graphical model framework makes it easier to describe new codes and their corresponding iterative decoding algorithms. For example, decoding a turbocode that has three constituent convolutional codes is just a matter of propagating probabilities in the corresponding graphical model.
Figure 6.8a shows the Bayesian network for a rate 1/2 turbocode. For a given information vector, the codeword consists of the concatenation of two constituent convolutional codewords, each of which is based based on a different permutation in the order of the information bits. The subnetwork indicated by a dashed loop is essentially the same as the network for the systematic convolutional code described above. The only difference is that every second LFSR output is left off, for a reason given below. The information bits are also fed into the upper convolutional encoder, but in permuted order. Every second LFSR output of the upper code is also left off. By leaving off every second LFSR output in both constituent codes, the total number of codeword bits is twice the number of information bits, so the rate is 1/2.
Once the channel output y for an encoded information vector is observed, an approximate inference method can be used to approximate P(uk|y) and


Page 150
Figure 6.8
(a) The Bayesian network for a K = 6, N = 12 rate 1/2 turbocode.
(b) The performance of a K = 65, 536 rate 1/2 turbocode using
18 iterations of turbodecoding.
perform approximate MAP bit decoding. Variational inference has been used to decode a similar class of "low-density parity-check" codes [MacKay 1997a] (see Section 6.2.5), but there it was found that probability propagation in the multiply-connected graphical model gave better results. Here, I focus on decoding by probability propagation. Figure 6.8b shows the performance of the probability propagation decoder for a K = 65, 536 rate 1/2 turbocode with a randomly drawn permuter. The software used to obtain these results is available by following pointers at http://mitpress.mit.edu/book-home.tcl?isbn=026206202X.
Each (identical) constituent convolutional code uses a 4-bit LFSR with polynomials (21/37)8. Although at low Eb/No the turbocode gives a BER that is significantly higher than the BER for uncoded transmission, the turbocode curve drops below a BER of 10-5 at less than 0.5 dB from Shannon's limit. Berrou and Glavieux suggest that for very low BER performance (say 10-10), the permuter should be designed to maximize dmin [Berrou and Glavieux 1996]. I have found that for BERs at or above 10-5, a randomly drawn permuter typically works fine.
Since the turbocode network is multiply-connected (there are cycles when the edge directions are ignored), we must specify a message-passing schedule in order to decode by probability propagation. That is, the order in which messages are passed can affect the final result as well as the rate of convergence


Page 151
Figure 6.9
The message-passing schedule corresponding
to the standard turbodecoding algorithm.
to a good decoding solution. Since the network is multiply-connected, we must also specify when to stop passing messages, since otherwise they would propagate indefinitely. Figure 6.9 shows how messages are passed up to the end of the first iteration of turbodecoding. First, messages are passed from the channel output variables (crossed vertices) to the state variables of both constituent codes. Assuming we are only interested in estimating P(u k|y), we can now ignore the channel output variables and the codeword variables. The simplified network with buffered messages waiting to be sent is shown in the upper-right picture in Figure 6.9.
Next, messages are passed from the information variables to the state variables of one of the constituent codes. This chain is processed in the forward-backward manner and then messages are propagated to the information variables. Messages are then passed to the state variables of the other constituent code (these are called ''extrinsic information'' in Berrou and Glavieux [1996]).


Page 152
Once the second chain has been processed in the forward-backward manner, messages are propagated back to the information variables, as shown in the lower-right picture in Figure 6.9. This completes the first iteration of turbodecoding. Messages are then propagated from the information variables back to the first constituent code chain, and so on. The series of 16 pictures outlined by a dashed rectangle in Figure 6.9 shows how messages are passed during one complete iteration of turbodecoding. (Note that after the first iteration, there aren't any buffered messages in the first picture within the dashed rectangle. The buffered messages in this picture are due to the initial observations.)
6.2.5 Serially-concatenated convolutional codes, low-density parity-check codes, and product codes
It turns out that many of the iterative decoding algorithms for a variety of codes can be viewed as probability propagation in the corresponding graphical models for the codes [Frey and Kschischang 1997]. Figure 6.10a shows the Bayesian network for a serially-concatenated convolutional code [Benedetto and Montorsi 1996b]. The information bits are first encoded using the upper convolutional code. The generated codeword bits x1 are then permuted and fed into a second convolutional encoder, whose output bits x2 are transmitted over the channel. The iterative decoding algorithm introduced by Benedetto and Montorsi [ 1996a] was presented without reference to any of the literature on probability propagation. However, their iterative decoding algorithm is in fact probability propagation in the corresponding Bayesian network. After observing the channel output y, the decoder propagates messages from y to the lower chain. Then, messages are propagated forward and backward along the lower chain before being passed to the upper chain. The upper chain is processed and then messages are passed back to the lower chain, and so on.
The analytic upper bounds on the error rate for maximum-likelihood decoded serially-concatenated convolutional codes are superior to those for turbocodes [Benedetto et al. 1997]. However, it is not clear that these theoretical bounds are of any practical value. First, the bounds are based on the average performance over all possible permuters. Suppose that on average 1 in every 1000 permuters gives a very poor code that when ML-decoded gives a BER of 0.1. Further, suppose that the other permuters give codes that when ML-decoded give BERs of 10-10. If we randomly pick a permuter, we are very likely to get a code that gives a BER of 10-10 However, the average performance over all permuters is 0.001 × 0.1 + 0.999 × 10-10 » 10-4. In this way, the average performance over permuters can be misleading. Second, since ML decoding is intractable, in practice we must use a suboptimal decoder, such as probability propagation. Even if the ML-curve for one code is superior to


Page 153
Figure 6.10
The Bayesian networks for (a) a K = 3, N = 9 rate 1/3 serially-concatenated
convolutional code; (b) a K = 2, N = 6 rate 1/3 low-density parity-check code;
and (c) a K = 3, N = 9 rate 1/3 product code. In the latter two figures, the
channel output variables have been left out for visual clarity.
that of another code, the performance of the practical iterative decoder may be inferior.
It is suggested by Benedetto et al. [1997] that for short block lengths (say, K = 200) serially-concatenated convolutional codes give better performance than turbocodes, when iterative decoding is used. However, for short block lengths, it is not at all clear that either of these codes performs better than sequential decoding [Lin and Costello 1983] with a convolutional code with large memory.
Figure 6.10b shows the factor graph for a low-density parity-check code [Gallager 1963; Tanner 1981; MacKay and Neal 1996]. Unlike convolutional codes, turbocodes, and serially-concatenated convolutional codes, these codes do not have a simple directed generative structure that maps information bits to codeword bits. So, the factor graph is a more suitable graphical model for describing low-density parity-check codes. Notice that the relationship between the codeword bits x and the channel output signals y can still be signified with directed edges. Low-density parity-check codes were largely forgotten in the channel coding community for roughly 35 years, probably due to the computationally intensive encoder and decoder that Gallager proposed. However,


Page 154
it turns out that these codes have excellent theoretical performance [MacKay 1997a] and that the iterative decoder proposed by Gallager is in fact equivalent to probability propagation in the network shown above. In these codes, each parity-check function ¦i(×) evaluates to 1 if the codeword bits {x j}jÎQ i to which ¦i(×) is connected have even parity:
where Å is addition modulo 2 (XOR). The term "low-density" refers to the fact that each parity-check function is connected to very few codeword bits (a vanishing fraction, as N ® ¥). (Notice that since this network is parity-check oriented and does not show how an information vector leads to a codeword, it appears an encoder must use a pre-derived generator matrix and encode the K information bits in [image] time.) The iterative decoder passes messages between the parity-check functions and the codeword bit variables. In fact, each parity-check function can be represented by a simple convolutional code network with one bit of memory. As a result, the decoder is simpler than the iterative decoder for turbocodes. However, it appears they do not come as close to Shannon's limit as do turbocodes for rates of 1/3 and 1/2 [MacKay and Neal 1996].
Figure 6.10c shows the factor graph for a product code. In this graph, the parity-check functions check rows and columns of codeword bits. Recently proposed iterative decoders for product codes [Lodge et al. 1993; Hagenauer, Offer and Papke 1996] can be viewed as probability propagation in the corresponding networks. As with the low-density parity-check code, the decoder iteratively passes messages between the parity-check functions and the codeword bit variables.
6.3 "A code by any other network would not decode as sweetly"
In the previous section, I presented the graphical models for a variety of codes whose iterative decoding algorithms can be viewed as probability propagation in the corresponding graphical models. Can we use this perspective to propose new codes and derive new iterative decoders? Partly, the answer is "yes". However, we cannot expect to obtain good results simply by tossing the ingredients of a graphical model into a bag and shaking. First, we want the resulting code to give excellent performance if ML decoding is used. Second, we want the resulting code to give good results when decoded by probability propagation, which is only an approximation to maximum likelihood decoding.


Page 155
For example, since any linear binary code can be described by a set of parity-check equations, it may seem that a fruitful approach to getting closer to capacity is to simply find a good code (e.g., a random linear code), write down its parity-check equations, construct the corresponding graphical model, and then decode it using probability propagation. However, in general the parity-check network will be multiply-connected. Since probability propagation is only approximate in such networks, the performance of the decoder will depend heavily on which set of linearly independent equations is used. Graphical operations, such as grouping parity-check functions together, will also heavily influence the decoder's performance.
In order to obtain a good coding system, we need to simultaneously find a good code and a corresponding graphical model that gives good performance when decoded by probability propagation. Keeping these issues in mind, a wise approach to proposing new code networks is to incrementally generalize previous work.
6.4 Trellis-constrained codes (TCCs)
In this section, I present "trellis-constrained" codes [Frey and MacKay 1998b] that can be viewed as a trellis-based generalization of turbocodes, serially-concatenated convolutional codes, low-density parity-check codes, and product codes.
The term trellis was introduced by Forney [ 1973] and refers to a diagram that explicitly shows the values of a discrete state variable at each time step and the allowed state transitions. A trellis is more general than a LFSR, since in a trellis the state transitions and even the number of states may vary with time. (Also, a trellis can represent a nonlinear code.) Figure 6.1 la shows the trellis for the first 4 time steps of a rate 1/2 systematic recursive convolutional code with LFSR polynomials (5/7)8, i.e., x2k = u k and x2k+1 = (1 +D2)u k /1 + D + D2). The number of levels of the state variables (black discs) corresponds to the memory of the LFSR, and in this case there are 2 bits of memory. Each branch in the trellis indicates an allowed state transition, and the corresponding branch variable values (in this case the LFSR outputs x2k, x2k+1) are written beside each branch. Figure 6.11b shows the corresponding Bayesian network. In the Bayesian network, the branch variables are functions of the state alone, so each state variable has 8 levels instead of the 4 levels used in the trellis.
Figure 6.1 2a shows the factor graph for a trellis with 6 segments. For linear codes, any two adjacent state variables (large white discs) and the corresponding codeword bits (small white discs) must satisfy a linear set of equations, so the local transition function is represented by the small black discs with a


Page 156
Figure 6.11
(a) shows the trellis for a simple rate 1/2 systematic recursive
convolutional code. Each branch indicates an allowed state
transition, and the corresponding pair of output bits are written
beside the branch. (b) shows the corresponding Bayesian
network, which requires one extra bit of state so that the outputs
can be determined directly from the state variables.
"+" inside. Each function evaluates to 1 if its neighboring variables satisfy the local set of linear equations, and to 0 otherwise. The global function is equal to the product of the local functions. A given configuration of the codeword bits is a codeword if the global function evaluates to 1 for some configuration of the state variables.
Figures 6.12b and c show the factor graphs for a simple turbocode and a low-density parity-check code. In the latter case, each of the six trellises is a simple parity-check trellis that enforces even parity on its six codeword bits. These codes can be viewed as a set of codeword bits that must satisfy multiple trellis constraints. A codeword in a trellis-constrained code (TCC) must simultaneously be a codeword of multiple constituent trellises. If ¦t(x) is the constituent codeword indicator function for trellis t Î {1, . . . , T}, the global codeword indicator function is


Page 157
Figure 6.12
A codeword in a trellis-constrained code must simultaneously
be a codeword of all the constituent trellises, with the codeword
bits reordered. The factor graphs are shown for (a) a single
convolutional code; (b) a turbocode, where each parity bit is a
codeword bit of only one trellis; (c) a low-density parity-check
code; (d) a homogeneous trellis-constrained code; and (e) a
ring-connected trellis-constrained code. The small
unfilled discs represent codeword bits.
Each constituent indicator function is given by a product of the local functions within the corresponding trellis. Usually, the codeword bits interact with the constituent trellises through permuters that rearrange the order of the codeword bits. Figure 6.13 pictorially depicts this scheme.
For the turbocode in Figure 6.12b, there are two constituent functions, ¦1(×) and ¦2(×). ¦1(×) indicates that the upper row of codeword bits are valid output from a convolutional encoder with the middle row of codeword bits as input. ¦1(×) does not directly place any restrictions on the lower row of codeword bits, so it effectively only checks 3/4 of the codeword bits. ¦2(×) indicates that the lower row of codeword bits are valid output from a convolutional encoder with the middle row of codeword bits as input. In contrast to ¦1(×),¦2(×) does not place any restrictions on the upper row of codeword bits.
The rate R of a TCC is related to the rates of the constituent trellises R t, t = 1, . . . , T in a simple way. If n t is the fraction of codeword bits that


Page 158
Figure 6.13
A trellis-constrained code is the intersection of
the codes for multiple trellises, where the codeword
symbols of each constituent trellis code are permuted.
trellis t checks, then trellis t removes at most (1 - R t)n t N binary degrees of freedom from the code. It may remove a small number less if some of its constraint equations are linearly dependent on those given by other constituent trellises. For large, randomly generated permuters this effect is quite small, so it can be ignored when computing rates in the remainder of this section. (As a result, the actual rates may be slightly higher than the given rates.) The total number of binary degrees of freedom left after all trellis constraints are satisfied is [image] , so the rate of the TCC is
From this equation, it is easy to verify that the turbocode in Figure 6.12b has rate 1/2 and that the low-density parity-check code in Figure 6.12c also has rate 1/2. (Note that a k-bit parity-check trellis has rate (k - 1)/k.)
Unlike encoding turbocodes and serially-concatenated convolutional codes, encoding a general TCC takes quadratic time in N. In a general TCC, we can designate a subset of the codeword bits as the systematic bits of the entire code and then use Gaussian elimination to compute a generator matrix (once only). Using this method for encoding requires R(1 - R)N2 binary operations.
Decoding a TCC involves performing the forward-backward algorithm for each trellis and exchanging information between trellises in the fashion specified by the probability propagation algorithm. The constituent trellises may be processed in parallel or sequentially.
Next, I present two new families of TCCs and show that they perform in the same regime as do turbocodes and low-density parity-check codes.


Page 159
6.4.1 Homogeneous trellis-constrained codes
In a turbocode, the constituent trellises share only a systematic subset of their codeword bits. The other parity bits of each constituent encoder are not constrained by the other encoders. Figure 6.12d shows the factor graph for a simple homogeneous TCC with T = 2, in which all of the bits are constrained by each constituent trellis. From the general rate formula in (6.28), we see that the rate for a homogeneous turbocode is
where [image] .
One difference between the homogeneous TCC and the turbocode is that the rate of a homogeneous TCC decreases directly with the number of trellises T. In the simulations discussed below, I used T = 2 and R1 = R2 = Ravg = 3/4 to get R = 1/2. To obtain the same rate with T = 3 would require Ravg = 5/6. In contrast, the rate for a turbocode varies roughly inversely with T. A rate 1/2 turbocode with T = 3 can be obtained with R1 = R2 = R3 = 3/4. Another difference is that the permuter length of a homogeneous TCC is N, whereas for a turbocode, the permuter length is RN.
The trellises in a homogeneous TCC share all their bits, so we can't simply encode by dividing the bits in each constituent trellis into a systematic set and a parity set and running a linear-time encoding method for each trellis, as is possible in a turbocode. Instead, a previously computed generator matrix is applied to a previously selected systematic subset of codeword bits, which takes R(1 - R)N2 binary operations.
The iterative decoder processes each constituent trellis using the forward-backward algorithm, and passes "extrinsic information" between the trellises in the manner specified by the probability propagation algorithm. For two trellises, the decoding schedule is straightforward. For T > 2, different decoding schedules are possible. The trellises may be processed sequentially, in which case the current trellis uses the most recently computed probabilities produced by the other trellises. Alternatively, the trellises may be processed in parallel, in which case the current trellis uses the probabilities produced by the other trellises in the previous decoding iteration.
For the sake of gaining insight into these new compound codes and the behavior of their iterative decoders, decoding continues until a codeword is found or a large number of iterations have been performed. After each decoding iteration, the current bit-wise MAP estimates are used to determine whether a valid codeword has been found, in which case the iterative procedure is terminated. If 100 iterations are completed without finding a codeword, the block is


Page 160
labeled a decoding failure. Notice that given the factor graph of a code, determining that a codeword is valid is simply a matter of checking that all the local functions evaluate to 1.
Using Monte Carlo, I estimated the performance of an N = 131,072, T = 2 homogeneous TCC with R1 = R2 = 3/4, giving R = 1/2. (See Frey and MacKay [1998b] for a description of how the BER confidence intervals were computed.) Each rate 3/4 trellis was obtained by shortening every fifth bit of a rate 4/5 nonsystematic convolutional code with maximum dmin. (The generator polynomials for this code are given in Daut, Modestino and Wismer [1982] and are (32, 4, 22, 15, 17)octal.) Figure 6.14 shows the performance of this homogeneous TCC, relative to the turbocode introduced by Berrou et al. [1993] and the best rate 1/2, N = 65,389 low-density parity-check code published to date [MacKay and Neal 1996]. Although it does not perform as well as the turbocode, it performs significantly better than the low-density parity-check code. It seems likely that there is room for improvement here, since the set of generator polynomials that gave maximum dmin were used and this is probably not the best choice. (Keep in mind, however, that the performance of a homogeneous TCC is not necessarily governed by the same constituent trellis properties that govern the performance of a turbocode.) Of significant importance compared to turbocodes, I have observed that this homogeneous TCC has higher-weight low-weight codewords.
6.4.2 Ring-connected trellis-constrained codes
Figure 6.12e shows the factor graph for a simple ring-connected TCC with T = 3. This code can be viewed as a serially-concatenated convolutional code [Benedetto and Montorsi 1996b; Benedetto and Montorsi 1996a] in which some of the output bits are constrained to be equal to some of the input bits. The factor graph thus forms a ring of connected trellises. In the ring-connected TCC shown, each constituent trellis checks exactly 2/T of the codeword bits. From the general rate formula in (6.28), the rate for such a ring-connected turbocode is
Unlike turbocodes and homogeneous TCCs, for such a ring-connected TCC the rate does not depend on the number of constituent trellises T. However, the permuter lengths are 1/T relative to the block length.
For ring-connected TCCs, we cannot use the same type of encoding procedure that is used for serially-concatenated convolutional codes, since some of the "output" bits must match some of the "input" bits. As with the homogeneous TCC, we can encode with approximately R(1 - R)N2 binary operations


Page 161
Figure 6.14
The performances of a homogeneous TCC and a
ring-connected TCC compared to the best rate 1/2
turbocode and low-density parity-check code
performances published to date.
by applying a previously computed generator matrix to a previously selected systematic subset of codeword bits. However, if R t = 3/4, t = 1, . . . , T (giving R = 1/2), encoding computations can be saved in the following way. First, we pick 3N/2T systematic bits for trellis 1 and generate the corresponding N/2T parity bits, using a number of computations that is linear in N. Then, we work around the ring and for each trellis pick N/2T systematic bits and generate the corresponding N/2T parity bits. When all but two trellises have been encoded, the last two trellises are used to form a binary vector that when multiplied by a previously computed N/T ´ N/T binary matrix yields the final set of N/T parity bits. The computations for the last step dominate and require approximately N 2/T 2 binary operations. For R avg = 3/4, R = 1/2, T = 3, the first method described above takes N2/4 operations, whereas the second method takes N2/9 operations.
As with homogeneous TCCs, when T > 2 ring-connected TCCs can be iteratively decoded using a variety of schedules. In my simulations, the trellises are processed sequentially while passing probabilities around the ring in one direction. Iterative decoding continues until a valid codeword is found or until 100 iterations are completed.


Page 162
The performance was estimated for an N = 131,070, T = 3 ring-connected TCC with n1 = n2 = n3 = 2/3 an R1 = R2 = R3 = 3/4, giving R = 1/2. Each rate 3/4 trellis used the generator polynomials (12, 4, 17, 11)octal, which were found by trial and error. The constituent codeword bits were shared alternately with the two neighboring trellises. Figure 6.14 shows the performance of this ring-connected TCC. It performs significantly better than the homogeneous TCC and the low-density parity-check code and only 0.2 dB worse than the turbocode.
6.5 Decoding complexity of iterative decoders
The decoding complexities per iteration for low-density parity-check codes, turbocodes, and homogeneous TCCs vary as significantly as do their proximities to Shannon's limit. The decoding complexity for a low-density parity-check code is roughly WGL = 6It multiplies per codeword bit, where I is the average number of iterations required to find the correct codeword, and t is the average number of checks with which each codeword bit participates [MacKay and Neal 1996].
For turbocodes and homogeneous TCCs, most of the computations are spent processing the constituent trellises. Each section of a bi-proper trellis requires roughly 6 ´ 2v multiplies to process, where 2v is the number of states in the regular trellis. For a turbocode with rate R and nt constituent convolutional codes, there are N Rnt trellis sections in all, so that the decoding complexity for a turbocode is roughly WTC = 6RInt2v multiplies per codeword bit. For a homogeneous TCC, there are Nnt trellis sections in all, so that the decoding complexity is roughly WHTCC = 6Int2v multiplies per codeword bit.
For example, at E b/N0 = 1.3 dB, the t » 3 low-density parity-check code discussed in the previous section has I = 11.2 (David MacKay, personal communication), so WGL = 202 multiplies per codeword bit. The R = 1/2, nt = 2, v = 4 turbocode has I = 10.5, so WTC = 509 multiplies per codeword bit. The R = 1/2, nt = 2, v = 4 homogeneous TCC has I = 10.5, so WHTCC = 2016 multiplies per codeword bit.
6.6 Parallel iterative decoding
The decoding algorithm for low-density parity-check codes proposed by Gallager [1963] and later by MacKay and Neal [1996] is inherently a parallel algorithm. As described in Section 6.2.5, probability propagation in the Bayesian network for a low-density parity-check code consists of passing sets of messages back and forth between the codeword bits and the clamped parity-check


Page 163
variables. It turns out that the standard decoders for turbocodes and serially-concatenated convolutional codes are inherently serial. In this section, I consider a parallel message-passing schedule. Not only do these results shed light on the insensitivity of turbodecoding to the message-passing schedule, but they also introduce the possibility of highly parallel turbodecoding hardware.
6.6.1 Concurrent turbodecoding
If each chain in a turbocode is viewed as a single vertex, then turbodecoding can also be viewed as a ''parallel'' algorithm.2 However, if each chain in a turbocode is viewed at a refined level (e.g., Figure 6.8), then the standard turbodecoding algorithm is inherently serial. That is, when messages are passed as shown in Figure 6.9, most of the computations are used to compute messages that cannot be propagated in parallel.
Here, I consider concurrent turbodecoding in which messages are passed in a parallel fashion. One time step of concurrent turbodecoding consists of simultaneously passing messages in both directions on all graph edges in the Bayesian network for the code. (Although "concurrent" is not quite the right term for such a parallel algorithm, the term "parallel" is used in the other name for turbocodes, "parallel concatenated convolutional codes".) Notice that concurrent turbodecoding is not just a parallel implementation of standard turbodecoding. It is a different algorithm which may have different properties.
A naive approach to a hardware implementation of concurrent turbodecoding would be to build one simple processor for each vertex in the Bayesian network for a code. Of course, for reasonably long block lengths, a prohibitively large number of these processors would be needed for a fully parallel VLSI implementation of concurrent turbodecoding. In the following section, I empirically compare the time complexity of standard decoding with concurrent decoding, while ignoring practical implementation issues such as wiring complexity. In practice, a more space-efficient implementation (e.g., time-shared processors) would be used at some detriment to the computational efficiency.
6.6.2 Results
The code used for the simulations was a rate 1/2 K = 5,000, N = 10,000 turbocode with two constituent convolutional codes, each with generator polynomials (21/37)octal. The constituent chains were connected by a randomly selected permuter. Every second output of each constituent chain was punctured to get a rate of 1/2. For each of three values of E b/N0, the transmission of
2 When there is more than one chain in a turbocode, messages may be passed between chains in either a serial or parallel manner.


Page 164
Figure 6.15
Performance of standard (S) and concurrent (C)
turbodecoding when implemented on a serial
computer, for 3 values of E b/N0.
Figure 6.16
Performance of standard (S) and concurrent (C)
turbodecoding when implemented on a
parallel computer, for 3 values of E b/N0.
107 information bits was simulated, and the results are shown in Figures 6.15 and 6.16. Interestingly, for a given E b/N0 it appears that both algorithms converge to the same BER.
Figure 6.15 shows the BER versus the number of messages passed in the constituent chains, for standard and concurrent turbodecoding. (The computation of the messages passed in the constituent chain dominates the decoding time). The number of messages passed gives a good indication of decoding complexity on a serial computer. Not surprisingly, it is apparent that the stan-


Page 165
dard algorithm is better suited to serial implementation. It is interesting that for a given BER, the concurrent decoding algorithm is roughly 30 times slower than the standard algorithm when implemented on a serial computer.
Figure 6.16 shows the BER versus the number of time steps for the case where 10,000 processors are available for concurrent turbodecoding. We assume that pipelining is not used for standard turbodecoding. For a given BER, the concurrent decoding algorithm is roughly 850 times faster than the standard algorithm when implemented on a parallel computer. If one processor is used for each half iteration of 5 iterations of standard turbodecoding in a pipeline fashion, standard decoding can be sped up by a factor of only 10 (extra pipeline stages do not improve the BER). Concurrent turbodecoding is still 85 times faster.
6.7 Speeding up iterative decoding by detecting variables early
The excellent bit error rate performance of iterative probability propagation decoders is achieved at the expense of a computationally burdensome decoding procedure. In this section, I present a method called early detection [Frey and Kschischang 1998] that can be used to reduce the computational complexity of a variety of iterative decoders. Using a confidence criterion, some information symbols, state variables and codeword symbols are detected early on in the iterative decoding procedure. In this way, the complexity of further processing is reduced with a controllable increase in BER. I present an easily implemented instance of this algorithm, called trellis splicing, that can be used with turbodecoding. For a simulated system of this type, I obtain a reduction in computational complexity of up to a factor of four, relative to a turbodecoder that performs the fewest iterations needed to achieve the same BER.
6.7.1 Early detection
It seems plausible that in some cases when iteratively decoding, some parts of the codeword may be more easily decoded than other parts. Although different parts of a codeword are usually inter-dependent, for particular noise patterns the coupling between certain parts may be weak. In these cases, it makes sense that the decoder should spend more computations on "tough" parts, and fewer computations on "easy" parts. During decoding, those parts that are deemed to be successfully decoded are clamped. Decoding computations are then focused on the remaining parts.
In general, the computation time of an iteration decreases with the number of variables detected early. So, in order to obtain the greatest speed-up, the decoder should early-detect as many variables as possible. However, an


Page 166
overly aggressive early-detection criterion will lead to a high rate of erroneous decisions, spoiling the BER performance. In addition to this constraint, the early-detection criterion should be relatively simple, so that the overhead of ascertaining which variables ought to be detected early does not overshadow the reduction in the computational complexity of subsequent iterative decoding. Here, I consider applying a threshold to the soft decisions for information bits, codeword bits, and state variables during decoding.
The soft decisions used for iterative decoding can be represented as log-odds ratios that approximate the true a posteriori log-odds ratios. The log-odds ratio for an information symbol, state variable, or codeword symbol z at iteration i given the channel output y is
where [image] is the approximation to the a posteriori distribution P(z|y) produced at iteration i. I will let i be fractional when the meaning is clear. For example, in a turbodecoder with two constituent codes, i = 0.5 refers to quantities produced by processing the first constituent code for the first time. By monitoring the log-odds ratios, we hope to accurately determine when to detect variables early.
6.7.2 Early detection for turbocodes: Trellis splicing
Here, I consider early detection of information symbols only. Consider the simple two-state trellis shown in Figure 6.17a. Let u k be the random variable for the information bit in the kth section of the trellis, and let s k be the random variable for the state at the beginning of the kth section of the trellis. The edge in the kth section of the trellis that leaves state s k Î {0,1} in response to information bit u k Î {0,1} has an associated branch metric, [image] . These metrics are determined from the received signals and the a priori probabilities regarding the transmitted information bit values. (In a systematic code, the likelihoods for the noisy received information bits can be included in the a priori probabilities.) If p(y k|u k, s k) is the likelihood function for the kth received signal and P(uk) is the a priori probability for information bit uk, then [image] . The forward pass consists of computing the flows from these metrics in the forward direction. This results in a flow value ak(sk) for each state s k at each section k, k = 0 . . . K - 1, computed as [image] and [image]. The backward pass simply consists of a flow computation in the reverse direction in order to obtain a flow value


Page 167
Figure 6.17
Trellis splicing. (a) shows a two-state trellis with
edges accompanied by information bit labels and
metrics and with nodes accompanied by flows.
(b) and (c): If we know that information bit k + 1
has a value of 1, we can cut the corresponding
section out of the trellis and splice the trellis back
together, introducing new information bit labels
and new metrics for the connecting edges.
for each state at each section: , and [image]. These two types of flow are combined to obtain the a posteriori log-odds ratio that each information bit is 1 versus 0, given the received signal sequence y:


Page 168
The computational cost of each section in the forward-backward algorithm thus consists of the time spent computing the a's and b's for each state, as well as the time spent computing the a posteriori log-odds ratios. Although there are various useful techniques and approximations for decreasing this cost [Hagenauer, Offer and Papke 1996; Benedetto et al. 1996], such as the SOVA [Hagenauer, Offer and Papke 1996], I will define it as a basic computational unit, and refer to it as a trellis section operation.
Suppose that according to some early detection criterion, we decide that the value of information bit u k+1 is 1. (Here, I will consider early detection for information bits only.) As a consequence, the trellis simplifies to the one shown in Figure 6.17b. The trellis can be simplified further by multiplying out the path metrics, giving the trellis shown in Figure 6.17c. Note that not only have the path metrics changed, but also the transitions now correspond to different information bit values. In general, portions of the trellis corresponding to early-detected information bits can be cut away, and the remaining segments spliced together with new path metrics and new information bit edge labels. If the values of b information bits are known, the spliced trellis will be b sections shorter, leading to a computational savings of b section operations for each future forward-backward sweep.
In order to implement trellis splicing, an integer array must be used to determine the state transitions, (s k, u k) ® s k+1. Whereas in the original trellis this mapping is very regular, after trellis splicing it is usually not. (E.g., the information bits associated with the outgoing edges of the kth state in Figure 6.17c have opposite values compared to those in Figure 6.17a.) The use of this array slightly increases the computational complexity of each section operation. Also, the array must be modified each time a section is cut away. However, both of these computational costs are insignificant compared to the cost of the basic section operation. In the implementation of trellis splicing used for the experiments presented in Section 6.7.3, I found that the percentage of cpu time spent on trellis splicing was less than 6%. The integer array also requires extra memory. However, the total memory used actually decreases during decoding while using trellis splicing. When a single section is cut away, the memory liberated by the elimination of gs, as and bs more than makes up for the extra integer array memory introduced. Moreover, if sections adjacent to the first are cut away, the transition array is simply modified, so that the memory associated with the gs, as and bs of the adjacent sections is completely recovered.
6.7.3 Experimental results
I simulated trellis splicing at E b/N0 = 0.1,0.2 and 0.3 dB, for a rate 1/3 unpunctured turbocode that had 10,000 information bits, identical constituent en-


Page 169
Figure 6.18
BER performance of turbodecoding with and without early detection for
thresholds increasing by 2, for (a) E b/N0 = 0.1 dB, (b) E b/N0 = 0.2 dB, and
(c) E b/N0 = 0.3 dB. (c) shows the computational gain for turbodecoding
with early detection compared to turbodecoding without early detection
as a function of BER, for E b/N0 = 0. 1, 0.2, and 0.3 dB.
coders (21,37)octal, and a randomly drawn permuter. After each half-iteration of turbodecoding, the log-odds ratio of each information bit was compared with a threshold in order to decide whether or not the bit should be detected early. For each value of the threshold, the number of bit errors and the number of section operations were averaged over the transmission of 20,000 information vectors. Figs. 6.18a, b, and c show plots of BER versus average number of section operations per information bit decoded, for a variety of thresholds. The curves for turbodecoding without early detection are also shown. In these simulations, a fixed number of decoding iterations were performed for each block.


Page 170
For a given BER, the computational complexity of decoding can be reduced the most compared to standard turbodecoding by using the threshold that corresponds to the curve in each figure that bottoms out at the prespecified BER. Thus, the locus of points corresponding to the heels of the curves gives the optimal achievable BER-complexity performance. Using these curves, we can answer the question, "At a specified E b/N0 and BER, what is the computational gain obtained by using early detection compared to using fewer decoding iterations without early detection?" The locus of points described above is interpolated in Figure 6.1 8d which shows the computational gain as a function of BER for the different values of E b/N0. For all three values of E b/N0, the greatest computational gain is obtained near the minimal BER.


Page 171
7
Future Research Directions
My goal in this book was to present to the reader a unified graphical model framework for describing problems and developing algorithms in the areas of pattern classification, unsupervised learning, data compression, and channel coding. The previous four chapters have shown how graphical models and various inference and learning algorithms can be applied to problems in these areas. In this chapter, I hope to give the reader a sense for some of what is to come in the near future of graphical models research in machine learning and digital communication.
7.1 Modularity and abstraction
I basically know of two principles for treating complicated systems in simple ways; the first is the principle of modularity and the second is the principle of abstraction. I am an apologist for computational probability in machine learning, and particularly for graphical models and variational methods, because I believe that probability theory implements these two principles in deep and intriguing ways—namely through factorization and through averaging. Exploiting these two mechanisms as fully as possible seems to me to be the way forward in machine learning.
Michael I. Jordan
Massachusetts Institute of Technology, 1997
Both in machine learning and digital communication, the goal is to perform inference in learned or designed models. Graphical models describe how a probabilistic description of a system can be factored into local functions. Exact inference consists of abstracting variables correctly according to probability theory, by marginalization. In intractable models, the approximate inference algorithms reviewed in Chapter 2 make abstraction tractable by approximately averaging over variables, whether by ignoring cycles in the graph, by using samples, by using approximating distributions, or by using inverse models.
Surely the next (current?) leap forward in machine learning will take advantage of the modular structure that we know is present in natural data. Data-driven algorithms (e.g., nearest neighbor methods) ignore this structure.
In digital communication, the presence of modularity is clear in such areas as data networks. The need for modularity in channel coding is less obvious,


Page 172
since the codes that Shannon used to prove that capacity is achievable were in fact completely random, albeit impractical to decode. However, the excellent performance of trellis-constrained codes (e.g., turbocodes) indicates that near Shannon limit performance can be obtained using highly-structured codes. In fact, this modularity is necessary for the approximate inference algorithms to work so well.
I think one important direction for further research is to develop inference algorithms that can figure out what should not be abstracted. For example, in variational inference, a parametric distribution is fit to the hidden variables. The form of this distribution influences the degree to which abstraction takes place. How do we choose this form? We certainly choose it so that it can be optimized easily. What other considerations are important? In fact, the most appropriate form of the variational distribution will even depend on the current input. As another example, loopy probability propagation ignores the cycles in the model graph. However, a small number of cycles can be processed quite easily, so why not take into account only those cycles that are most important for the current inference?
7.2 Faster inference and learning
The two fundamental problems are first to choose a probabilistic model capable of learning interesting things, and second to find ways of implementing inference in that model—approximately, of course. I'd like to comment on the latter problem. I think all our favourite implementation methods—Gibbs sampling, overrelaxation, hybrid Monte Carlo, variational methods, EM, gradient descent—are all too creepy-crawly slow. Of course, they're really snappy compared with the algorithms they replaced, but I think inference and learning should be even quicker. The world isn't an adversary. It should be possible to solve many learning problems in a couple of iterations through a reasonable data set, rather than thousands. It may be too much to ask for a one-shot learning method, but maybe we should be looking for one-and-a-half-shot learning algorithms.
David J.C. MacKay
Cambridge University, 1997
By using any fast inference algorithm while learning, the generalized expectation maximization algorithm will favor models for which the approximate inference algorithm works well. However, taken to the extreme, this approach


Page 173
will produce useless models in which we can make fast (and meaningless) inferences. It seems that there is a need for new cost functions that favor fast, but good probabilistic inferences.
Currently, the speed of inference is indirectly controlled by the number of iterations of loopy probability propagation, the Monte Carlo sample size, the simplicity of the variational distribution, and the simplicity of the recognition network. Is there a way to formulate a Lyapunov function that directly measures closeness to the correct answer and the number of iterations needed to get there?
With regards to learning, how does a machine recognize when it is finished learning? If the relative improvement in the cost function is measured, learning may go on indefinitely. Simply stopping learning after a specified time has elapsed does not guarantee the time was used efficiently. In Monte Carlo Bayesian learning, how does the sampler know when to stop? If the sampler is given 10 minutes, how should it proceed to obtain the ''best'' possible sample?
7.3 Scaling up to the brain
The brain is amazingly good at converting images into internal representations of their real causes. This could be the result of a huge collection of co-adapted clever tricks, but I think it is far more likely to be the result of a clever, general-purpose learning algorithm that can efficiently suck structure out of the rich sensory input. The search for this learning algorithm is one of the most exciting areas of science and is also of great technological importance. Graphical models are currently the most promising direction, but we do not yet know how to fit very large, hierarchical models to data in a computationally tractable way. It feels as if we are quite close to a really powerful learning algorithm, but something is still missing: None of the current algorithms scale well and as a result we can only tackle toy problems.
Geoffrey E. Hinton
University of Toronto, 1997
Why do current approximate inference and learning algorithms not scale well? It is certainly true that for richly-connected neural networks, the time needed for exact inference scales exponentially with the number of units. We would like our approximate algorithms to scale roughly linearly with the number of units. Is this too much to ask for? I don't believe so, for scale is a double-edged sword. The greater the number of sensory inputs, the less ambi-


Page 174
guity there tends to be in the inference. There is usually one explanation for a visual scene—multimodal inferences such as the orientation of the Necker cube are rare in everyday life. It seems to me that we need to find approximate algorithms that can allocate resources more efficiently in large networks.
Another important research topic is the interaction between the model and the approximate inference algorithm during learning, and how this interaction can help the scaling problem. Different types of model will be best-suited to different types of approximate inference algorithm. In unsupervised learning, the generalized expectation maximization algorithm favors models for which the approximate inference algorithm works well. Which properties of this interaction influence performance under scaling, and how can these properties be used to improve inference?
7.4 Improving model structures
For me, estimating joint distributions is a bit like playing God. You can't do everything!
Vladimir Vapnik
AT & T - Research, 1997
The premise of using estimates of the joint distribution to classify patterns is that by learning the natural generative process that produces each class of pattern, our machine can do a good job of classifying a test pattern by identifying the most appropriate generative process. In fact, if our joint distribution estimates are correct, this method is optimal. However, modelling the full joint distribution over the input patterns in this way is much more difficult than modelling the conditional probability of the class identity given the input pattern (as is done, e.g., in Le Cun et al. [1989] and Vapnik [1998]). In the former case, the machine tries to learn all aspects of the data, whereas in the latter case the machine tries to focus on only those aspects that influence the classification probability. While modelling conditional probabilities is more direct, it also ignores possible generalizations that can be made from apparently irrelevant input features to features that are relevant to the classification task.
A valid concern about attempts to model the full joint distribution is that the suboptimal solutions based on limited data sets may perform worse than simpler conditional probability estimates. As with much of neural network research, we find inspiration from cognitive science: Real neural networks do try to model everything in a mostly unsupervised fashion, and they succeed quite remarkably. Additional inspiration comes from the field of Bayesian learning, where proponents argue that given enough computer power, we ought


Page 175
to use models that are as complex as we believe the data to be, and then average over model predictions:
For problems where we do not expect a simple solution, the proper Bayesian approach is therefore to use a model of a suitable type that is as complex as we can afford computationally, regardless of the size of the training set.
Radford M. Neal [ 1996]
Of course, "Bayesian Bayesian networks" will require significant computational resources and so we will need to find appropriate models. What are the right models?
To date, research on learning in intractable graphical models with hidden variables has focused mainly on "vanilla" layered architectures of units with identical activation functions (e.g., sigmoid units, noisy-or combination functions). Clearly, these models are not most appropriate for a wide variety of real-world learning problems. It appears the time is right for a "killer application" of a latent-variable graphical model that will undoubtably use a well-thought out model structure. For example, it is probably a good idea to build in whatever physics are known, and leave some underlying physical parameters to be estimated or averaged over.
7.5 Iterative decoding
To me, it is clear that turbocodes are the most significant development in channel coding since Shannon's original 1948 paper. Of course, they weren't invented in a vacuum and build on ideas that came before, but nevertheless they are a truly remarkable discovery. Berrou et al. are primarily superb experimentalists, and they left plenty of work for the theorists to do. Some of the theorists have tackled the job of explaining why the classical turbocode construction should yield near-Shannon limit performance with optimal decoding. Some have found ways of modifying the classical construction to obtain near-Shannon limit performance with manageable complexity in a variety of situations (e.g., serial concatenation, product codes, "self"-concatenation, etc.).
However, right from the beginning, a few of us were vexed by the experimentally verified but theoretically unexplained iterative turbodecoding algorithm, which is, in my opinion, the key innovation of the French group. Turbodecoding is a low-complexity


Page 176
close approximation to the exponentially hard problem of exact decoding of turbocodes. While no one yet fully understands just why turbodecoding works as well as it does, we now know that there is a close connection between turbodecoding and many other algorithms known in the cryptography (Baum-Welch), information theory (Gallager-Tanner-MacKay-Wiberg), signal processing (backward-forward, FFT), statistics (expectation maximization), and artificial intelligence (probability propagation, Shafer-Shenoy theory) communities. All of these algorithms share the common feature of being distributed message passing schemes on graphs of various sorts. (I am at work on a tutorial paper with S. Aji as coauthor on the "generalized distributive law," which explains these connections.) However, unless the graph is a tree, these algorithms will not be exact, but (at best) are good approximations. Thus, turbodecoding and the Gallager-Tanner-MacKay-Wiberg algorithm are not fully understood, since the underlying graphs have cycles.
My own work is focused mainly on trying to provide a theoretically sound explanation for why these algorithms "work" when there are cycles. There have been small successes. The work of my group (Aji, Horn, and me) and that done independently by Weiss has pretty much cleared up the "one cycle" case. This is interesting but not directly relevant to turbodecoding since for classical turbocodes, the underlying graph is a bipartite graph with many cycles. (On the other hand it is directly relevant to the decoding of tail-biting codes.)
However, I believe that the success of turbodecoding will ultimately be found to be a special case of a theorem that says that the generalized distributive law, when applied iteratively to bipartite graphs with many cycles (and probably a broader class of graphs with cycles) using the sum-product semiring, gives, with high probability, decisions (though not conditional probabilities) on the hidden variables that are the same as those given by the brute force approach. I predict these theorems will be discovered in the next few years. (Of course, I hope we at Caltech will be the discoverers, but I am less confident of that!)
Robert J. McEliece
California Institute of Technology, 1997


Page 177
7.6 Iterative decoding in the real world
Real-world channels are not Gaussian. They have nonlinearities, bursty noise, multiple users, multiple reflections, delay constraints, and so on. An important direction for future research in digital communication is to extend the success of trellis-constrained codes such as turbocodes and low-density parity-check codes to real-world channels. One appealing approach is to build a graphical model that takes into account some of the effects mentioned above. Nonlinear channels can be modelled using (possibly parameterized) nonlinear conditional probabilities in the graphical model; bursty noise can be taken into account by a two-state ("bursty" versus "nonbursty") Markov model for the channel; each user in a multiple user system can attach extra graphical models that account for the other users; and so on. Optimal decoding is then a matter of inference in the total graphical model, and suboptimal decoding is accomplished using an approximate inference technique.
7.7 Unification
I think a great deal of perspective can be obtained by describing known algorithms in terms of inference and learning in graphical models. For example, we already know of many algorithms that are essentially probability propagation in appropriate graphical models: the forward-backward algorithm, Kalman filtering, extended Kalman filtering, the fast-Fourier transform, Gallager's iterative decoder and turbodecoding. What other connections can be made between known algorithms and the other inference techniques for graphical models? As these connections become clear, we can explore the use of different inference algorithms. How well does variational inference work for decoding turbocodes? Can we use Helmholtz machine recognition networks to assist decoding? How well does iterative probability propagation work for learning intractable graphical models? What new inference algorithms can be applied to new problems that can be described in terms of graphical models? I believe that in answering these questions and similar ones, we are sure to discover new fundamental connections and new powerful algorithms for inference and learning.


Page 179
References
Aji, S. M. and McEliece, R. J. (1997). A general algorithm for distributing information in a graph. In Proceedings of IEEE International Symposium on Information Theory.
Amari, S.-I. (1985). Differential-Geometrical Methods in Statistics, Lecture Notes in Statistics vol. 28. Springer, New York NY.
Amari, S.-I., Cichocki, A., and Yang, H. (1996). A new learning algorithm for blind signal separation. In Touretzky, D. S., Mozer, M. C., and Hasselmo, M. E., editors, Advances in Neural Information Processing Systems 8. MIT Press, Cambridge MA.
Bahl, L. R., Cocke, J., Jelinek, F., and Raviv, J. (1974). Optimal decoding of linear codes for minimizing symbol error rate. IEEE Transactions on Information Theory, 20:284-287.
Baum, L. E. and Petrie, T. (1966). Statistical inference for probabilistic functions of finite state markov chains. Annals of Mathematical Statistics, 37:1559-1563.
Becker, S. and Hinton, G. E. (1992). A self-organizing neural network that discovers surfaces in random-dot stereograms. Nature, 355:161-163.
Bell, A. J. and Sejnowski, T. J. (1995). An information maximization approach to blind separation and blind deconvolution. Neural Computation, 7:1129-1159.
Benedetto, S., Divsalar, D., Montorsi, G., and Pollara, F. (1996). Soft-ouput decoding algorithms in iterative decoding of parallel concatenated convolutional codes. Submitted to IEEE International Conference on Communications.
Benedetto, S. and Montorsi, G. (1996a). Iterative decoding of serially concatenated convolutional codes. Electronics Letters, 32:1186-1188.
Benedetto, S. and Montorsi, G. (1996b). Serial concatenation of block and convolutional codes. Electronics Letters, 32:887-888.


Page 180
Benedetto, S., Montorsi, G., Divsalar, D., and Pollara, F. (1997). Serial concatenation of interleaved codes: Performance analysis, design, and iterative decoding. To appear in IEEE Transactions on Information Theory.
Berrou, C. and Glavieux, A. (1996). Near optimum error correcting coding and decoding: Turbo-codes. IEEE Transactions on Communications, 44:1261-1271.
Berrou, C., Glavieux, A., and Thitimajshima, P. (1993). Near Shannon limit error-correcting coding and decoding: Turbo codes. In Proceedings of the IEEE International Conference on Communications.
Bishop, C. M., Svensén, M., and Williams, C. K. I. (1997). Gtm: the generative topographic mapping. To appear in Neural Computation.
Blahut, R. E. (1990). Digital Transmission of Information. Addison-Wesley Pub. Co., Reading MA.
Breiman, L., Friedman, J. H., Olshen, R. A., and Stone, C. J. (1984). Classification and regression trees. Wadsworth, Blemont CA.
Calderbank, A. R. and Sloane, N. J. A. (1987). New trellis codes based on lattices and cosets. IEEE Transactions on Information Theory, 33:177.
Chandler, D. (1987). Introduction to Modern Statistical Mechanics. Oxford University Press, New York NY. Chow, C. K. (1957). An optimum character recognition system using decision functions. IRE Transactions on Electronic Computing, 6:247-254.
Comon, P., Jutten, C., and Herault, J. (1991). Blind separation of sources. Signal Processing, 24:11-20.
Cooper, G. F. (1990). The computational complexity of probabilistic inference using Bayesian belief networks. Artificial Intelligence, 42:393-405.
Cover, T. M. and Thomas, J. A. (1991). Elements of Information Theory. John Wiley & Sons, New York NY.
Dagum, P. (1993). Approximating probabilistic inference in Bayesian belief networks is NP-hard. Artificial Intelligence, 60:141-153.
Dagum, P. and Chavez, R. M. (1993). Approximating probabilistic inference in Bayesian belief networks. IEEE Transactions on Pattern Analysis and Machine Intelligence, 15(3):246-255.


Page 181
Daut, D. G., Modestino, J. W., and Wismer, L. D. (1982). New short constraint length convolutional code constructions for selected rational rates. IEEE Transactions on Information Theory, 28(5):794-800.
Dayan, P. and Hinton, G. E. (1996). Varieties of Helmholtz machine. Neural Networks, 9:1385-1403.
Dayan, P., Hinton, G. E., Neal, R. M., and Zemel, R. S. (1995). The Helmholtz machine. Neural Computation, 7:889-904.
Dayan, P. and Zemel, R. S. (1995). Competition and multiple cause models. Neural Computation, 7:565-579.
Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Proceedings of the Royal Statistical Society, B-39: 1-38.
Devroye, L. (1986). Nonuniform Random Variate Generation. SpringerVerlag, New York NY.
Divsalar, D. and Pollara, F. (1995). Turbo-codes for PCS applications. In Proceedings of the International Conference on Communications, pages 54-59.
Duane, S., Kennedy, A. D., Pendleton, B. J., and Roweth, D. (1987). Hybrid Monte Carlo. Physical Letters B, 195:216-222.
Duda, R. O. and Hart, P. E. (1973). Pattern Classification and Scene Analysis. John Wiley, New York NY.
Fletcher, R. (1987). Practical methods of optimization. John Wiley & Sons, New York NY.
Foldiak, P. (1990). Forming sparse representations by local anti-hebbian learning. Biological Cybernetics, 64:165-170.
Forney, Jr., G. D. (1973). The Viterbi algorithm. Proceedings of the IEEE, 61 (3):268-277.
Forney, Jr., G. D. (1988). Coset codes - Part I: Introduction and geometrical classification. IEEE Transactions on Information Theory, 34:1123.
Forney, Jr., G. D. (1997). Approaching AWGN channel capacity with coset codes and multilevel coset codes. Submitted to IEEE Transactions on Information Theory.


Page 182
Frey, B. J. (1997a). Continuous sigmoidal belief networks trained using slice sampling. In Mozer, M. C., Jordan, M. I., and Petsche, T., editors, Advances in Neural Information Processing Systems 9. MIT Press, Cambridge MA. Available at http://www.cs.utoronto.ca/~frey.
Frey, B. J. (1997b). Variational inference for continuous sigmoidal Bayesian networks. In Sixth International Workshop on Artificial Intelligence and Statistics.
Frey, B. J. and Hinton, G. E. (1996). Free energy coding. In Storer, J. A. and Cohn, M., editors, Proceedings of the Data Compression Conference 1996. IEEE Computer Society Press. Available at http://www.cs.utoronto.ca/~frey.
Frey, B. J. and Hinton, G. E. (1997). Efficient stochastic source coding and an application to a Bayesian network source model. The Computer Journal, 40(2/3).
Frey, B. J. and Hinton, G. E. (1998). Variational learning in non-linear Gaussian belief networks. To appear in Neural Computation; available at http://www.cs.utoronto.ca/~frey.
Frey, B. J. and Kschischang, F. R. (1997). Probability propagation and iterative decoding. In Proceedings of the 34th Allerton Conference on Communication, Control and Computing 1996. Available at http://www.cs.utoronto.ca/~frey.
Frey, B. J. and Kschischang, F. R. (1998). Early detection and trellis splicing: Reduced-complexity iterative decoding. IEEE Journal on Selected Areas in Communications, 16:153-159.
Frey, B. J., Kschischang, F. R., Loeliger, H. A., and Wiberg, N. (1998). Factor graphs and algorithms. In Proceedings of the 35th Allerton Conference on Communication, Control and Computing 1997.
Frey, B. J. and MacKay, D. J. C. (1998a). A revolution: Belief propagation in graphs with cycles. In Jordan, M. I., Kearns, M. I., and Solla, S. A., editors, Advances in Neural Information Processing Systems 10. MIT Press. Available at http://www.cs.utoronto.ca/~frey.
Frey, B. J. and MacKay, D. J. C. (1998b). Trellis-constrained codes. In Proceedings of the 35th Allerton Conference on Communication, Control and Computing 1997.


Page 183
Gallager, R. G. (1963). Low-Density Parity-Check Codes. MIT Press, Cambridge MA.
Geman, S. and Geman, D. (1984). Stochastic relaxation, Gibbs distribution and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 6:721-741.
Ghahramani, Z. and Jordan, M. I. (1998). Factorial hidden Markov models. Machine Learning, 29:245-273.
Gilks, W. R. and Wild, P. (1992). Adaptive rejection sampling for Gibbs sampling. Applied Statistics, 41:337-348.
Goodman, L. A. (1970). The multivariate analysis of qualitative data: Interaction among multiple classifications. Journal of the American Statistical Association, 65:226-256.
Hagenauer, J., Offer, E., and Papke, L. (1996). Iterative decoding of binary block and convolutional codes. IEEE Transactions on Information Theory, 42(2):429-445.
Hammersley, J. M. and Handscomb, D. C. (1964). Monte Carlo Methods. Chapman and Hall, London England.
Heckerman, D. and Geiger, D. (1995). Learning Bayesian networks: a unification for discrete and Gaussian domains. In Besnard, P. and Hanks, S., editors, Proceedings of the Eleventh Conference on Uncertainty in Artificial Intelligence. Morgan Kaufmann.
Hinton, G. E., Dayan, P., Frey, B. J., and Neal, R. M. (1995). The wake-sleep algorithm for unsupervised neural networks. Science, 268:1158-1161.
Hinton, G. E. and Ghahramani, Z. (1997). Generative models for discovering sparse distributed representations. Philosophical Transactions of the Royal Society of London B, 352:1177-1190.
Hinton, G. E. and Sejnowski, T. J. (1986). Learning and relearning in Boltzmann machines. In Rumelhart, D. E. and McClelland, J. L., editors, Parallel Distributed Processing: Explorations in the Microstructure of Cognition, volume I, pages 282-317. MIT Press, Cambridge MA.
Hinton, G. E. and Zemel, R. S. (1994). Autoencoders, minimum description length and Helmholtz free energy. In Cowan, J. D., Tesauro, G., and Alspector, J., editors, Advances in Neural Information Processing Systems 6. Morgan Kaufmann.


Page 184
Hofmann, R. and Tresp, V. (1996). Discovering structure in continuous variables using Bayesian networks. In Touretzky, D., Mozer, M., and Hasselmo, M., editors, Advances in Neural Information Processing Systems 8. MIT Press.
Huffman, D. A. (1952). A method for the construction of minimum redundancy codes. Proceedings of the Institute of Radio Engineers, 40:1098-1101.
Imai, H. and Hirakawa, S. (1977). A new multilevel coding method using error-correcting codes. IEEE Transactions on Information Theory, 23:371-377. Correction, Nov. 1977, p. 784.
Jaakkola, T. and Jordan, M. I. (1997). A variational approach to Bayesian logistic regression models and their extensions. In Sixth International Workshop on Artificial Intelligence and Statistics.
Jaakkola, T., Saul, L. K., and Jordan, M. I. (1996). Fast learning by bounding likelihoods in sigmoid type belief networks. In Touretzky, D. S., Mozer, M. C., and Hasselmo, M. E., editors, Advances in Neural Information Processing Systems 8. MIT Press.
Jordan, M. I. (1995). Why the logistic function? A tutorial discussion on probabilities and neural networks. Technical Report Computational Cognitive Science 9503, MIT, Cambridge MA.
Kalos, M. H. and Whitlock, P. A. (1986). Monte Carlo Methods, Volume I: Basics. John Wiley, New York NY.
Kinderman, R. and Snell, J. L. (1980). Markov Random Fields and Their Applications. American Mathematical Society, Providence USA.
Kschischang, F. R. and Frey, B. J. (1998). Iterative decoding of compound codes by probability propagation in graphical models. IEEE Journal on Selected Areas in Communications, 16:219-230. Available at http://www.cs.utoronto.ca/~frey.
Lauritzen, S. L. (1996). Graphical Models. Oxford University Press, New York NY.
Lauritzen, S. L., Dawid, A. P., Larsen, B. N., and Leimer, H. G. (1990). Independence properties of directed Markov fields. Networks, 20:491-505.
Lauritzen, S. L. and Spiegelhalter, D. J. (1988). Local computations with probabilities on graphical structures and their application to expert systems. Journal of the Royal Statistical Society B, 50:157-224.


Page 185
Lauritzen, S. L. and Wermuth, N. (1989). Graphical models for associations between variables, some of which are qualitative and some quantitative. Annals of Statistics, 17:31-57.
Le Cun, Y., Boser, B., Denker, J. S., Henderson, D., Howard, R. E., Hubbard, W., and Jackel, L. D. (1989). Back-propagation applied to handwritten zip code recognition. Neural Computation, 1:541-551.
Lee, E. A. and Messerschmitt, D. G. (1994). Digital Communication. Kluwer Academic Publishers, Norwell MA.
Lin, S. and Costello, Jr., D. J. (1983). Error Control Coding: Fundamentals and Applications. Prentice-Hall Inc., Englewood Cliffs NJ.
Lodge, J., Young, R., Hoeher, P., and Hagenauer, J. (1993). Separable MAP 'filters' for the decoding of product and concatenated codes. In Proceedings of IEEE International Conference on Communications, pages 1740-1745.
MacKay, D. J. C. (1995). Bayesian neural networks and density networks. Nuclear Instruments and Methods in Physics Research, 354:73-80.
MacKay, D. J. C. (1997a). Good codes based on very sparse matrices. Submitted to IEEE Transactions on Information Theory.
MacKay, D. J. C. (1997b). Maximum likelihood and covariant algorithms for independent component analysis. Unpublished manuscript available at http://wol.ra.phy.cam.ac.uk/mackay.
MacKay, D. J. C. (1998). Information Theory, Inference and Learning Algorithms. Book in preparation, currently available at http://wol.ra.phy.cam.ac.uk/mackay.
MacKay, D. J. C. and Neal, R. M. (1995). Good codes based on very sparse matrices. In Boyd, C., editor, Cryptography and Coding. 5th IMA Conference, number 1025 in Lecture Notes in Computer Science, pages 100-111. Springer, Berlin Germany.
MacKay, D. J. C. and Neal, R. M. (1996). Near Shannon limit performance of low density parity check codes. Electronics Letters, 32(18): 1645-1646. Reprinted in Electronics Letters, vol. 33, March 1997, 457-458.
McCullagh, P. and Nelder, J. A. (1983). Generalized Linear Models. Chapman and Hall, London England.


Page 186
McEliece, R. J., MacKay, D. J. C., and Cheng, J. F. (1998). Turbo-decoding as an instance of Pearl's 'belief propagation' algorithm. IEEE Journal on Selected Areas in Communications, 16.
Meng, X. L. and Rubin, D. B. (1992). Recent extensions of the EM algorithm (with discussion). In Bernardo, J. M., Berger, J. O., Dawid, A. P., and Smith, A. F. M., editors, Bayesian Statistics 4. Clarendon Press, Oxford England.
Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., and Teller, E. (1953). Equation of state calculation by fast computing machines. Journal of Chemical Physics, 21:1087-1092.
Movellan, J. R. and McClelland, J. L. (1992). Learning continuous probability distributions with symmetric diffusion networks. Cognitive Science, 17:463-496.
Neal, R. M. (1992). Connectionist learning of belief networks. Artificial Intelligence, 56:71-113.
Neal, R. M. (1993). Probabilistic inference using Markov chain Monte Carlo methods. Unpublished manuscript available over the internet by ftp at ftp://ftp.cs.utoronto.ca/pub/radford/review.ps.Z.
Neal, R. M. (1996). Bayesian Learning for Neural Networks. Springer-Verlag, New York NY.
Neal, R. M. (1997). Markov chain Monte Carlo methods based on ''slicing'' the density function. University of Toronto, Dept. of Statistics Technical Report No. 9722.
Neal, R. M. and Hinton, G. E. (1993). A new view of the EM algorithm that justifies incremental and other variants. Unpublished manuscript available over the internet by ftp at ftp://ftp.cs.utoronto.ca/pub/radford/em.ps.Z
Pearl, J. (1986). Fusion, propagation, and structuring in belief networks. Artificial Intelligence, 29:241-288.
Pearl, J. (1987). Evidential reasoning using stochastic simulation of causal models. Artificial Intelligence, 32:245-257.
Pearl, J. (1988). Probabilistic Reasoning in Intelligent Systems. Morgan Kaufmann, San Mateo CA.


Page 187
Peterson, C. and Anderson, J. R. (1987). A mean field theory learning algorithm for neural networks. Complex Systems, 1:995-1019.
Potamianos, G. G. and Goutsias, J. K. (1993). Partition function estimation of Gibbs random field images using Monte Carlo simulations. IEEE Transactions on Information Theory, 39:1322-1332.
Rabiner, L. (1989). A tutorial on hidden Markov models and selected applications in speech recognition. Proceedings of the IEEE, 77:257-286.
Rasmussen, C. E. (1996). Evaluation of Gaussian Processes and Other Methods for Non-Linear Regression. Department of Computer Science, University of Toronto, Toronto Canada. Doctoral dissertation (ftp://ftp.cs.toronto.edu/pub/carl/thesis.ps.gz).
Rasmussen, C. E., Neal, R. M., Hinton, G. E., van Camp, D., Revow, M., Ghahramani, Z., Kustra, R., and Tibshirani, R. (1996). The DELVE Manual. University of Toronto, Toronto Canada. http://www.cs.utoronto.ca/~delve.
Ripley, B. D. (1987). Stochastic Simulation. John Wiley, New York NY.
Rissanen, J. (1989). Stochastic Complexity in Statistical Inquiry. World Scientific, Singapore.
Rissanen, J. and Langdon, G. G. (1976). Arithmetic coding. IBM Journal of Research and Development, 23:149-162.
Saul, L. K., Jaakkola, T., and Jordan, M. I. (1996). Mean field theory for sigmoid belief networks. Journal of Artificial Intelligence Research, 4:61-76.
Saund, E. (1995). A multiple cause mixture model for unsupervised learning. Neural Computation, 7:51-71.
Schubert, L. K. (1976). Extending the expressive power of semantic networks. Artificial Intelligence, 7:163-198.
Shannon, C. E. (1948). A mathematical theory of communication. Bell System Technical Journal, 27:379-423, 623-656.
Sheykhet, I. I. and Simkin, B. Y. (1990). Monte Carlo method in the theory of solutions. Computer Physics Reports, 12:67-133.


Page 188
Smyth, P., Heckerman, D., and Jordan, M. I. (1997). Probabilistic independence networks for hidden markov probability models. Neural Computation, 9(2):227-270.
Spiegelhalter, D. J. (1986). Probabilistic reasoning in predictive expert systems. In Kanal, L. N. and Lemmer, J. F., editors, Uncertainty in Artificial Intelligence, pages 47-68. North Holland, Amsterdam.
Spiegelhalter, D. J. (1990). Fast algorithms for probabilistic reasoning in influence diagrams, with applications in genetics and expert systems. In Oliver, R. M. and Smith, J. Q., editors, Influence Diagrams, Belief Nets, and Decision Analysis, pages 361-384. John Wiley & Sons, New York NY.
Spiegelhalter, D. J. and Lauritzen, S. L. (1990). Sequential updating of conditional probabilities on directed graphical structures. Networks, 20:579-605.
Tanner, R. M. (1981). A recursive approach to low complexity codes. IEEE Transactions on Information Theory, 27:533-547.
Tibshirani, R. (1992). Principal curves revisited. Statistics and Computing, 2:183-190.
Ungerboeck, G. (1982). Channel coding with multilevel/phase signals. IEEE Transactions on Information Theory, 28(1).
Vapnik, V. (1998). Statistical Learning Theory. John Wiley, New York NY.
Viterbi, A. J. and Omura, J. K. (1979). Principles of Digital Communication and Coding. McGraw-Hill, New York NY.
Vorobev, N. N. (1962). Consistent families of measures and their extensions. Theory of Probability and Applications, 7:147-163.
Wachsmann, U. and Huber, J. (1995). Power and bandwidth efficient digital communication using turbo-codes in multilevel codes. European Transactions on Telecommunications, 6(5):557-567.
Wallace, C. S. (1990). Classification by minimum-message-length inference. In S. G. Akl, et. al., editor, Advances in Computing and Information—ICCI 1990, number 468 in Lecture Notes in Computer Science. Springer, Berlin Germany.


Page 189
Wiberg, N. (1996). Codes and Decoding on General Graphs. Department of Electrical Engineering, Linköping University, Linköping Sweden. Doctoral dissertation.
Wiberg, N., Loeliger, H.-A., and Kötter, R. (1995). Codes and iterative decoding on general graphs. European Transactions on Telecommunications, 6:513-525.
Wicker, S. (1995). Error Control Systems for Digital Communications and Storage. Prentice-Hall Inc., Englewood Cliffs NJ.
Witten, I. H., Neal, R. M., and Cleary, J. G. (1987). Arithmetic coding for data compression. Communications of the ACM, 30:520-540.
Woods, W. A. (1975). What's in a link? Foundations for semantic networks. In Bobrow, D. and Collins, A., editors, Representation and understanding, pages 35-72. Academic Press, New York NY.
Wright, S. (1921). Correlation and causation. Journal of Agricultural Research, 20:557-585.
Zemel, R. E. (1993). A minimum description length framework for unsupervised learning. Department of Computer Science, University of Toronto, Toronto Canada. Doctoral dissertation.
Zhang, J. (1993). The mean field theory in EM procedures for blind Markov random field image restoration. IEEE Transactions on Image Processing, 2:27-40.


Page 191
Index
Symbols
B Q||P, 67
C, 130
D Q||p, 45
E b, 136
E b/N 0, 136
F Q||P, 45
H(.), 46
N0, 131
P, 130
R, 135
d(.,.), 18
Å, 135
dmin, 144
ak, 11
ck, 12
dk, 12
nk, 12
, 64
(.), 114
(.), 115
, 15
A
absorbing messages, 34
acyclic graph, 11
additive white Gaussian noise, 131
algebraic decoding, 139
ancestral ordering, 13, 110
ancestral simulation, 15, 38
approximate inference
Gibbs sampling, 39
Helmholtz machine, 49
Monte Carlo, 38
slice sampling, 42
variational inference, 43
arithmetic coding, 6, 109
using the ancestral ordering, 110
autoregressive network, 59, 60
ensemble prediction, 63
MAP estimation, 60
scaled priors, 62
AWGN, 131
B
band-limited channel, 129
bandwidth of channel, 129
bars problem, 23
image reconstruction, 93
noisy, 89
true difficulty, 95
Bayes' decision theory, 3
Bayes' rule, 3
Bayesian data compression with bits-back coding, 128
Bayesian network, 11
BER, 136
binary antipodal signalling, 135
binary code, 134
binary phase-shift keying, 135
binary sigmoidal Bayesian network, 22, 89
binary signalling, 135
bit error rate, 136
bits-back coding, 111
Bayesian data compression, 128
bits-back coding algorithm, 119
without auxiliary data, 121
blocked path, 16
buffered messages, 34
burglar alarm problem, 31
Helmholtz machine, 52
Monte Carlo inference, 31
variational inference, 46
C
capacity, 137


Page 192
capacity of an AWGN channel, 132
chain graph, 9
chain rule, 13
channel
band-limited, 129
bandwidth, 129
capacity, 132, 137
output waveform, 129
usage, 130
channel coding, 5, 7, 129
capacity, 132
decoder, 5, 130
encoder, 5, 130
Gaussian noise, 131
information, 130
linear binary code, 134
low-density parity-check code, 152
rate, 135
receiver, 129
Shannon's limit, 130, 149
signal constellation, 133
transmitter, 129
trellis-constrained code, 155
turbocode, 149
channel model, 7
children, 12
classification, 3
decision theory, 3
communication over a noisy channel, 129
conditional independence, 8, 13, 15
consistent method, 55
continuous sigmoidal Bayesian network, 97
continuous stereo disparity problem, 105
convolutional code, 17, 143
creating messages, 34
D
data compression, 5, 109
bits-back coding, 111
bits-back coding algorithm, 119
bits-back coding rate, 114
bits-back coding without auxiliary data, 121
decoder, 5
encoder, 5
decision theory, 3
decoding
algebraic, 139
convolutional code, 145
Hamming code, iteratively, 141
maximum likelihood, 140
probabilistic, 140
structure of iterative decoder, 154
symbol MAP, 140
symbol maximum a posteriori, 140
turbocode, iteratively, 149
delta rule, 71
DELVE results on classification, 85
dependency separation, 15
descendents, 12
digital communication, 5
disparity problem, 105
duplicating variables, 35
E
earthquake, 31
EM algorithm, 65
ensemble prediction
in latent variable model, 79
in the autoregressive network, 63
entropy, 46
erasure, 140
error-correcting code, 7
error-correcting coding, 129
estimation in latent variable models, 65
exact inference, 27, 37
expectation maximization algorithm, 65
generalized, 66
explaining away, 31


Page 193
F
factor graph, 9
fantasies in the Helmholtz machine, 75
flexible model, 57
free energy estimation, 71
function-to-variable message, 31
G
Gallager low-density parity-check code, 152
Gaussian noise, 131
generalization, 56
generalized EM
by Gibbs sampling, 70
by probability propagation, 70
by the wake-sleep algorithm, 74
by variational inference, 72
generalized expectation maximization algorithm, 66
generalized forward-backward algorithm, 29, 32
generator polynomials, 145
Gibbs sampling, 39, 105
graphical model, 8, 9
grouping variables, 35
H
Hamming code, 140
handwritten digits
classification, 80
examples, 55
Helmholtz machine, 49, 74
hidden variable model, 64, 68
ensemble prediction, 79
homogeneous trellis-constrained code, 159
I
i.i.d., 5
inference, 8, 27
information, 130
initializing a network, 34
iterative decoding complexity, 162
decoder structure, 154
early detection, 165
Hamming code, 141
in parallel, 162
turbocode, 149
L
latent variable model, 64, 68
ensemble prediction, 79
learning
by Gibbs sampling, 70
by probability propagation, 70
by the wake-sleep algorithm, 74
by variational inference, 72
learning in latent variable models, 65
likelihood estimation, 71
likelihood bound, 66
linear binary code, 134
linear code, 134
linear feedback shift register, 18, 145
local function, 10
log-likelihood estimation, 71
log-likelihood bound, 66
log-linear models, 1
logistic autoregressive network, 60
logistic function, 22
logistic regression, 3, 22, 60
low-density parity-check code, 152
M
MAP symbol decoding, 140
Markov chain, 15
Markov random field, 9, 10
maximum a posteriori symbol decoding, 140
maximum likelihood decoding, 140
message
absorption, 34
creation, 34
(continued on next page)


Page 194
message (continued from previous page)
function-to-variable message, 31
initialization, 34
passing, 34
propagation, 34
variable-to-function message, 31
message passing, 27, 34
messages
buffered, 34
schedule, 34
Monte Carlo, 15
Monte Carlo inference, 38
multiple-cause model, 68
multiple-model classification, 58
N
nearest neighbor classifier, 55
network initialization, 34
NLGBN, 102
variational inference, 103
noisy bars problem, 89
non-linear Gaussian Bayesian network, 102
rectified unit, 102
variational inference, 103
nondescendents, 12
NP-hard, 37
Nyquist sampling, 129
P
parallel concatenated convolutional code, 149
parameter estimation, 21
by Gibbs sampling, 70
by probability propagation, 70
by the wake-sleep algorithm, 74
by variational inference, 72
parameter estimation in latent variable models, 65
parameterized Bayesian network, 21
parents, 11
passing messages, 34
path blockage, 16
pattern classification, 3, 55
pattern classification results, 85
permutation map, 13
polynomials for a convolutional code, 145
potential, 10
probabilistic decoding, 140
probabilistic inference, 8, 27
probabilistic structure, 1, 8
probability propagation, 27, 32, 34
probability:structure, 1
product code, 152
propagating messages, 34
propagating probabilities, 34
R
rate of a code, 135
rate-distortion curve for binary-input channel, 138
receiver in channel coding, 129
recognition model, 49
rectified unit, 102
regularizer, 21
repetition code, 135
ring-connected trellis-constrained code, 160
S
scheduled messages, 34
serially-concatenated convolutional code, 152
Shannon's limit, 6, 7, 109, 130, 149
shaping, 133
shift problem, 105
sigmoidal Bayesian network, 22, 89
signal constellation, 133
signal-to-noise ratio, 136
signalling waveform, 129
simulation, 15
singly-connected graphical model, 27
sleep-phase learning, 75
slice sampling, 42, 99, 100
source coding, 5, 109
(continued on next page)


Page 195
source coding (continued from previous page)
bits-back coding, 111
bits-back coding algorithm, 119
bits-back coding rate, 114
bits-back coding without auxiliary data, 121
decoder, 5
encoder, 5
state variables, 7
stereo disparity problem, 105
sum-product algorithm, 28
supervised learning, 89
symbol error rate, 7
T
transmission power, 130
transmitter in channel coding, 129
trellis, 155
trellis-constrained code, 155
homogeneous, 159
low-density parity-check code, 156
ring-connected, 160
turbocode, 156
turbocode, 8, 20, 149
turbodecoding, 20
U
unobserved data, 64
unsupervised learning, 4, 89
V
variable-to-function message, 31
variational distribution, 46
variational inference, 22, 43, 103
variational learning, 103
W
wake-phase learning, 75
wake-sleep algorithm, 74
word error rate, 7


Page i

Graphical Models for Machine Learning and Digital Communication


Page ii
Adaptive Computation and Machine Learning
Thomas G. Dietterich, Editor
Christopher M. Bishop, David Heckerman, Michael I. Jordan, and Michael J. Kearns, Associate Editors
Bioinfomatics: The Machine Learning Approach, Pierre Baldi and Søren Brunak
Reinforcement Learning: An Introduction, Richard S. Sutton and Andrew G. Barto
Graphical Models for Machine Learning and Digital Communication, Brendan J. Frey


Page iii

Graphical Models for Machine Learning and Digital Communication
Brendan J. Frey
A Bradford Book
The MIT Press
Cambridge, Massachusetts
London, England


Page iv

Ó 1998 Massachusetts Institute of Technology
All rights reserved. No part of this book may be reproduced in any form by any electronic or mechanical means (including photocopying, recording, or information storage and retrieval) without permission in writing from the publisher.
This book was set in Times Roman by the author using TEX and was printed and bound in the United States of America.
Library of Congress Cataloging-in-Publication Data
Frey, Brendan J.
Graphical models for machine learning and digital communication/Brendan J. Frey.
p. cm.—(Adaptive computation and machine learning)
Includes bibliographical references and index.
ISBN 0-262-06202-X (alk. paper)
1. Machine learning. 2. Digital communications. 3. Graph theory. I. Title. II. Series.
Q325.5.F74 [image] 1998
006.3' 1—dc21 [image] [image] 98-13414
CIP


Page ix
Series Foreword
The goal of building systems that can adapt to their environments and learn from their experience has attracted researchers from many fields, including computer science, engineering, mathematics, physics, neuroscience, and cognitive science. Out of this research has come a wide variety of learning techniques that have the potential to transform many industrial and scientific fields. Recently, several research communities have begun to converge on a common set of issues surrounding supervised, unsupervised, and reinforcement learning problems. The MIT Press Series on Adaptive Computation and Machine Learning seeks to unify the many diverse strands of machine learning research and to foster high quality research and innovative applications.
This book by Brendan Frey is a perfect illustration of the convergence of these fields around a common set of mathematical and computational tools. Frey studies hierarchical probabilistic models, algorithms for learning and reasoning with them, and applications of these models to problems in supervised learning, unsupervised learning, data compression, and perhaps most surprisingly, error-control coding. The recent discovery that certain algorithms for decoding error-correcting codes are performing belief propagation over Bayesian belief networks is thrilling. Not only does it provide yet another connection between information theory and machine learning, but it also suggests new approaches to the design of improved error-correcting codes and to the design and analysis of learning algorithms. Frey provides the reader with a beautifully written introduction to graphical models and their associated algorithms. Then he takes us on a tour of the research frontier and of his recent work, which has advanced this frontier significantly. The book concludes with a thought-provoking assessment of future directions in this exciting research area.
THOMAS G. DIETTERICH


Page v
In memory of Granny Frey


Page vii
Contents
Series Foreword
ix
Preface
xi
1 Introduction
1
1.1 A probabilistic perspective
2
1.2 Graphical models: Factor graphs, Markov random fields and Bayesian belief networks
8
1.3 Organization of this book
25
2 Probabilistic Inference in Graphical Models
27
2.1 Exact inference using probability propagation (the sum-product algorithm)
27
2.2 Monte Carlo inference: Gibbs sampling and slice sampling
38
2.3 Variational inference
43
2.4 Helmholtz machines
49
3 Pattern Classification
55
3.1 Bayesian networks for pattern classification
58
3.2 Autoregressive networks
59
3.3 Estimating latent variable models using the EM algorithm
64
3.4 Multiple-cause networks
68
3.5 Classification of handwritten digits
80
4 Unsupervised Learning
89
4.1 Extracting structure from images using the wake-sleep algorithm
89
4.2 Simultaneous extraction of continuous and categorical structure
96
4.3 Nonlinear Gaussian Bayesian networks (NLGBNs)
102
5 Data Compression
109
5.1 Fast compression with Bayesian networks
110
5.2 Communicating extra information through the codeword choice
111
5.3 Relationship to maximum likelihood estimation
117
5.4 The "bits-back" coding algorithm
119


Page viii
5.5 Experimental results
123
5.6 Integrating over model parameters using bits-back coding
128
6 Channel Coding
129
6.1 Review: Simplifying the playing field
131
6.2 Graphical models for error correction: Turbocodes, low-density parity-check codes and more
139
6.3 "A code by any other network would not decode as sweetly"
154
6.4 Trellis-constrained codes (TCCs)
155
6.5 Decoding complexity of iterative decoders
162
6.6 Parallel iterative decoding
162
6.7 Speeding up iterative decoding by detecting variables early
165
7 Future Research Directions
171
7.1 Modularity and abstraction
171
7.2 Faster inference and learning
172
7.3 Scaling up to the brain
173
7.4 Improving model structures
174
7.5 Iterative decoding
175
7.6 Iterative decoding in the real world
177
7.7 Unification
177
References
179
Index
191


Page xi
Preface
What is this book about?
A variety of problems in machine learning and digital communication deal with complex but structured natural or artificial systems. Natural patterns that we wish to automatically classify are a consequence of a hierarchical causal physical process. Learning about the world in which we live requires that we extract useful sensory features and abstract concepts and then form a model for how these interact. Universal data compression involves estimating the probability distribution of a data source, which is often produced by some natural hierarchical process. Error-correcting codes used in telephone modems and deep-space communication consist of electrical signals that are linked together in a complex fashion determined by the designed code and the physical nature of the communication channel. Not only are these tasks characterized by complex structure, but they also contain random elements. Graphical models such as Bayesian belief networks and Markov random fields provide a way to describe the relationships between random variables in a complex stochastic system.
In this book, I use graphical models as an overarching framework to describe and solve problems in the areas of pattern classification, unsupervised learning, data compression, and channel coding. This book covers research I did while in my doctoral program at the University of Toronto. Rather than being a textbook, this book is a treatise that covers several leading-edge areas of research in machine learning and digital communication.
The book begins with a review of graphical models and algorithms for inferring probabilities in these models, including the probability propagation algorithm, Markov chain Monte Carlo (Gibbs sampling), variational inference, and Helmholtz machines. I then turn to the practical problem of learning models for pattern classification. Results on the classification of handwritten digits show that Bayesian network pattern classifiers outperform other standard methods, such as the k-nearest neighbor method.
In the area of unsupervised learning, I show that with just a simple local delta-rule learning algorithm (the "wake-sleep" algorithm) Bayesian network—neural network hybrids can learn hierarchical structure from images. When Bayesian networks with hidden variables are used as source models for data compression, an exponentially large number of codewords are asso-


Page xii
ciated with each input pattern. However, it turns out that the code can still be used efficiently, if a new technique called ''bits-back coding'' is used.
The current best error-correcting decoding algorithms (e.g., turbodecoding) are instances of "probability propagation" in various graphical models. These new schemes are rapidly closing the gap between the performances of practical channel coding systems and Shannon's fifty-year-old channel coding limit. The graphical model framework exposes the similarities between these codes and leads the way to a new class of "trellis-constrained codes" which also operate close to Shannon's limit.
The book concludes with suggestions by myself and other experts for future directions of research in graphical models for machine learning and digital communication.
Nomenclature
Vectors, matrices, high-dimensional matrices and sets of variables are written in boldface type. (Vectors are usually written in lower-case type.) Sets are quite different from vectors, but this abuse of notation permits set operations (e.g., "Í", "\") while at the same time permitting cardinal access to the set members (e.g., weighted sum of the elements via indexing). I use curly braces {. . .} to write the elements in a set or vector of variables. {z k} is the set containing a singleton variable zk whereas [image] = {z1, . . . , z k} Extra labels on variables usually appear as upper-case type in superscripts (e.g., q v whereas vector, matrix, and high-dimensional matrix indices usually appear as subscripts (e.g., [image]) For example, we can write the following with respect to the set of parameters [image]. Some types of index (notably training case indices) appear as superscripts in parentheses (e.g., v(t)).
Probability mass functions are usually written in upper-case italics (e.g. P(·), Q(·)) whereas probability density functions are usually written in lower case italics (e.g. p(·), q(·)) The distribution is identified by the random variable, so the distribution P(v) is different from P(x). Also, to keep the formulas short, the symbols for a random variable and its value are usually the same. So, P( u k |y) sometimes refers to the probability that the random variable U k takes on the value u k and at other times refers to the set of probabilities corresponding to the values that U k can take on. In cases where a random variable and its value must be distinguishable, I write an assignment. So, [image] means [image] A distribution subscripted with “r” refers to the correct, or "real" distribution. For example, if P(v) is a model distribution, we hope that P(v) » Pr(v). E[·] is an expectation with respect to the stated distribution.


Page xiii
I use log(·) for the natural logarithm and logx(·) for the logarithm to the base x.
Acknowledgments
The material in this book is a consequence of the interactions I have had with several excellent researchers. I thank my doctoral advisor Geoffrey Hinton for his guidance; I greatly value his open-mindedness, his inspirational discussions, and his honest criticisms. I also thank Radford Neal, whose creativity and persistent pursuit of precision has certainly enhanced my research and this book. I am very grateful to Frank Kschischang for valuable discussions on error-correcting coding and for introducing me to the information theory community. Glenn Gulak was very helpful in suggesting that there might be similarities between my research on machine learning and the turbodecoding algorithm for error-correcting codes. I also greatly appreciate recent energetic collaborations with David MacKay, who I find is lots of fun to work with.
Thanks also go to the following people for valuable conversations: Peter Dayan, Zoubin Ghahramani, Tommi Jaakkola, Michael Jordan, Carl Rasmussen, Virginia de Sa, Brian Sallans, Lawrence Saul, and Chris Williams.
My love goes to my wife Utpala for being supportive of my interest in research and to my children Shardul and Sarth for motivating me to get some work done. Last, but not least, I also thank my parents John and Cecelia for providing me with an environment conducive to curiosity and exploration, while I was growing up.
I was financially supported by a Natural Sciences and Engineering Research Council 1967 Science and Engineering Scholarship, a Walter C. Sumner Memorial Fellowship, and a University of Toronto Open Fellowship. My work was financially supported by grants from the Natural Sciences and Engineering Research Council, the Institute for Robotics and Intelligent Systems, and the Arnold and Mabel Beckman Foundation.