Factorial Hidden Markov Models.
Hidden Markov models (HMMs) have proven to be one of the most widely used tools for learning probabilistic models of time series data. In an HMM, information about the past is conveyed through a single discrete variablethe hidden state. We discuss a generalization of HMMs in which this state is factored into multiple state variables and is therefore represented in a distributed manner. We describe an exact algorithm for inferring the posterior probabilities of the hidden state variables given the observations, and relate it to the forwardbackward algorithm for HMMs and to algorithms for more general graphical models. Due to the combinatorial nature of the hidden state representation, this exact algorithm is intractable. As in other intractable systems, approximate inference can be carried out using Gibbs sampling or variational methods. Within the variational framework, we present a structured approximation in which the the state variables are decoupled, yielding a tractable algorithm for learning the parameters of the model. Empirical comparisons suggest that these approximations are efficient and provide accurate alternatives to the exact methods. Finally, we use the structured approximation to model Bachs chorales and show that factorial HMMs can capture statistical structure in this data set which an unconstrained HMM cannot.
Introduction
Due to its flexibility and to the simplicity and efficiency of its parameter estimation
algorithm, the hidden Markov model (HMM) has emerged as one of the basic statistical
tools for modeling discrete time series, finding widespread application in the
areas of speech recognition (Rabiner & Juang, 1986) and computational molecular
biology (Krogh, Brown, Mian, Sj-olander, & Haussler, 1994). An HMM is essentially
a mixture model, encoding information about the history of a time series in
the value of a single multinomial variable-the hidden state-which can take on
one of K discrete values. This multinomial assumption supports an efficient parameter
estimation algorithm-the Baum-Welch algorithm-which considers each
of the K settings of the hidden state at each time step. However, the multinomial
assumption also severely limits the representational capacity of HMMs. For exam-
ple, to represent bits of information about the history of a time sequence, an
HMM would need distinct states. On the other hand, an HMM with a
distributed state representation could achieve the same task with binary state
Z. GHAHRAMANI AND M.I. JORDAN
variables (Williams & Hinton, 1991). This paper addresses the problem of constructing
efficient learning algorithms for hidden Markov models with distributed
state representations.
The need for distributed state representations in HMMs can be motivated in two
ways. First, such representations let the model automatically decompose the state
space into features that decouple the dynamics of the process that generated the
data. Second, distributed state representations simplify the task of modeling time
series that are known a priori to be generated from an interaction of multiple,
loosely-coupled processes. For example, a speech signal generated by the superposition
of multiple simultaneous speakers can be potentially modeled with such an
architecture.
Williams and Hinton (1991) first formulated the problem of learning in HMMs
with distributed state representations and proposed a solution based on deterministic
learning. 1 The approach presented in this paper is similar to Williams
and Hinton's in that it can also be viewed from the framework of statistical mechanics
and mean field theory. However, our learning algorithm is quite different
in that it makes use of the special structure of HMMs with a distributed state
representation, resulting in a significantly more efficient learning procedure. Anticipating
the results in Section 3, this learning algorithm obviates the need for
the two-phase procedure of Boltzmann machines, has an exact M step, and makes
use of the forward-backward algorithm from classical HMMs as a subroutine. A
different approach comes from Saul and Jordan (1995), who derived a set of rules
for computing the gradients required for learning in HMMs with distributed state
spaces. However, their methods can only be applied to a limited class of architectures

Hidden Markov models with distributed state representations are a particular
class of probabilistic graphical model (Pearl, 1988; Lauritzen & Spiegelhalter, 1988),
which represent probability distributions as graphs in which the nodes correspond
to random variables and the links represent conditional independence relations.
The relation between hidden Markov models and graphical models has recently
been reviewed in Smyth, Heckerman and Jordan (1997). Although exact probability
propagation algorithms exist for general graphical models (Jensen, Lauritzen, &
Olesen, 1990), these algorithms are intractable for densely-connected models such
as the ones we consider in this paper. One approach to dealing with this issue is
to utilize stochastic sampling methods (Kanazawa et al., 1995). Another approach,
which provides the basis for algorithms described in the current paper, is to make
use of variational methods (cf. Saul, Jaakkola, & Jordan, 1996).
In the following section we define the probabilistic model for factorial HMMs
and in Section 3 we present algorithms for inference and learning. In Section 4 we
describe empirical results comparing exact and approximate algorithms for learning
on the basis of time complexity and model quality. We also apply factorial HMMs
to a real time series data set consisting of the melody lines from a collection of
chorales by J. S. Bach. We discuss several generalizations of the probabilistic model
FACTORIAL HIDDEN MARKOV MODELS 247
Y t+1
Y t-1
Y
S (2)
Y
S (2)
Y t+1
S (2)
Y t-1
(a) (b)

Figure

1. (a) A directed acyclic graph (DAG) specifying conditional independence relations for
a hidden Markov model. Each node is conditionally independent from its non-descendants given
its parents. (b) A DAG representing the conditional independence relations in a factorial HMM
with underlying Markov chains.
in Section 5, and we conclude in Section 6. Where necessary, details of derivations
are provided in the appendixes.
2. The probabilistic model
We begin by describing the hidden Markov model, in which a sequence of observations
modeled by specifying a probabilistic relation
between the observations and a sequence of hidden states fS t g, and a Markov
transition structure linking the hidden states. The model assumes two sets of conditional
independence relations: that Y t is independent of all other observations and
states given S t , and that S t is independent of
Markov property). Using these independence relations, the joint probability for the
sequence of states and observations can be factored as
Y
The conditional independencies specified by equation (1) can be expressed graphically
in the form of Figure 1 (a). The state is a single multinomial random variable
that can take one of K discrete values, S t Kg. The state transition
probabilities, are specified by a K \Theta K transition matrix. If the observations
are discrete symbols taking on one of D values, the observation probabilities
can be fully specified as a K \Theta D observation matrix. For a continuous
observation vector, P (Y t jS t ) can be modeled in many different forms, such as a
Gaussian, a mixture of Gaussians, or even a neural network. 2
In the present paper, we generalize the HMM state representation by letting the
state be represented by a collection of state variables
Z. GHAHRAMANI AND M.I. JORDAN
each of which can take on K (m) values. We refer to these models as factorial
hidden Markov models, as the state space consists of the cross product of these state
variables. For simplicity, we will assume that K although all the
results we present can be trivially generalized to the case of differing K (m) . Given
that the state space of this factorial HMM consists of all K M combinations of the
t variables, placing no constraints on the state transition structure would result
in a K M \Theta K M transition matrix. Such an unconstrained system is uninteresting
for several reasons: it is equivalent to an HMM with K M states; it is unlikely
to discover any interesting structure in the K state variables, as all variables are
allowed to interact arbitrarily; and both the time complexity and sample complexity
of the estimation algorithm are exponential in M .
We therefore focus on factorial HMMs in which the underlying state transitions
are constrained. A natural structure to consider is one in which each state variable
evolves according to its own dynamics, and is a priori uncoupled from the other
state variables:
Y
A graphical representation for this model is presented in Figure 1 (b). The transition
structure for this system can be represented as M distinct K \Theta K matrices.
Generalizations that allow coupling between the state variables are briefly discussed
in Section 5.
As shown in Figure 1 (b), in a factorial HMM the observation at time step t can
depend on all the state variables at that time step. For continuous observations,
one simple form for this dependence is linear Gaussian; that is, the observation Y t
is a Gaussian random vector whose mean is a linear function of the state variables.
We represent the state variables as K \Theta 1 vectors, where each of the K discrete
values corresponds to a 1 in one position and 0 elsewhere. The resulting probability
density for a D \Theta 1 observation vector Y t is
ae
oe
where
Each W (m) matrix is a D \Theta K matrix whose columns are the contributions to the
means for each of the settings of S (m)
, C is the D \Theta D covariance matrix, 0 denotes
matrix transpose, and j \Delta j is the matrix determinant operator.
One way to understand the observation model in equations (4a) and (4b) is to
consider the marginal distribution for Y t , obtained by summing over the possible
states. There are K settings for each of the M state variables, and thus there
FACTORIAL HIDDEN MARKOV MODELS 249
are K M possible mean vectors obtained by forming sums of M columns where one
column is chosen from each of the W (m) matrices. The resulting marginal density
of Y t is thus a Gaussian mixture model, with K M Gaussian mixture components
each having a constant covariance matrix C. This static mixture model, without
inclusion of the time index and the Markov dynamics, is a factorial parameterization
of the standard mixture of Gaussians model that has interest in its own right (Zemel,
1993; Hinton & Zemel, 1994; Ghahramani, 1995). The model we are considering in
the current paper extends this model by allowing Markov dynamics in the discrete
state variables underlying the mixture. Unless otherwise stated, we will assume the
Gaussian observation model throughout the paper.
The hidden state variables at one time step, although marginally independent,
become conditionally dependent given the observation sequence. This can be determined
by applying the semantics of directed graphs, in particular the d-separation
criterion (Pearl, 1988), to the graphical model in Figure 1 (b). Consider the Gaussian
model in equations (4a)-(4b). Given an observation vector Y t , the posterior
probability of each of the settings of the hidden state variables is proportional to the
probability of Y t under a Gaussian with mean - t . Since - t is a function of all the
state variables, the probability of a setting of one of the state variables will depend
on the setting of the other state variables. 3 This dependency effectively couples all
of the hidden state variables for the purposes of calculating posterior probabilities
and makes exact inference intractable for the factorial HMM.
3. Inference and learning
The inference problem in a probabilistic graphical model consists of computing
the probabilities of the hidden variables given the observations. In the context
of speech recognition, for example, the observations may be acoustic vectors and
the goal of inference may be to compute the probability for a particular word or
sequence of phonemes (the hidden state). This problem can be solved efficiently
via the forward-backward algorithm (Rabiner & Juang, 1986), which can be shown
to be a special case of the Jensen, Lauritzen, and Olesen (1990) algorithm for
probability propagation in more general graphical models (Smyth et al., 1997). In
some cases, rather than a probability distribution over hidden states it is desirable
to infer the single most probable hidden state sequence. This can be achieved via
the Viterbi (1967) algorithm, a form of dynamic programming that is very closely
related to the forward-backward algorithm and also has analogues in the graphical
model literature (Dawid, 1992).
The learning problem for probabilistic models consists of two components: learning
the structure of the model and learning its parameters. Structure learning is a
topic of current research in both the graphical model and machine learning communities
(e.g. Heckerman, 1995; Stolcke & Omohundro, 1993). In the current paper we
deal exclusively with the problem of learning the parameters for a given structure.
Z. GHAHRAMANI AND M.I. JORDAN
3.1. The EM algorithm
The parameters of a factorial HMM can be estimated via the expectation maximization
(EM) algorithm (Dempster, Laird, & Rubin, 1977), which in the case of
classical HMMs is known as the Baum-Welch algorithm (Baum, Petrie, Soules, &
Weiss, 1970). This procedure iterates between a step that fixes the current parameters
and computes posterior probabilities over the hidden states (the E step) and
a step that uses these probabilities to maximize the expected log likelihood of the
observations as a function of the parameters (the M step). Since the E step of EM
is exactly the inference problem as described above, we subsume the discussion of
both inference and learning problems into our description of the EM algorithm for
factorial HMMs.
The EM algorithm follows from the definition of the expected log likelihood of
the complete (observed and hidden) data:
Q(OE new
log
where Q is a function of the parameters OE new , given the current parameter estimate
OE and the observation sequence fY t g. For the factorial HMM the parameters
of the model are
consists of computing Q. By expanding (5)
using equations (1)-(4b), we find that Q can be expressed as a function of three
types of expectations over the hidden state variables: hS (m)
t i, and
t i, where h\Deltai has been used to abbreviate E f\DeltajOE; fY t gg. In the HMM
notation of Rabiner and Juang (1986), hS (m)
corresponds to fl t , the vector of
state occupation probabilities, hS (m)
corresponds to - t , the K \Theta K matrix of
state occupation probabilities at two consecutive time steps, and hS (m)
t i has
no analogue when there is only a single underlying Markov model. The M step uses
these expectations to maximize Q as a function of OE new . Using Jensen's inequality,
Baum, Petrie, Soules & Weiss (1970) showed that each iteration of the E and M
steps increases the likelihood, P (fY t gjOE), until convergence to a (local) optimum.
As in hidden Markov models, the exact M step for factorial HMMs is simple
and tractable. In particular, the M step for the parameters of the output model
described in equations (4a)-(4b) can be found by solving a weighted linear regression
problem. Similarly, the M steps for the priors, - (m) , and state transition matrices,
are identical to the ones used in the Baum-Welch algorithm. The details
of the M step are given in Appendix A. We now turn to the substantially more
difficult problem of computing the expectations required for the E step.
3.2. Exact inference
Unfortunately, the exact E step for factorial HMMs is computationally intractable.
This fact can best be shown by making reference to standard algorithms for prob-
FACTORIAL HIDDEN MARKOV MODELS 251
abilistic inference in graphical models (Lauritzen & Spiegelhalter, 1988), although
it can also be derived readily from direct application of Bayes rule. Consider the
computations that are required for calculating posterior probabilities for the factorial
HMM shown in Figure 1 (b) within the framework of the Lauritzen and
Spiegelhalter algorithm. Moralizing and triangulating the graphical structure for
the factorial HMM results in a junction tree (in fact a chain) with
cliques of size M+1. The resulting probability propagation algorithm has time complexity
O(TMK M+1 ) for a single observation sequence of length T . We present a
forward-backward type recursion that implements the exact E step in Appendix B.
The naive exact algorithm which consists of translating the factorial HMM into an
equivalent HMM with K M states and using the forward-backward algorithm, has
time complexity O(TK 2M ). Like other models with multiple densely-connected
hidden variables, this exponential time complexity makes exact learning and inference
intractable.
Thus, although the Markov property can be used to obtain forward-backward-
like factorizations of the expectations across time steps, the sum over all possible
configurations of the other hidden state variables within each time step is unavoid-
able. This intractability is due inherently to the cooperative nature of the model:
for the Gaussian output model, for example, the settings of all the state variables
at one time step cooperate in determining the mean of the observation vector.
3.3. Inference using Gibbs sampling
Rather than computing the exact posterior probabilities, one can approximate them
using a Monte Carlo sampling procedure, and thereby avoid the sum over exponentially
many state patterns at some cost in accuracy. Although there are many
possible sampling schemes (for a review see Neal, 1993), here we present one of the
simplest-Gibbs sampling (Geman & Geman, 1984). For a given observation sequence
fY t g, this procedure starts with a random setting of the hidden states fS t g.
At each step of the sampling process, each state vector is updated stochastically
according to its probability distribution conditioned on the setting of all the other
state vectors. The graphical model is again useful here, as each node is conditionally
independent of all other nodes given its Markov blanket, defined as the set of
children, parents, and parents of the children of a node. To sample from a typical
state variable S (m)
t we only need to examine the states of a few neighboring nodes:
t sampled from P (S (m)
Sampling once from each of the TM hidden variables in the model results in a
new sample of the hidden state of the model and requires O(TMK) operations.
The sequence of overall states resulting from each pass of Gibbs sampling defines
a Markov chain over the state space of the model. Assuming that all probabilities
are bounded away from zero, this Markov chain is guaranteed to converge to the
Z. GHAHRAMANI AND M.I. JORDAN
posterior probabilities of the states given the observations (Geman & Geman, 1984).
Thus, after some suitable time, samples from the Markov chain can be taken as
approximate samples from the posterior probabilities. The first and second-order
statistics needed to estimate hS (m)
are collected using
the states visited and the probabilities estimated during this sampling process are
used in the approximate E step of EM. 4
3.4. Completely factorized variational inference
There also exists a second approximation of the posterior probability of the hidden
states that is both tractable and deterministic. The basic idea is to approximate the
posterior distribution over the hidden variables P (fS t gjfY t g) by a tractable distribution
Q(fS t g). This approximation provides a lower bound on the log likelihood
that can be used to obtain an efficient learning algorithm.
The argument can be formalized following the reasoning of Saul, Jaakkola, and
Jordan (1996). Any distribution over the hidden variables Q(fS t g) can be used to
define a lower bound on the log likelihood
log
log
Q(fS t g) log
where we have made use of Jensen's inequality in the last step. The difference
between the left-hand side and the right-hand side of this inequality is given by the
Kullback-Leibler divergence (Cover & Thomas, 1991):
Q(fS t g) log
The complexity of exact inference in the approximation given by Q is determined
by its conditional independence relations, not by its parameters. Thus, we can chose
Q to have a tractable structure-a graphical representation that eliminates some
of the dependencies in P . Given this structure, we are free to vary the parameters
of Q so as to obtain the tightest possible bound by minimizing (6).
We will refer to the general strategy of using a parameterized approximating distribution
as a variational approximation and refer to the free parameters of the
distribution as variational parameters. To illustrate, consider the simplest variational
approximation, in which the state variables are assumed independent given
the observations This distribution can be written as
FACTORIAL HIDDEN MARKOV MODELS 253
S (2)
S (2)
S (2)
S (2)
S (2)
S (2)
(a) (b)

Figure

2. (a) The completely factorized variational approximation assumes that all the state variables
are independent (conditional on the observation sequence). (b) The structured variational
approximation assumes that the state variables retain their Markov structure within each chain,
but are independent across chains.
Y
Y
Q(S (m)
The variational parameters,
t g, are the means of the state variables, where,
as before, a state variable S (m)
t is represented as a K-dimensional vector with a 1
in the k th position and 0 elsewhere, if the m th Markov chain is in state k at time t.
The elements of the vector ' (m)
therefore define the state occupation probabilities
for the multinomial variable S (m)
t under the distribution Q:
Q(S (m)
Y
t;k
t;k
This completely factorized approximation is often used in statistical physics, where
it provides the basis for simple yet powerful mean field approximations to statistical
mechanical systems (Parisi, 1988).
To make the bound as tight as possible we vary ' separately for each observation
sequence so as to minimize the KL divergence. Taking the derivatives of (6) with
respect to ' (m)
t and setting them to zero, we obtain the set of fixed point equations
(see

Appendix

C) defined by
new
ae
Y (m)
oe
where ~
Y (m)
t is the residual error in Y t given the predictions from all the state
variables not including m:
~
Y (m)
'6=m
Z. GHAHRAMANI AND M.I. JORDAN
\Delta (m) is the vector of diagonal elements of W (m) 0
C 'f\Deltag is the softmax
operator, which maps a vector A into a vector B of the same size, with elements
and log P (m) denotes the elementwise logarithm of the transition matrix P (m) .
The first term of (9a) is the projection of the error in reconstructing the observation
onto the weights of state vector m-the more a particular setting of a
state vector can reduce this error, the larger its associated variational parameter.
The second term arises from the fact that the second order correlation hS (m)
evaluated under the variational distribution is a diagonal matrix composed of the
elements of ' (m)
t . The last two terms introduce dependencies forward and backward
in time. 5 Therefore, although the posterior distribution over the hidden variables is
approximated with a completely factorized distribution, the fixed point equations
couple the parameters associated with each node with the parameters of its Markov
blanket. In this sense, the fixed point equations propagate information along the
same pathways as those defining the exact algorithms for probability propagation.
The following may provide an intuitive interpretation of the approximation being
made by this distribution. Given a particular observation sequence, the hidden
state variables for the M Markov chains at time step t are stochastically coupled.
This stochastic coupling is approximated by a system in which the hidden variables
are uncorrelated but have coupled means. The variational or "mean-field" equations
solve for the deterministic coupling of the means that best approximates the
stochastically coupled system.
Each hidden state vector is updated in turn using (9a), with a time complexity
of O(TMK 2 ) per iteration. Convergence is determined by monitoring the KL
divergence in the variational distribution between successive time steps; in practice
convergence is very rapid (about 2 to 10 iterations of (9a)). Once the fixed point
equations have converged, the expectations required for the E step can be obtained
as a simple function of the parameters (equations (C.6)-(C.8) in Appendix C).
3.5. Structured variational inference
The approximation presented in the previous section factors the posterior probability
such that all the state variables are statistically independent. In contrast to
this rather extreme factorization, there exists a third approximation that is both
tractable and preserves much of the probabilistic structure of the original system. In
this scheme, the factorial HMM is approximated by M uncoupled HMMs as shown
in

Figure

(b). Within each HMM, efficient and exact inference is implemented
via the forward-backward algorithm. The approach of exploiting such tractable
substructures was first suggested in the machine learning literature by Saul and
Jordan (1996).
FACTORIAL HIDDEN MARKOV MODELS 255
Note that the arguments presented in the previous section did not hinge on the
the form of the approximating distribution. Therefore, more structured variational
approximations can be obtained by using more structured variational distributions
Q. Each such Q provides a lower bound on the log likelihood and can be used to
obtain a learning algorithm.
We write the structured variational approximation as
Y
Q(S (m)
Y
Q(S (m)
where ZQ is a normalization constant ensuring that Q integrates to one and
Q(S (m)
Y
Q(S (m)
Y
t;k
t;k
Y
t;k
Y
t;k
where the last equality follows from the fact that S (m)
is a vector with a 1 in one position
and 0 elsewhere. The parameters of this distribution are
t gthe
original priors and state transition matrices of the factorial HMM and a time-varying
bias for each state variable. Comparing equations (11a)-(11c) to equation
(1), we can see that the K \Theta 1 vector h (m)
t plays the role of the probability of
an observation (P (Y t jS t ) in (1)) for each of the K settings of S (m)
t . For example,
Q(S (m)
1jOE) corresponds to having an observation at time
that under state S (m)
1;j .
Intuitively, this approximation uncouples the M Markov chains and attaches to
each state variable a distinct fictitious observation. The probability of this fictitious
observation can be varied so as to minimize the KL divergence between Q and P .
Applying the same arguments as before, we obtain a set of fixed point equations
for h (m)
t that minimize KL(QkP ), as detailed in Appendix D:
h (m) new
ae
Y (m)
oe
where \Delta (m) is defined as before, and where we redefine the residual error to be
~
Y (m)
'6=m
Z. GHAHRAMANI AND M.I. JORDAN
The parameter h (m)
t obtained from these fixed point equations is the observation
probability associated with state variable S (m)
t in hidden Markov model m. Using
these probabilities, the forward-backward algorithm is used to compute a new set
of expectations for hS (m)
t i, which are fed back into (12a) and (12b). This algorithm
is therefore used as a subroutine in the minimization of the KL divergence.
Note the similarity between equations (12a)-(12b) and equations (9a)-(9b) for the
completely factorized system. In the completely factorized system, since hS (m)
t , the fixed point equations can be written explicitly in terms of the variational
parameters. In the structured approximation, the dependence of hS (m)
t i on h (m)
is computed via the forward-backward algorithm. Note also that (12a) does not
contain terms involving the prior, - (m) , or transition matrix, P (m) . These terms
have cancelled by our choice of approximation.
3.6. Choice of approximation
The theory of the EM algorithm as presented in Dempster et al. (1977) assumes
the use of an exact E step. For models in which the exact E step is intractable,
one must instead use an approximation like those we have just described. The
choice among these approximations must take into account several theoretical and
practical issues.
Monte Carlo approximations based on Markov chains, such as Gibbs sampling,
offer the theoretical assurance that the sampling procedure will converge to the
correct posterior distribution in the limit. Although this means that one can come
arbitrarily close to the exact E step, in practice convergence can be slow (especially
for multimodal distributions) and it is often very difficult to determine how close
one is to convergence. However, when sampling is used for the E step of EM, there
is a time tradeoff between the number of samples used and the number of EM
iterations. It seems wasteful to wait until convergence early on in learning, when
the posterior distribution from which samples are drawn is far from the posterior
given the optimal parameters. In practice we have found that even approximate
steps using very few Gibbs samples (e.g. around ten samples of each hidden
variable) tend to increase the true likelihood.
Variational approximations offer the theoretical assurance that a lower bound on
the likelihood is being maximized. Both the minimization of the KL divergence in
the E step and the parameter update in the M step are guaranteed not to decrease
this lower bound, and therefore convergence can be defined in terms of the bound.
An alternative view given by Neal and Hinton (1993) describes EM in terms of the
negative free energy, F , which is a function of the parameters, OE, the observations,
Y , and a posterior probability distribution over the hidden variables, Q(S):
where EQ denotes expectation over S using the distribution Q(S). The exact E
step in EM maximizes F with respect to Q given OE. The variational E steps used
FACTORIAL HIDDEN MARKOV MODELS 257
here maximize F with respect to Q given OE, subject to the constraint that Q is
of a particular tractable form. Given this view, it seems clear that the structured
approximation is preferable to the completely factorized approximation since it
places fewer constraints on Q, at no cost in tractability.
4. Experimental results
To investigate learning and inference in factorial HMMs we conducted two experi-
ments. The first experiment compared the different approximate and exact methods
of inference on the basis of computation time and the likelihood of the model obtained
from synthetic data. The second experiment sought to determine whether
the decomposition of the state space in factorial HMMs presents any advantage in
modeling a real time series data set that we might assume to have complex internal
structure-Bach's chorale melodies.
4.1. Experiment 1: Performance and timing benchmarks
Using data generated from a factorial HMM structure with M underlying Markov
models with K states each, we compared the time per EM iteration and the training
and test set likelihoods of five models:
ffl HMM trained using the Baum-Welch algorithm;
ffl Factorial HMM trained with exact inference for the E step, using a straight-forward
application of the forward-backward algorithm, rather than the more
efficient algorithm outlined in Appendix B;
ffl Factorial HMM trained using Gibbs sampling for the E step with the number
of samples fixed at 10 samples per variable; 6
ffl Factorial HMM trained using the completely factorized variational approxima-
tion; and
ffl Factorial HMM trained using the structured variational approximation.
All factorial HMMs consisted of M underlying Markov models with K states each,
whereas the HMM had K M states. The data were generated from a factorial HMM
structure with M state variables, each of which could take on K discrete values.
All of the parameters of this model, except for the output covariance matrix, were
sampled from a uniform [0; 1] distribution and appropriately normalized to satisfy
the sum-to-one constraints of the transition matrices and priors. The covariance
matrix was set to a multiple of the identity matrix
The training and test sets consisted of 20 sequences of length 20, where the observable
was a four-dimensional vector. For each randomly sampled set of parameters, a
separate training set and test set were generated and each algorithm was run once.
Z. GHAHRAMANI AND M.I. JORDAN
Fifteen sets of parameters were generated for each of the four problem sizes. Algorithms
were run for a maximumof 100 iterations of EM or until convergence, defined
as the iteration k at which the log likelihood L(k), or approximate log likelihood if
an approximate algorithm was used, satisfied
At the end of learning, the log likelihoods on the training and test set were computed
for all models using the exact algorithm. Also included in the comparison
was the log likelihood of the training and test sets under the true model that generated
the data. The test set log likelihood for N observation sequences is defined
as
log P (Y (n)
obtained by maximizing the log likelihood
over a training set that is disjoint from the test set. This provides a measure
of how well the model generalizes to a novel observation sequence from the same
distribution as the training data.
Results averaged over 15 runs for each algorithm on each of the four problem sizes
(a total of 300 runs) are presented in Table 1. Even for the smallest problem size
the standard HMM with K M states suffers from overfitting:
the test set log likelihood is significantly worse than the training set log likelihood.
As expected, this overfitting problem becomes worse as the size of the state space
increases; it is particularly serious for
For the factorial HMMs, the log likelihoods for each of the three approximate
EM algorithms were compared to the exact algorithm. Gibbs sampling appeared
to have the poorest performance: for each of the three smaller size problems its
log likelihood was significantly worse than that of the exact algorithm on both the
training sets and test sets (p ! 0:05). This may be due to insufficient sampling.
However, we will soon see that running the Gibbs sampler for more than 10 samples,
while potentially improving performance, makes it substantially slower than the
variational methods. Surprisingly, the Gibbs sampler appears to do quite well on
the largest size problem, although the differences to the other methods are not
statistically significant.
The performance of the completely factorized variational approximation was not
statistically significantly different from that of the exact algorithm on either the
training set or the test set for any of the problem sizes. The performance of the
structured variational approximation was not statistically different from that of the
exact method on three of the four problem sizes, and appeared to be better on one of
the problem sizes 0:05). Although this result may be a fluke
arising from random variability, there is another more interesting (and speculative)
explanation. The exact EM algorithm implements unconstrained maximization of
F , as defined in section 3.6, while the variational methods maximize F subject to
a constrained distribution. These constraints could presumably act as regularizers,
reducing overfitting.
There was a large amount of variability in the final log likelihoods for the models
learned by all the algorithms. We subtracted the log likelihood of the true generative
model from that of each trained model to eliminate the main effect of the randomly
sampled generative model and to reduce the variability due to training and test
sets. One important remaining source of variance was the random seed used in
FACTORIAL HIDDEN MARKOV MODELS 259

Table

1. Comparison of the factorial HMM on four problems of varying size. The negative log
likelihood for the training and test set, plus or minus one standard deviation, is shown for each
problem size and algorithm, measured in bits per observation (log likelihood in bits divided by
NT ) relative to the log likelihood under the true generative model for that data set. 7 True is
the true generative model (the log likelihood per symbol is defined to be zero for this model by
our measure); HMM is the hidden Markov model with K M states; Exact is the factorial HMM
trained using an exact E step; Gibbs is the factorial HMM trained using Gibbs sampling; CFVA
is the factorial HMM trained using the completely factorized variational approximation; SVA is
the factorial HMM trained using the structured variational approximation.
M K Algorithm Training Test
HMM 1.19 \Sigma 0.67 2.29 \Sigma 1.02
Exact 0.88 \Sigma 0.80 1.05 \Sigma 0.72
Gibbs 1.67 \Sigma 1.23 1.78 \Sigma 1.22
CFVA 1.06 \Sigma 1.20 1.20 \Sigma 1.11
SVA 0.91 \Sigma 1.02 1.04 \Sigma 1.01
HMM 0.76 \Sigma 0.67 9.81 \Sigma 2.55
Exact 1.02 \Sigma 1.04 1.26 \Sigma 0.99
Gibbs 2.21 \Sigma 0.91 2.50 \Sigma 0.87
CFVA 1.24 \Sigma 1.50 1.50 \Sigma 1.53
Exact 2.29 \Sigma 1.19 2.51 \Sigma 1.21
Gibbs 3.25 \Sigma 1.17 3.35 \Sigma 1.14
CFVA 1.73 \Sigma 1.34 2.07 \Sigma 1.74
Exact 4.23 \Sigma 2.28 4.49 \Sigma 2.24
Gibbs 3.63 \Sigma 1.13 3.95 \Sigma 1.14
CFVA 4.85 \Sigma 0.68 5.14 \Sigma 0.69
260 Z. GHAHRAMANI AND M.I. JORDAN
iterations of EM
-log
likelihood
(bits)
(a)
iterations of EM
-log
likelihood
(bits)
(b)
iterations of EM
-log
likelihood
(bits)
(c)
iterations of EM
-log
likelihood
(bits)
(d)

Figure

3. Learning curves for five runs of each of the four learning algorithms for factorial HMMs:
(a) exact; (b) completely factorized variational approximation; (c) structured variational approx-
imation; and (d) Gibbs sampling. A single training set sampled from the
size was used for all these runs. The solid lines show the negative log likelihood per observation
(in bits) relative to the true model that generated the data, calculated using the exact algorithm.
The circles denote the point at which the convergence criterion was met and the run ended. For
the three approximate algorithms, the dashed lines show an approximate negative log likelihood. 8
each training run, which determined the initial parameters and the samples used in
the Gibbs algorithm. All algorithms appeared to be very sensitive to this random
seed, suggesting that different runs on each training set found different local maxima
or plateaus of the likelihood (Figure 3). Some of this variability could be eliminated
by explicitly adding a regularization term, which can be viewed as a prior on the
parameters in maximuma posteriori parameter estimation. Alternatively, Bayesian
(or ensemble) methods could be used to average out this variability by integrating
over the parameter space.
The timing comparisons confirm the fact that both the standard HMM and the exact
are extremely slow for models with large state spaces (Fig-
FACTORIAL HIDDEN MARKOV MODELS 261
Time/iteration
HMM

Figure

4. Time per iteration of EM on a Silicon Graphics R4400 processor running Matlab.
ure 4). Gibbs sampling was slower than the variational methods even when limited
to ten samples of each hidden variable per iteration of EM. Since one pass of the
variational fixed point equations has the same time complexity as one pass of Gibbs
sampling, and since the variational fixed point equations were found to converge
very quickly, these experiments suggest that Gibbs sampling is not as competitive
time-wise as the variational methods. The time per iteration for the variational
methods scaled well to large state spaces.
4.2. Experiment 2: Bach chorales
Musical pieces naturally exhibit complex structure at many different time scales.
Furthermore, one can imagine that to represent the "state" of the musical piece
at any given time it would be necessary to specify a conjunction of many different
features. For these reasons, we chose to test whether a factorial HMM would provide
an advantage over a regular HMM in modeling a collection of musical pieces.
The data set consisted of discrete event sequences encoding the melody lines of
J. S. Bach's Chorales, obtained from the UCI Repository for Machine Learning
originally discussed in Conklin and Witten
(1995). Each event in the sequence was represented by six attributes, described
in

Table

2. Sixty-six chorales, with 40 or more events each, were divided into a
training set (30 chorales) and a test set (36 chorales). Using the first set, hidden
Markov models with state space ranging from 2 to 100 states were trained until
convergence (30 \Sigma 12 steps of EM). Factorial HMMs of varying sizes (K ranging
from 2 to 6; M ranging from 2 to were also trained on the same data. To
262 Z. GHAHRAMANI AND M.I. JORDAN

Table

2. Attributes in the Bach chorale data set. The key
signature and time signature attributes were constant over the
duration of the chorale. All attributes were treated as real
numbers and modeled as linear-Gaussian observations (4a).
Attribute Description Representation
pitch pitch of the event int [0; 127]
fermata event under fermata? binary
st start time of event int (1/16 notes)
dur duration of event int (1/16 notes)
approximate the E step for factorial HMMs we used the structured variational ap-
proximation. This choice was motivated by three considerations. First, for the size
of state space we wished to explore, the exact algorithms were prohibitively slow.
Second, the Gibbs sampling algorithm did not appear to present any advantages
in speed or performance and required some heuristic method for determining the
number of samples. Third, theoretical arguments suggest that the structured approximation
should in general be superior to the completely factorized variational
approximation, since more of the dependencies of the original model are preserved.
The test set log likelihoods for the HMMs, shown in Figure 5 (a), exhibited the
typical U-shaped curve demonstrating a trade-off between bias and variance (Ge-
man, Bienenstock, & Doursat, 1992). HMMs with fewer than 10 states did not
predict well, while HMMs with more than 40 states overfit the training data and
therefore provided a poor model of the test data. Out of the 75 runs, the highest
test set log likelihood per observation was \Gamma9:0 bits, obtained by an HMM with
hidden states. 9
The factorial HMM provides a more satisfactory model of the chorales from three
points of view. First, the time complexity is such that it is possible to consider
models with significantly larger state spaces; in particular, we fit models with up to
1000 states. Second, given the componential parametrization of the factorial HMM,
these large state spaces do not require excessively large numbers of parameters relative
to the number of data points. In particular, we saw no evidence of overfitting
even for the largest factorial HMM as seen in Figures 5 (c) & (d). Finally, this
approach resulted in significantly better predictors; the test set likelihood for the
best factorial HMM was an order of magnitude larger than the test set likelihood
for the best HMM, as Figure 5 (d) reveals.
While the factorial HMM is clearly a better predictor than a single HMM, it
should be acknowledged that neither approach produces models that are easily
interpretable from a musicological point of view. The situation is reminiscent of
that in speech recognition, where HMMs have proved their value as predictive
models of the speech signal without necessarily being viewed as causal generative
models of speech. A factorial HMM is clearly an impoverished representation of
FACTORIAL HIDDEN MARKOV MODELS 263
log
likelihood
(bits)
(a)
(b)
Size of state space
(c)
(d)

Figure

5. Test set log likelihood per event of the Bach chorale data set as a function of number of
states for (a) HMMs, and factorial HMMs with (b)
dashed line) and line). Each symbol represents a single run; the lines indicate the
mean performances. The thin dashed line in (b)-(d) indicates the log likelihood per observation
of the best run in (a). The factorial HMMs were trained using the structured approximation. For
both methods the true likelihood was computed using the exact algorithm.
musical structure, but its promising performance as a predictor provides hope that
it could serve as a step on the way toward increasingly structured statistical models
for music and other complex multivariate time series.
5. Generalizations of the model
In this section, we describe four variations and generalizations of the factorial HMM.
5.1. Discrete observables
The probabilistic model presented in this paper has assumed real-valued Gaussian
observations. One of the advantages arising from this assumption is that the
conditional density of a D-dimensional observation, P (Y t jS (1)
t ), can be
compactly specified through M mean matrices of dimension D \Theta K, and one D \Theta D
covariance matrix. Furthermore, the M step for such a model reduces to a set of
weighted least squares equations.
The model can be generalized to handle discrete observations in several ways.
Consider a single D-valued discrete observation Y t . In analogy to traditional HMMs,
the output probabilities could be modeled using a matrix. However, in the case of a
factorial HMM, this matrix would have D \Theta K M entries for each combination of the
state variables and observation. Thus the compactness of the representation would
be entirely lost. Standard methods from graphical models suggest approximating
such large matrices with "noisy-OR" (Pearl, 1988) or "sigmoid" (Neal, 1992) models
of interaction. For example, in the softmax model, which generalizes the sigmoid
model to D ? 2, P (Y t jS (1)
t ) is multinomial with mean proportional to
264 Z. GHAHRAMANI AND M.I. JORDAN
exp
. Like the Gaussian model, this specification is again com-
pact, using M matrices of size D \Theta K. (As in the linear-Gaussian model, the W (m)
are overparametrized since they can each model the overall mean of Y t , as shown in


Appendix

A.) While the nonlinearity induced by the softmax function makes both
the E step and M step of the algorithm more difficult, iterative numerical methods
can be used in the M step whereas Gibbs sampling and variational methods can
again be used in the E step (see Neal, 1992; Hinton et al., 1995; and Saul et al.,
1996, for discussions of different approaches to learning in sigmoid networks).
5.2. Introducing couplings
The architecture for factorial HMMs presented in Section 2 assumes that the underlying
Markov chains interact only through the observations. This constraint can
be relaxed by introducing couplings between the hidden state variables (cf. Saul &
Jordan, 1997). For example, if S (m)
t depends on S (m)
equation (3) is
replaced by the following factorization
Y
Similar exact, variational, and Gibbs sampling procedures can be defined for this
architecture. However, note that these couplings must be introduced with caution,
as they may result in an exponential growth in parameters. For example, the above
factorization requires transition matrices of size K 2 \Theta K. Rather than specifying
these higher-order couplings through probability transition matrices, one can introduce
second-order interaction terms in the energy (log probability) function. Such
terms effectively couple the chains without the number of parameters incurred by
a full probability transition matrix. In the graphical model formalism these correspond
to symmetric undirected links, making the model a chain graph. While
the Jensen, Lauritzen and Olesen (1990) algorithm can still be used to propagate
information exactly in chain graphs, such undirected links cause the normalization
constant of the probability distribution-the partition function-to depend on the
coupling parameters. As in Boltzmann machines (Hinton & Sejnowski, 1986), both
a clamped and an unclamped phase are therefore required for learning, where the
goal of the unclamped phase is to compute the derivative of the partition function
with respect to the parameters (Neal, 1992).
5.3. Conditioning on inputs
Like the hidden Markov model, the factorial HMM provides a model of the unconditional
density of the observation sequences. In certain problem domains, some of
the observations can be better thought of as inputs or explanatory variables, and
FACTORIAL HIDDEN MARKOV MODELS 265
the others as outputs or response variables. The goal, in these cases, is to model
the conditional density of the output sequence given the input sequence. In machine
learning terminology, unconditional density estimation is unsupervised while
conditional density estimation is supervised.
Several algorithms for learning in hidden Markov models that are conditioned on
inputs have been recently presented in the literature (Cacciatore & Nowlan, 1994;
Bengio & Frasconi, 1995; Meila & Jordan, 1996). Given a sequence of input vectors
g, the probabilistic model for an input-conditioned factorial HMM is
Y
\Theta
Y
Y
The model depends on the specification of P (Y t jS (m)
which are conditioned both on a discrete state variable and on a (possibly con-
tinuous) input vector. The approach used in Bengio and Frasconi's Input Output
HMMs (IOHMMs) suggests modeling P (S (m)
separate neural
networks, one for each setting of S (m)
. This decomposition ensures that a valid
probability transition matrix is defined at each point in input space if a sum-to-one
constraint (e.g., softmax nonlinearity) is used in the output of these networks.
Using the decomposition of each conditional probability into K networks, inference
in input-conditioned factorial HMMs is a straightforward generalization of the
algorithms we have presented for factorial HMMs. The exact forward-backward
algorithm in Appendix B can be adapted by using the appropriate conditional
probabilities. Similarly, the Gibbs sampling procedure is no more complex when
conditioned on inputs. Finally, the completely factorized and structured approximations
can also be generalized readily if the approximating distribution has a
dependence on the input similar to the model's. If the probability transition structure
not decomposed as above, but has a complex dependence
on the previous state variable and input, inference may become considerably more
complex.
Depending on the form of the input conditioning, the Maximization step of learning
may also change considerably. In general, if the output and transition probabilities
are modeled as neural networks, the M step can no longer be solved exactly
and a gradient-based generalized EM algorithm must be used. For log-linear
models, the M step can be solved using an inner loop of iteratively reweighted
least-squares (McCullagh & Nelder, 1989).
5.4. Hidden Markov decision trees
An interesting generalization of factorial HMMs results if one conditions on an
input X t and orders the M state variables such that S (m)
t depends on S (n)
t for
266 Z. GHAHRAMANI AND M.I. JORDAN
S (2)
S (2)
S (2)

Figure

6. The hidden Markov decision tree.

Figure

6). The resulting architecture can be seen as a probabilistic
decision tree with Markovian dynamics linking the decision variables. Consider how
this probabilistic model would generate data at the first time step, Given
the top node S (1)
can take on K values. This stochastically partitions X
space into K decision regions. The next node down the hierarchy, S (2), subdivides
each of these regions into K subregions, and so on. The output Y 1 is generated from
the input X 1 and the K-way decisions at each of the M hidden nodes. At the next
time step, a similar procedure is used to generate data from the model, except that
now each decision in the tree is dependent on the decision taken at that node in the
previous time step. Thus, the "hierarchical mixture of experts" architecture (Jordan
Jacobs, 1994) is generalized to include Markovian dynamics for the decisions.
Hidden Markov decision trees provide a useful starting point for modeling time
series with both temporal and spatial structure at multiple resolutions. We explore
this generalization of factorial HMMs in Jordan, Ghahramani, and Saul (1997).
6. Conclusion
In this paper we have examined the problem of learning for a class of generalized
hidden Markov models with distributed state representations. This generalization
provides both a richer modeling tool and a method for incorporating prior structural
information about the state variables underlying the dynamics of the system
generating the data. Although exact inference in this class of models is generally
intractable, we provided a structured variational approximation that can be computed
tractably. This approximation forms the basis of the Expectation step in an
EM algorithm for learning the parameters of the model. Empirical comparisons
to several other approximations and to the exact algorithm show that this approximation
is both efficient to compute and accurate. Finally, we have shown that
FACTORIAL HIDDEN MARKOV MODELS 267
the factorial HMM representation provides an advantage over traditional HMMs in
predictive modeling of the complex temporal patterns in Bach's chorales.


Appendix

A
The M step
The M step equations for each parameter are obtained by setting the derivatives
of Q with respect to that parameters to zero. We start by expanding Q using
equations (1)-(4b):
tr
C
tr
(log P (m) )hS (m)
log Z; (A.1)
where tr is the trace operator for square matrices and Z is a normalization term
independent of the states and observations ensuring that the probabilities sum to
one.
Setting the derivatives of Q with respect to the output weights to zero, we obtain
a linear system of equations for the W
0: (A.2)
Assuming Y t is a D\Theta1 vector, let S t be the MK \Theta1 vector obtained by concatenating
the S (m) vectors, and W be the D \Theta MK matrix obtained by concatenating the
W (m) matrices (of size D \Theta K). Then solving (A.2) results in
where y is the Moore-Penrose pseudo-inverse. Note that the model is overparameterized
since the D \Theta 1 means of each of the W (m) matrices add up to a single mean.
Using the pseudo-inverse removes the need to explicitly subtract this overall mean
from each W (m) and estimate it separately as another parameter.
To estimate the priors, we solve @Q=@- subject to the constraint that
they sum to one, obtaining
268 Z. GHAHRAMANI AND M.I. JORDAN
Similarly, to estimate the transition matrices we solve @Q=@P subject to
the constraint that the columns of P (m) sum to one. For element (i;
new
Finally, the re-estimation equations for the covariance matrix can be derived by
taking derivatives with respect to C \Gamma1
The first term arises from the normalization for the Gaussian density function: Z is
proportional to jCj T=2 and @jCj=@C Substituting (A.2) and re-organizing
we get
For equations reduce to the Baum-Welch re-estimation equations
for HMMs with Gaussian observables. The above M step has been presented for
the case of a single observation sequence. The extension to multiple sequences is
straightforward.


Appendix

Exact forward-backward algorithm
Here we specify an exact forward-backward recursion for computing the posterior
probabilities of the hidden states in a factorial HMM. It differs from a straightforward
application of the forward-backward algorithm on the equivalent K M state
HMM, in that it does not depend on a K M \Theta K M transition matrix. Rather, it
makes use of the independence of the underlying Markov chains to sum over M
transition matrices of size K \Theta K.
Using the notation fY - g r
t to mean the observation sequence Y
ff (1)
ff (M)
FACTORIAL HIDDEN MARKOV MODELS 269
Then we obtain the forward recursions
and
At the end of the forward recursions, the likelihood of the observation sequence is
the sum of the K M elements in ff T .
Similarly, to obtain the backward recursions we define
from which we obtain
The posterior probability of the state at time t is obtained by multiplying ff t and
This algorithm can be shown to be equivalent to the Jensen, Lauritzen and Olesen
algorithm for probability propagation in graphical models. The probabilities
are defined over collections of state variables corresponding to the cliques in the
equivalent junction tree. Information is passed forwards and backwards by summing
over the sets separating each neighboring clique in the tree. This results in
forward-backward-type recursions of order O(TMK M+1 ).
Using the ff t , fi t , and fl t quantities, the statistics required for the E step are
Z. GHAHRAMANI AND M.I. JORDAN


Appendix

Completely factorized variational approximation
Using the definition of the probabilistic model given by equations (1)-(4b), the
posterior probability of the states given an observation sequence can be written as
Z
where Z is a normalization constant ensuring that the probabilities sum to one and
Similarly, the probability distribution given by the variational approximation (7)-
(8) can be written as
expf\GammaH
where
log ' (m)
Using this notation, and denoting expectation with respect to the variational distribution
using angular brackets h\Deltai, the KL divergence is
Three facts can be verified from the definition of the variational approximation:
diagf' (m)
FACTORIAL HIDDEN MARKOV MODELS 271
where diag is an operator that takes a vector and returns a square matrix with
the elements of the vector along its diagonal, and zeros everywhere else. The KL
divergence can therefore be expanded to
log ' (m)
C
tr
C
trf' (m)
log P (m)
Taking derivatives with respect to ' (m)
t , we obtain
log ' (m)
C
\Gamma(log
where \Delta (m) is the vector of diagonal elements of W (m) 0
C c is a term
arising from log ZQ , ensuring that the ' (m)
t sum to one. Setting this derivative
equal to 0 and solving for ' (m)
t gives equation (9a).


Appendix

Structured approximation
For the structured approximation, HQ is defined as
log h (m)
Using (C.2), we write the KL divergence as
tr
C
tr
C diag
log Z: (D.2)
Z. GHAHRAMANI AND M.I. JORDAN
Since KL is independent of - (m) and P (m) , the first thing to note is that these
parameters of the structured approximation remain equal to the equivalent parameters
of the true system. Now, taking derivatives with respect to log h (n)
- , we
get
@ log h (n)
log h (m)
'6=m
C
@ log h (n)
The last term, which we obtained by making use of the fact that
@ log ZQ
@ log h (n)
cancels out the first term. Setting the terms inside the brackets in (D.3) equal to
zero yields equation (12a).

Acknowledgments

We thank Lawrence Saul for helpful discussions and Geoffrey Hinton for support.
This project was supported in part by a grant from the McDonnell-Pew Foundation,
by a grant from ATR Human Information Processing Research Laboratories, by a
gift from Siemens Corporation, and by grant N00014-94-1-0777 from the Office of
Naval Research. Zoubin Ghahramani was supported by a grant from the Ontario
Information Technology Research Centre.
Notes
1. For related work on inference in distributed state HMMs, see Dean and Kanazawa (1989).
2. In speech, neural networks are generally used to model P (S t jY t ); this probability is converted
to the observation probabilities needed in the HMM via Bayes rule.
3. If the columns of W (m) and W (n) are orthogonal for every pair of state variables, m and n, and
C is a diagonal covariance matrix, then the state variables will no longer be dependent given
the observation. In this case there is no "explaining away": each state variable is modeling the
variability in the observation along a different subspace.
4. A more Bayesian treatment of the learning problem, in which the parameters are also considered
hidden random variables, can be handled by Gibbs sampling by replacing the "M step"
with sampling from the conditional distribution of the parameters given the other hidden
variables (for example, see Tanner and Wong, 1987).
5. The first term is replaced by log - (m) for the second term does not appear for
6. All samples were used for learning; that is, no samples were discarded at the beginning of the
run. Although ten samples is too few to even approach convergence, it provides a run-time
roughly comparable to the variational methods. The goal was to see whether this "impatient"
Gibbs sampler would be able to compete with the other approximate methods.
FACTORIAL HIDDEN MARKOV MODELS 273
7. Lower values suggest a better probabilistic model: a value of one, for example, means that
it would take one bit more than the true generative model to code each observation vector.
Standard deviations reflect the variation due to training set, test set, and the random seed of
the algorithm. Standard errors on the mean are a factor of 3.8 smaller.
8. For the variational methods these dashed lines are equal to minus the lower bound on the
log likelihood, except for a normalization term which is intractable to compute and can vary
during learning, resulting in the apparent occasional increases in the bound.
9. Since the attributes were modeled as real numbers, the log likelihoods are only a measure of
relative coding cost. Comparisons between these likelihoods are meaningful, whereas to obtain
the absolute cost of coding a sequence, it is necessary to specify a discretization level.
10.This is analogous to the fully-connected Boltzmann machine with N units (Hinton & Sejnowski,
1986), in which every binary unit is coupled to every other unit using O(N 2 ) parameters, rather
than the O(2 N ) parameters required to specify the complete probability table.



--R

A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains.
An input-outputHMM architecture
Mixtures of controllers for jump linear and non-linear plants
Multiple viewpoint systems for music prediction.
Elements of information theory.
Applications of a general propagation algorithm for probabilistic expert systems.
A model for reasoning about persistence and causation.

Maximum likelihood from incomplete data via the EM algorithm.
Neural networks and the bias/variance dilemma.
Stochastic relaxation
Factorial learning and the EM algorithm.

Learning and relearning in Boltzmann machines.

Bayesian updating in recursive graphical models by local computations.

Hierarchical mixtures of experts and the EM algorithm.
Neural Computation
Stochastic simulation algorithms for dynamic probabilistic networks.
Hidden Markov models in computational biology: Applications to protein modeling
Local computations with probabilities on graphical structures and their application to expert systems.
Generalized linear models.
Learning fine motion by Markov mixtures of experts.
UCI Repository of machine learning databases
Connectionist learning of belief networks.
Probabilistic inference using Markov chain Monte Carlo methods (Technical Report CRG-TR-93-1)
A new view of the EM algorithm that justifies incremental and other variants.
Statistical field theory.
Probabilistic reasoning in intelligent systems: Networks of plausible inference.
CA: Morgan Kaufmann.
An Introduction to hidden Markov models.
Mixed memory Markov models.

Mean Field Theory for Sigmoid Belief Networks.
Journal of Artificial Intelligence Research
Boltzmann chains and hidden Markov models.
Exploiting tractable substructures in Intractable networks.

Probabilistic independence networks for hidden Markov probability models.
Hidden Markov model induction by Bayesian model merging.
The calculation of posterior distributions by data augmentation (with discussion).
bounds for convolutional codes and an asymptotically optimal decoding algorithm.
Mean field networks that learn to discriminate temporally distorted strings.
A minimum description length framework for unsupervised learning.

Received
--TR
Probabilistic reasoning in intelligent systems: networks of plausible inference
A model for reasoning about persistence and causation
Learning and relearning in Boltzmann machines
Elements of information theory
Connectionist learning of belief networks
Neural networks and the bias/variance dilemma
Hierarchical mixtures of experts and the EM algorithm
Probabilistic independence networks for hidden Markov probability models
Hidden Markov Model} Induction by Bayesian Model Merging
A minimum description length framework for unsupervised learning

--CTR
P. Xing , Michael I. Jordan , Stuart Russell, Graph partition strategies for generalized mean field inference, Proceedings of the 20th conference on Uncertainty in artificial intelligence, p.602-610, July 07-11, 2004, Banff, Canada
Ricardo Silva , Jiji Zhang , James G. Shanahan, Probabilistic workflow mining, Proceeding of the eleventh ACM SIGKDD international conference on Knowledge discovery in data mining, August 21-24, 2005, Chicago, Illinois, USA
Andrew Howard , Tony Jebara, Dynamical systems trees, Proceedings of the 20th conference on Uncertainty in artificial intelligence, p.260-267, July 07-11, 2004, Banff, Canada
Raul Fernandez , Rosalind W. Picard, Modeling drivers' speech under stress, Speech Communication, v.40 n.1-2, p.145-159, April
Robert A. Jacobs , Wenxin Jiang , Martin A. Tanner, Factorial hidden Markov models and the generalized backfitting algorithm, Neural Computation, v.14 n.10, p.2415-2437, October 2002
Terry Caelli , Andrew McCabe , Garry Briscoe, Shape tracking and production using hidden Markov models, Hidden Markov models: applications in computer vision, World Scientific Publishing Co., Inc., River Edge, NJ, 2001
Agnieszka Betkowska , Koichi Shinoda , Sadaoki Furui, Robust speech recognition using factorial HMMs for home environments, EURASIP Journal on Applied Signal Processing, v.2007 n.1, p.10-10, 1 January 2007
Yunhua Hu , Hang Li , Yunbo Cao , Dmitriy Meyerzon , Qinghua Zheng, Automatic extraction of titles from general documents using machine learning, Proceedings of the 5th ACM/IEEE-CS joint conference on Digital libraries, June 07-11, 2005, Denver, CO, USA
Yunhua Hu , Hang Li , Yunbo Cao , Li Teng , Dmitriy Meyerzon , Qinghua Zheng, Automatic extraction of titles from general documents using machine learning, Information Processing and Management: an International Journal, v.42 n.5, p.1276-1293, September 2006
Tony Jebara , Risi Kondor , Andrew Howard, Probability Product Kernels, The Journal of Machine Learning Research, 5, p.819-844, 12/1/2004
Fine , Yoram Singer , Naftali Tishby, The Hierarchical Hidden Markov Model: Analysis and Applications, Machine Learning, v.32 n.1, p.41-62, July 1998
Charles Sutton , Khashayar Rohanimanesh , Andrew McCallum, Dynamic conditional random fields: factorized probabilistic models for labeling and segmenting sequence data, Proceedings of the twenty-first international conference on Machine learning, p.99, July 04-08, 2004, Banff, Alberta, Canada
Wang , Nan-Ning Zheng , Yan Li , Ying-Qing Xu , Heung-Yung Shum, Learning kernel-based HMMs for dynamic sequence synthesis, Graphical Models, v.65 n.4, p.206-221, July
Jie Tang , Hang Li , Yunbo Cao , Zhaohui Tang, Email data cleaning, Proceeding of the eleventh ACM SIGKDD international conference on Knowledge discovery in data mining, August 21-24, 2005, Chicago, Illinois, USA
Cen Li , Gautam Biswas, A Bayesian approach for structural learning with hidden Markov models, Scientific Programming, v.10 n.3, p.201-219, August 2002
Sophie Deneve, Bayesian spiking neurons i: Inference, Neural Computation, v.20 n.1, p.91-117, January 2008
Lawrence K. Saul , Michael I. Jordan, Mixed Memory Markov Models: Decomposing Complex Stochastic Processes as Mixtures of Simpler Ones, Machine Learning, v.37 n.1, p.75-87, Oct. 1999
Yong Cao , Petros Faloutsos , Frdric Pighin, Unsupervised learning for speech motion editing, Proceedings of the ACM SIGGRAPH/Eurographics symposium on Computer animation, July 26-27, 2003, San Diego, California
Hung H. Bui , Svetha Venkatesh , Geoff West, Tracking and surveillance in wide-area spatial environments using the abstract hidden Markov model, Hidden Markov models: applications in computer vision, World Scientific Publishing Co., Inc., River Edge, NJ, 2001
R. Anderson , Pedro Domingos , Daniel S. Weld, Relational Markov models and their application to adaptive web navigation, Proceedings of the eighth ACM SIGKDD international conference on Knowledge discovery and data mining, July 23-26, 2002, Edmonton, Alberta, Canada
Tom ingliar , Milo Hauskrecht, Noisy-OR Component Analysis and its Application to Link Analysis, The Journal of Machine Learning Research, 7, p.2189-2213, 12/1/2006
Martin V. Butz, Kernel-based, ellipsoidal conditions in the real-valued XCS classifier system, Proceedings of the 2005 conference on Genetic and evolutionary computation, June 25-29, 2005, Washington DC, USA
Ying Wu , Thomas S. Huang, Robust Visual Tracking by Integrating Multiple Cues Based on Co-Inference Learning, International Journal of Computer Vision, v.58 n.1, p.55-71, June 2004
Michael I. Jordan , Zoubin Ghahramani , Tommi S. Jaakkola , Lawrence K. Saul, An Introduction to Variational Methods for Graphical Models, Machine Learning, v.37 n.2, p.183-233, Nov.1.1999
Andrea Torsello , Antonio Robles-Kelly , Edwin R. Hancock, Discovering Shape Classes using Tree Edit-Distance and Pairwise Clustering, International Journal of Computer Vision, v.72 n.3, p.259-285, May       2007
Charles Sutton , Andrew McCallum , Khashayar Rohanimanesh, Dynamic Conditional Random Fields: Factorized Probabilistic Models for Labeling and Segmenting Sequence Data, The Journal of Machine Learning Research, 8, p.693-723, 5/1/2007
H. Attias, Independent factor analysis, Neural Computation, v.11 n.4, p.803-851, May 15, 1999
Jinhai Cai , Zhi-Qiang Liu, Hidden Markov Models with Spectral Features for 2D Shape Recognition, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.23 n.12, p.1454-1458, December 2001
John Binder , Daphne Koller , Stuart Russell , Keiji Kanazawa, Adaptive Probabilistic Networks with Hidden Variables, Machine Learning, v.29 n.2-3, p.213-244, Nov./Dec. 1997
Xiangdong An , Dawn Jutla , Nick Cercone, Privacy intrusion detection using dynamic Bayesian networks, Proceedings of the 8th international conference on Electronic commerce: The new e-commerce: innovations for conquering current barriers, obstacles and limitations to conducting successful business on the internet, August 13-16, 2006, Fredericton, New Brunswick, Canada
Akio Utsugi, Ensemble of Independent Factor Analyzers with Application to Natural Image Analysis, Neural Processing Letters, v.14 n.1, p.49-60, August 2001
Zoubin Ghahramani, An introduction to hidden Markov models and Bayesian networks, Hidden Markov models: applications in computer vision, World Scientific Publishing Co., Inc., River Edge, NJ, 2001
Cristian Sminchisescu , Atul Kanaujia , Dimitris Metaxas, Conditional models for contextual human motion recognition, Computer Vision and Image Understanding, v.104 n.2, p.210-220, November 2006
Hichem Snoussi , Ali Mohammad-Djafari, Bayesian Unsupervised Learning for Source Separation with Mixture of Gaussians Prior, Journal of VLSI Signal Processing Systems, v.37 n.2-3, p.263-279, June-July 2004
Inna Stainvas , David Lowe, A Generative Probabilistic Oriented Wavelet Model for Texture Segmentation, Neural Processing Letters, v.17 n.3, p.217-238, June
Russell Greiner , Christian Darken , N. Iwan Santoso, Efficient reasoning, ACM Computing Surveys (CSUR), v.33 n.1, p.1-30, March 2001
mean field theory;bayesian networks;hidden markov models;graphical models;time series;EM algorithm