Editorial Changes

"Footnotes" section added to end of the paper to collect together all of the footnotes that originally appeared at the foot of various pages. Hyperlinks to these footnotes have also been added to the body of the paper.

Notation "\(TraditionalForm\`■\_\(■\ ■\)\)" changed to "\(TraditionalForm\`■\_\(■, ■\)\)" (and, analogously, higher rank notations) throughout the paper.

"ie" changed to "i.e." throughout the paper.

TYPO"Metropolis N, Rosenbluth A W, Rosenbluth M N and Teller A H" changed to "Metropolis N, Rosenbluth A W, Rosenbluth M N, Teller A H and Teller E".

TYPOTYPOTYPOTYPOTYPOTYPOTYPOTYPOTYPOTYPOTYPOEquation reformatted to make it more readable.

TYPO".<footnote link>." changed to "<footnote link>.".

TYPO"higher order interactions" changed to "higher order interactions.".

TYPO"∂G/∂(Subscript[λ, i])" changed to "∂(G/∂)(Subscript[λ, i])" (twice).

TYPO"∂G/∂(Subscript[λ, i])" changed to "∂(G/∂)(Subscript[λ, i])".

TYPOTYPOFigure moved to a more appropriate place in the text.

TYPO"1 ≤ t ≤ τ" changed to "1 ≤ t < τ".

TYPO"(∑_(h=0)^((Subscript[N, h]) - 1))" changed to "\(TraditionalForm\`∑\+\(h\_1, h\_2, ⋯, h\_T = 0\)\%\(N\_h - 1\)\)" (twice).

TYPO"(∑_(x=0)^((Subscript[N, x]) - 1))" changed to "\(TraditionalForm\`∑\+\(x\_1, x\_2, ⋯, x\_T = 0\)\%\(N\_x - 1\)\)", and "(∑_(h=0)^((Subscript[N, h]) - 1))" changed to "\(TraditionalForm\`∑\+\(h\_1, h\_2, ⋯, h\_T = 0\)\%\(N\_h - 1\)\)".

TYPO"∂G/∂(Subscript[λ, i])" changed to "∂(G/∂)(Subscript[λ, i])", and "equation (12))" changed to "equation (12)".

TYPOTYPO"(∑_(x=0)^((Subscript[N, x]) - 1))" changed to "\(TraditionalForm\`∑\+\(x\_1, x\_2, ⋯, x\_T = 0\)\%\(N\_x - 1\)\)", and "(∑_(h=0)^((Subscript[N, h]) - 1))" changed to "\(TraditionalForm\`∑\+\(h\_1, h\_2, ⋯, h\_T = 0\)\%\(N\_h - 1\)\)" (twice).

TYPO"∂G/∂(Subscript[λ, i, j]^((hh)))" changed to "((∂G)/(∂(Subscript[λ, i, j]^((hh)))))", and "∂G/∂(Subscript[λ, i, j]^((hx)))" changed to "((∂G)/(∂(Subscript[λ, i, j]^((hx)))))".

TYPO"(Subscript[D, i, j]^((hx)))(τ)\({∏\+\(t = 1\)\%\(τ - 1\)S\_\(h\_t, x\_t, h\_\(t + 1\)\)\%\((+)\)}\)" changed to "(Subscript[D, i, j]^((hh)))(τ)\({∏\+\(t = 1\)\%\(τ - 1\)S\_\(h\_t, x\_t, h\_\(t + 1\)\)\%\((+)\)}\)".

TYPO"∂G/∂(Subscript[λ, i, j]^((hh)))" changed to "((∂G)/(∂(Subscript[λ, i, j]^((hh)))))".

TYPOTYPO"(1/(Z(x)))(∑_(x∈training set))" changed to "(∑_(x∈training set))(1/(Z(x)))".

TYPO"∂G/∂(Subscript[λ, i, j]^((hx)))" changed to "((∂G)/(∂(Subscript[λ, i, j]^((hx)))))".

TYPO"1/(Subscript[S, i, j]^((hx)))" changed to "(1/(Subscript[S, i, j]^((hx))))".

TYPO"∂G/∂(Subscript[λ, i, 0]^((hh)))" changed to "((∂G)/(∂(Subscript[λ, i, 0]^((hh)))))", and "∂G/∂(Subscript[λ, i, 0]^((hx)))" changed to "((∂G)/(∂(Subscript[λ, i, 0]^((hx)))))".

TYPO"an tandem" changed to "in tandem".

TYPO"(e^(-(Subscript[λ, i, j]^b)))" changed to "(e^(-(Subscript[((λ^′)), i, j]^b)))".

TYPO"[2]" changed to "[2].".

TYPO"log det(2π ((C')^(-1))(s))" changed to "log det((1/(2  π)) ((C')^(-1))(s))".

TYPO"+constant" appended to the right hand side of the inequality.

The Gibbs Machine applied to hidden Markov model problems. Part 1: Basic theory

Stephen P Luttrell
Pattern Processing Principles
SP4 division, RSRE
Malvern, WORCS, WR14 3PS

This paper appeared as SP4 Research Note, No. 99, 10 October 1989.

Copyright © Controller HMSO, London, 1989.

Abstract

I show how a hidden Markov model can be expressed as a Gibbs distribution. I review my Gibbs distribution training algorithm (a "Gibbs Machine [9] rather than a "Boltzmann Machine" [11]), which I use to perform gradient ascent on the relative entropy between the Gibbs distribution and the data distribution. I demonstrate how this reduces to elementary matrix computations of exactly the same form as encountered in the Baum-Welch re-estimation method [2]. Although this toy problem is amenable to Baum-Welch re-estimation, the same cannot be said of non-tree-like Markov models. In such cases I propose that a hybrid Baum-Welch/Gibbs Machine optimisation scheme should be used.

Bayesian information processing has recently enjoyed a belated [7] increase in popularity in diverse fields. Gibbs distribution (or, equivalently, Markov random field) models are now used extensively in image modelling, where Monte Carlo techniques such as simulated annealing can be used to locate approximately the maxima of posterior probabilities. Hidden Markov models are used extensively in speech modelling, where dynamic programming can be used to discover the maxima of posterior probabilities. The Bayesian approach to inference is robust in the sense that it provides the only consistent means of information processing [5].

Adaptive training of Markov-type models has also undergone a recent interest in popularity, in line with the resurgence of interest in "neural networks"FOOTNOTE. Markov random field model training has received rather patchy attention in the literature. The best example is the celebrated Boltzmann machine training algorithm [1], which I developed into the Gibbs machine training algorithm [9] - this can be used to train any type of Markov model with hidden variables, arbitrary interactions, and arbitrary state spaces, by hill climbing on a relative entropy criterion. On the other hand, hidden Markov model training has received much systematic attention over the years [2, 8], especially in the field of speech modelling. In this case, a re-estimation method can be used to train the model to maximise the same relative entropy criterion that is used for Markov random field training.

Gibbs machine hill-climbing and Markov model re-estimation are two sides of the same coin - both maximise a relative entropy criterion. The problem is to unify these two approaches. Thus I will formulate a standard hidden Markov model as a Gibbs distribution (and vice-versa), then show that the Gibbs distribution training algorithm can be expressed in terms of the same forwards/backwards computations that arise during Baum-Welch re-estimation. This demonstrates the detailed connection between the two approaches for the simplest case of interest - the hidden Markov model. The two approaches can then be applied in a hybrid fashion to more complicated Markov model training problems.

A particularly desirable advantage of the Gibbs distribution approach is its ability to build up a Markov model systematically using n-point interactions between the model variables. This permits one to construct a "cheap and cheerful" representation of the mutual correlations between n variables by using functions of only 2 variables. The simplest model of interest where this trick can be used to advantage is the standard hidden Markov model with direct causal influence between successive emitted symbols. This type of model has a wide range of applications in fields as diverse as radar modelling (correlated clutter) and speech modelling (co-articulation).

Gibbs distributions

In this section 1 give a resumeé of Gibbs distributions and the Hammersley-Clifford expansion. I strongly recommend that the interested reader refers to [3] for further information.

. The Gibbs distribution

Denote the state of a system as x and its probability of occurrence as P(x). This probability can always be written in the form

P(x)≡(1/Z)(e^(-U(x)))
Z≡∫dx (e^(-U(x)))

where U(x)=-log(P(x))+constant. This parameterisation of P(x) is called a Gibbs distributionFOOTNOTE.

. The Hammersley-Clifford expansion

A Gibbs distribution can always be written using a U(x) of the form [3]

U(x)=(∑_i) (Subscript[G, i]^((1)))((Subscript[x, i])) 
              (Subscript[x, i]) + (∑_(i < j))(Subscript[G, i, j]^((2)))((Subscript[x, i]), 
                (Subscript[x, j])) (Subscript[x, i]) 
              (Subscript[x, j]) + (∑_(i < j < k))(Subscript[G, i, j, k]^((3)))((Subscript[x, i]), 
                (Subscript[x, j]), (Subscript[x, k])) (Subscript[x, i]) (Subscript[x, j]) (Subscript[x, k]) + ⋯

where I assume that the (Subscript[x, i]) each have (Subscript[N, x]) possible states, thus (Subscript[x, i]) ∈ {0, 1, 2, ⋯, 
          (Subscript[N, x]) - 1}. This is known as the Hammersely-Clifford expansion. In the Hammersley-Clifford expansion every n-tuple (n ∈ {1, 2, 3, ⋯}) of components extracted from x has its own unique piece of U(x), except for those n-tuples that have one or more components equal to zero. The Harnmersley-Clifford expansion provides a systematic expansion of U(x) by order of statistical properties: 1-tuples, 2-tuples, etcFOOTNOTE. Although equation () demonstrates a dislike for the zero state, the important point is that the parameterisation of U(x) has exactly the right number of parameters in the right places. This could be achieved in other less easily stated ways.

A special case of the Hammersley-Clifford expansion is used to express U(x) for a Boltzmann machine [1]. Thus the components (Subscript[x, i]) of the Boltzmann machine state x are restricted to be binary variables ((Subscript[x, i]) ∈ {0, 1}), and only the 1-tuple and 2-tuple terms of the Hammersley-Clifford expansion are used, whence

(Subscript[U, BM])(x)=(∑_i) (Subscript[a, i]) 
              (Subscript[x, i]) + (∑_(i < j))(Subscript[w, i, j]) (Subscript[x, i]) (Subscript[x, j])

Note that since (Subscript[x, i]) ∈ {0, 1} the coefficients (Subscript[a, i]) and (Subscript[w, i, j]) need not have an explicit dependence on x.

I shall usually represent the Hammersley-Clifford expansion in the compact form

U(x)=(∑_i)(Subscript[λ, i])(Subscript[U, i])(x)

where the (Subscript[U, i])(x) are pieces taken from the right hand side of the Hammersley-Clifford expansion, and the (Subscript[λ, i]) are weighting coefficients. I shall also use the following representation of the Hammersley-Clifford expansion

U(x)=\(∑\+i\(∑\+\(a = 1\)\%\(N\_x - 1\)\(λ\_\(i, 
                    a\)\%\((1)\)\) δ\_\(x\_i, a\)\) + ∑\+\(i < j\)\(∑\+\(a, b = 1\)\%\(N\_x - 1\)\(λ\_\(i, j, a, 
                    b\)\%\((2)\)\) \(δ\_\(x\_i, 
                    a\)\) δ\_\(x\_j, b\)\) + ∑\+\(i < j < k\)\(∑\+\(a, b, c = 1\)\%\(N\_x - 1\)\(λ\_\(i, j, k, a, b, 
                    c\)\%\((3)\)\) \(δ\_\(x\_i, 
                    a\)\) \(δ\_\(x\_j, 
                    b\)\) δ\_\(x\_k, c\)\) + ⋯\)

where (Subscript[δ, m, n]) is a Kronecker delta. This expression provides a separate contribution to U(x) from each individual n-tuple state, and thus provides very fine control over the form of U(x).

. Hidden variables

The system state x might be only part of the state of a larger system whose joint state is X=(x,h), where h is the state of a set of unobserved, or hidden, variables. In this case, I shall call x the visible variables, to emphasise the distinction between x and h. Write P(x,h) as a Gibbs distribution U(x,h)

P(x,h)≡(1/Z)(e^(-U(x,h)))
Z≡∫dx dh (e^(-U(x,h)))

whence P(x) is a marginal probability given by

P(x)≡∫dh P(x,h)

Strictly speaking, the Boltzmann machine is a hidden variables model, so (Subscript[U, BM])(x) should be replaced by (Subscript[U, BM])(x,h), where

(Subscript[U, BM])(x,h)≡(∑_i) (Subscript[a, i]^((x))) (Subscript[x, i])+(∑_(i < j))(Subscript[w, i, j]^((xx)))(Subscript[x, i])(Subscript[x, j])+(∑_i) (Subscript[a, i]^((h))) (Subscript[h, i])+(∑_(i < j))(Subscript[w, i, j]^((hh)))(Subscript[h, i])(Subscript[h, j])+(∑_(i < j))(Subscript[w, i, j]^((hx)))(Subscript[h, i])(Subscript[x, j])

where the (Subscript[h, i]) are both coupled to each other and to the (Subscript[x, i]). A particular advantage of models with hidden variables is the ability of hidden variables (together with their interactions with the visible variables) to capture high order correlations amongst the visible variables. In a sense, all models have hidden variables because one observes superficial quantities in experiments (call these the visible variables), and then one seeks an underlying explanation (call these the hidden variables) for those observationsFOOTNOTE.

The Gibbs machine training algorithm

In this section I shall present a Gibbs distribution form of the Boltzmann machine training algorithm [1]. 1 originally presented this derivation in [9], and called it the Gibbs machine (training algorithm)FOOTNOTE.

. The relative entropy optimisation criterion

I wish to build a model Q(x) of the training set P(x), so the appropriate cost function to maximise is relative entropy G defined as

G≡-∫dx P(x)log[((P(x))/(Q(x)))]≤0

with G=0 iff P(x)=Q(x) ∀ x. G can be interpreted as follows. Suppose one draws a large number M of samples of x at random from the training set P(x), then (e^(M G)) is the probability that these M samples could have been drawn at random from the model Q(x). Thus (e^G) is the per-sample probability that Q(x) produces samples that look as if they had come from P(x)FOOTNOTE.

I shall construct the model Q(x) as a Gibbs distribution with hidden variables. Thus

Q(x)≡(1/Z)∫dh (e^(-λ.U(x,h)))

whence

G=∫dx P(x) log[(1/Z)∫dh (e^(-λ.U(x,h)))]+constant

The constant term depends only on the true system probability P(x) over which one has no control.

. Gradient of the relative entropy

I shall now derive a gradient ascent scheme for maximising G with respect to the coefficient vector λ, which will lead to the Gibbs machine training algorithm [9]. Thus U(x,h)

∂(G/∂)(Subscript[λ, i])=∫dx P(x) [(1/Z)∫(dx)'dh (Subscript[U, i])(x',h)(e^(-λ.U(x',h)))-((∫dh (Subscript[U, i])(x,h)(e^(-λ.U(x,h))))/(∫dh (e^(-λ.U(x,h)))))]
=∫dx P(x) [∫(dx)'dh (Subscript[U, i])(x',h)Q(x',h)-((∫dh (Subscript[U, i])(x,h)Q(x,h))/(Q(x)))]
=∫dx P(x) [(Subscript[\(⟨U\_i⟩\), Q(x,h)])-(Subscript[\(⟨U\_i⟩\), Q(h|x)])]
=(Subscript[⟨(Subscript[U, i])⟩, Q(x,h)])-(Subscript[⟨(Subscript[U, i])⟩, Q(h|x)P(x)])

The notation \(TraditionalForm\`⟨⋯⟩\) is defined as follows

(Subscript[\(⟨\([\)\(⋯\)\(]\)⟩\), ρ(X)])≡∫dX ρ(X) [⋯]

where ρ(X) is any probability measure. In equation () ∂(G/∂)(Subscript[λ, i]) is the difference between the averages of (Subscript[U, i]) obtained under two distinct conditions: the model Gibbs distribution Q(x,h), and a constrained Gibbs distribution Q(h|x)P(x). In Boltzmann machine terminology the Q(x,h)-average is called free because the system state is described by an unconstrained Gibbs distribution, whereas the Q(h|x)P(x)-average is called clamped because part of the system state x is described in a way that is determined externally by P(x) whilst the rest of the state h is described by a conditional Gibbs distribution. The beauty of the expression for ∂(G/∂)(Subscript[λ, i]) is that it allows one to optimise all of the Gibbs distribution parameters (Subscript[λ, i]) - not just those concerned directly with the observed part x of the system state. This powerful approach thus allows one to build hidden variable models Q(x,h) whose marginal distribution Q(x) approximates (in the G measure sense) the training set distribution P(x).

Practical implementations of the ∂(G/∂)(Subscript[λ, i]) computation use a Monte Carlo technique (in the form of the Metropolis algorithmFOOTNOTE [10]) to provide samples from Q(x,h) and Q(h|x), whilst the training set provides samples from P(x). There are many variations on a theme here. For instance, the Boltzmann machine is usually trained with some sort of simulated annealing to circumvent the tendency of the Metropolis algorithm to get trapped sampling in the vicinity of local maxima of Q(x,h) (or Q(h|x)). It is important to recognise that the Metropolis algorithm and simulated annealing are nothing more than algorithmic tricksFOOTNOTE.

Hidden Markov model and Gibbs distribution parameterisations

In this section I shall derive the Gibbs distribution formulation of hidden Markov model training. My derivation cures major problems that would arise if the approach discussed in [4] were to be attemptedFOOTNOTE. My approach is similar to, but more general than, that given in [11]FOOTNOTE.

. Hidden Markov model as a Gibbs distribution

A simple (first order) Markov chain model has a joint state space defined as h≡((Subscript[h, 1]), (Subscript[h, 2]), ⋯, (Subscript[h, T])) with a set of conditional probabilities Q((Subscript[h, t + 1]) | (Subscript[h, t])), t = 1, 2, ⋯, 
      T. Note that I use the hidden variable notation h for convenience, and I assume that the conditional probabilities are not dependent on the "time" index t. The probability that joint

Q(h)={(∏_(t = 1)^(T - 1))Q((Subscript[h, t + 1]) | (Subscript[h, t]))} Q(
            (Subscript[h, 1]))

where Q((Subscript[h, 1])) is needed to prime the conditional probability chain. This may, or may not, be equal to the equilibrium state probability.

A hidden Markov model may easily be obtained by forming a joint state space defined as (x,h) and a further set of conditional probabilities Q((Subscript[x, t]) | (Subscript[h, t])), t = 1, 2, ⋯, T. Thus each (Subscript[h, t]) can be thought of as producing (Subscript[h, t + 1]) and (Subscript[x, t]), where each is produced independently of the other (given that (Subscript[h, t]) is known). The probability that joint state (x,h) occurs is given by

Q(x,h)={(∏_(t = 1)^T) Q((Subscript[x, t]) | (Subscript[h, t]))}Q(h)

The expression in equation () is the starting point for standard derivations of hidden Markov model re-estimation. Invariably, one chooses the (Subscript[h, t]) to be discrete, but chooses the (Subscript[x, t]) to be either discrete or continuous. Without loss of generalityFOOTNOTE, I choose both to be discrete with (Subscript[x, t]) ∈ {0, 1, 2, ⋯, 
          (Subscript[N, x]) - 1} and (Subscript[h, t]) ∈ {0, 1, 2, ⋯, 
          (Subscript[N, h]) - 1}, 1 < t < T.

A Gibbs distribution formulation of Q(x,h) is then

Q(x,h)=(e^(-U(x,h)))
U(x,h)=-log[
              Q ((Subscript[h, 1]))] - (∑_(t = 1)^(T - 1))log[
              Q ((Subscript[h, t + 1]) | (Subscript[h, t]))] - (∑_(t = 1)^T) log[
              Q ((Subscript[x, t]) | (Subscript[h, t]))]

This is a tautology, because I have merely placed the conditional probability inside an exponential function. I have made this replacement to demonstrate that the standard hidden Markov model formulation maps directly onto the language of Gibbs distributions.

. Gibbs distribution as a hidden Markov model

A general Gibbs distribution with hidden variables takes the form shown in equation (). I shall parameterise this in the style of equation () to yield

λ.U(x,h)=(∑_(t = 1)^T)[\(∑\+\(i = 1\)\%\(N\_x - 1\)\(λ\_i\%\((x)\)\) δ\_\(x\_t, i\)\)+\(∑\+\(i = 1\)\%\(N\_h - 1\)\(λ\_i\%\((h)\)\) δ\_\(h\_t, i\)\)+(∑_(i = 1)^((Subscript[N, h]) - 1))(∑_(j = 1)^((Subscript[N, x]) - 1))(Subscript[λ, i, j]^((hx)))(Subscript[δ, (Subscript[h, t]), 
                          i])(Subscript[δ, (Subscript[x, t]), j])]+(∑_(t = 1)^(T - 1))(∑_(i, j = 1)^((Subscript[N, h]) - 1))(Subscript[λ, i, j]^((hh)))(Subscript[δ, (Subscript[h, t]), 
                    i])(Subscript[δ, (Subscript[h, t + 1]), j])

where I retain only the 1-point and 2-point terms to correspond to the hidden Markov model in §.. Note that I have assumed time shift invariance, so the λ-coefficients do not depend on T. This use of symmetry to simplify the model leads to a "spread network" [4] in which the model structure is repeated identically at each time slice. Figure  shows a piece of a hidden Markov model corresponding to equation (), with its 1-point terms represented as small circles, and its 2-point terms represented as ovals each of which encloses a pair of 1-point terms.

A hidden Markov model with all 1-point and 2-point terms shown.

It is both instructive and theoretically convenient to reparameterise equation () to make clear its relationship to equation (). Diagrammatically we must extract pieces of Figure  that look like conditional probabilities. I show these in Figure . Figure (a) corresponds to the log[Q((Subscript[h, t + 1]) | (Subscript[h, t]))] piece of equation (), whereas Figure (b) corresponds to the log[Q((Subscript[x, t]) | (Subscript[h, t]))] piece.

The conditional probability pieces of a hidden Markov model.

Now that it is clear from Figure  what I intend to do to equation (), 1 can proceed with confidence to rearrange it as follows

λ.U(x,h)=\(∑\+\(i = 1\)\%\(N\_h - 1\)\(λ\_i\%\((h)\)\) δ\_\(h\_1, i\)\)+(∑_(t = 1)^(T - 1))(∑_(i = 0)^((Subscript[N, h]) - 1))(∑_(j = 1)^((Subscript[N, h]) - 1))((Subscript[λ, j]^((h)))+(Subscript[λ, i, j]^((hh)))(1 - (Subscript[δ, i, 0])))(Subscript[δ, (Subscript[h, t]), 
                      i])(Subscript[δ, (Subscript[h, t + 1]), j])+(∑_(t = 1)^T)(∑_(i = 0)^((Subscript[N, h]) - 1))(∑_(j = 1)^((Subscript[N, x]) - 1))((Subscript[λ, j]^((x)))+(Subscript[λ, i, j]^((hx)))(1 - (Subscript[δ, i, 0])))(Subscript[δ, (Subscript[h, t]), i])(Subscript[δ, (Subscript[x, t]), 
                      j])

Redefine the notation as follows

(Subscript[λ, j]^((h)))+(Subscript[λ, i, j]^((hh)))(1 - (Subscript[δ, i, 0]))⟶(Subscript[λ, i, j]^((hh)))
(Subscript[λ, j]^((x)))+(Subscript[λ, i, j]^((hx)))(1 - (Subscript[δ, i, 0]))⟶(Subscript[λ, i, j]^((hx)))
(Subscript[λ, i]^((h)))⟶(Subscript[λ, 0, i]^((hh)))

whence

λ.U(x,h)=(∑_(i = 1)^((Subscript[N, h]) - 1))(Subscript[λ, 0, i]^((hh)))(Subscript[δ, (Subscript[h, 1]), 
                  i])+(∑_(t = 1)^(T - 1))(∑_(i = 0)^((Subscript[N, h]) - 1))(∑_(j = 1)^((Subscript[N, h]) - 1))(Subscript[λ, i, j]^((hh)))(Subscript[δ, (Subscript[h, t]), 
                      i])(Subscript[δ, (Subscript[h, t + 1]), j])+(∑_(t = 1)^T)(∑_(i = 0)^((Subscript[N, h]) - 1))(∑_(j = 1)^((Subscript[N, x]) - 1))(Subscript[λ, i, j]^((hx)))(Subscript[δ, (Subscript[h, t]), 
                      i])(Subscript[δ, (Subscript[x, t]), j])

Equation () may now be compared directly with equation () term for term.

We can verify that the total number of parameters is the same in each case. In equation () there are two conditional probability matrices to parameterise. Q((Subscript[h, t + 1]) | (Subscript[h, t])) requires (Subscript[N, h]) ((Subscript[N, h]) - 1) parameters ((Subscript[N, h]) separate conditional probability, each requiring (Subscript[N, h]) - 1 parameters). Similarly, Q((Subscript[x, t]) | (Subscript[h, t])) requires (Subscript[N, h]) ((Subscript[N, x]) - 1) parameters. Thus the total number of parameters is (Subscript[N, h]^2) + (Subscript[N, h]) (Subscript[N, x]) - 2  (Subscript[N, h]). On the other hand, in equation () there are two matrices of Gibbs distribution parameters (Subscript[λ, i, j]^((hh))), 0 ≤ i ≤ (Subscript[N, h]) - 1, 1 ≤ j ≤ (Subscript[N, h]) - 1 and (Subscript[λ, i, j]^((hx))) 0≤i≤(Subscript[N, h]) - 1, 1 ≤ j ≤ (Subscript[N, x]) - 1. The number of parameters in the Gibbs distribution and the conditional probability formalisms is the same. However, in general, the interpretation of the parameters is not in one to one correspondenceFOOTNOTE.

There remains the problem of reconciling the -log[Q((Subscript[h, 1]))] term in equation () with the (∑_i)(Subscript[λ, 0, i]^((hh)))(Subscript[δ, (Subscript[h, 1]), i]) term in equation (). We are free to choose Q((Subscript[h, 1])) at will as an initial condition in the conditional probability formalism, but (Subscript[λ, 0, i]^((hh))) is already determined as a special case of the (Subscript[λ, i, j]^((hh)))-matrix by the form of the Hanimersley-Clifford expansion. This is not surprising since we assume time shift invariance in the Hammersley-Clifford expansion, but we explicitly violate this invariance in the conditional probability approach by treating t = 1 as a special case (i.e. an initial condition).

For some applications (such as radar clutter modelling) one should preserve explicit time shift invariance, whilst in other applications (such as speech modelling) one might wish to impose explicit initial and final conditions. I shall restrict my attention to the initial condition (Subscript[h, 1]) = 0, thus Q((Subscript[h, 1])) = (Subscript[δ, (Subscript[h, 1]), 0]). The -log[Q((Subscript[h, 1]))] term in equation () and the (∑_i)(Subscript[λ, 0, i]^((hh)))(Subscript[δ, (Subscript[h, 1]), i]) term in equation () may thus be ornitted provided that we clamp (Subscript[h, 1]) = 0 wherever else it appearsFOOTNOTE.

Training a hidden Markov model as a Gibbs distribution

In this section 1 shall derive the detail of the Gibbs machine training algorithm as applied to a hidden Maxkov model. The training scheme turns out to have a similar structure to the Baum-Welch re-estimation scheme [2].

. Summing over histories

I shall henceforth use the expression for λ.U(x,h) given in equation (). Usually one needs to run a Metropolis algorithm in order to estimate the averages in equation (), which leads to long training times for the λ-coefficients. However, in the case of an hidden Markov model this may be circumvented by performing an explicit sum over states of the Markov chain [11]. 1 shall call this a "sum over histories"FOOTNOTE.

In order to develop the sum over histories I shall introduce some convenient notation. Define two matrices (Subscript[S, (Subscript[h, t]), (Subscript[h, t + 1])]^((hh))) and (Subscript[S, (Subscript[h, t]), (Subscript[x, t])]^((hx))) as follows

(Subscript[S, (Subscript[h, t]), (Subscript[h, t + 1])]^((hh)))≡(e^(-(Subscript[λ, (Subscript[h, t]), (Subscript[h, t + 1])]^((hh)))(1 - (Subscript[δ, (Subscript[h, t + 1]), 0]))))
(Subscript[S, (Subscript[h, t]), (Subscript[x, t])]^((hx)))≡(e^(-(Subscript[λ, (Subscript[h, t]), (Subscript[x, t])]^((hx)))(1 - (Subscript[δ, (Subscript[x, t]), 0]))))

Now define forwards and backwards matrices (S^((+))) and (S^((-))) respectively as follows

(Subscript[S, (Subscript[h, t]), (Subscript[x, t]), (Subscript[h, t + 1])]^((+)))≡(Subscript[S, (Subscript[h, t]), (Subscript[x, t])]^((hx)))(Subscript[S, (Subscript[h, t]), (Subscript[h, t + 1])]^((hh)))
(Subscript[S, (Subscript[h, t]), (Subscript[h, t + 1]), (Subscript[x, t + 1])]^((-)))≡(Subscript[S, (Subscript[h, t + 1]), (Subscript[x, t + 1])]^((hx)))(Subscript[S, (Subscript[h, t]), (Subscript[h, t + 1])]^((hh)))

The forwards and backwards pieces of a hidden Markov model.

In Figure (a) and Figure (b) I show diagrammatic representations of (S^((+))) and (S^((-))) respectively. Note how Figure (a) and Figure (b) are two different ways of uniting Figure (a) with Figure (b).

With these definitions I can now write Q(x,h) as

Q(x,h)=(1/Z)\({∏\+\(t = 1\)\%\(τ - 1\)S\_\(h\_t, x\_t, h\_\(t + 1\)\)\%\((+)\)}\)(Subscript[S, (Subscript[h, τ]), (Subscript[x, τ])]^((hx)))\({∏\+\(t = τ\)\%\(T - 1\)S\_\(h\_t, h\_\(t + 1\), x\_\(t + 1\)\)\%\((-)\)}\)

In equation () 1 have split the product of contributions at t = τ, and used forwards matrices for 1 ≤ t < τ and backwards matrices for τ ≤ t ≤ T - 1. There is also an extra term at t = τ. which completes the product. In Figure  I show the decomposition of equation () in diagrammatic form. It is clear how the pieces fit together to build up the whole hidden Markov model chain of Figure , less the 1-point term corresponding to the initial condition which I have discarded.

A forwards/backwards decomposed hidden Markov model.

The marginal probability Q(x) may now be expressed as

Q(x)≡(∑_((Subscript[h, 1]), (Subscript[h, 2]), ⋯, 
                (Subscript[h, T]) = 0)^((Subscript[N, h]) - 1))Q(x,h)=(1/Z)(∑_((Subscript[h, 1]), (Subscript[h, 2]), ⋯, 
                (Subscript[h, T]) = 0)^((Subscript[N, h]) - 1))\({∏\+\(t = 1\)\%\(T - 1\)S\_\(h\_t, x\_t, h\_\(t + 1\)\)\%\((+)\)}\)(Subscript[S, (Subscript[h, T]), (Subscript[x, T])]^((hx)))

An equivalent expression using (S^((-))) can be constructed. Z may be evaluated by summing the numerator over x (note that (∑_x)Q(x)=1 is guaranteed), thus

Z≡(∑_((Subscript[x, 1]), (Subscript[x, 2]), ⋯, 
              (Subscript[x, T]) = 0)^((Subscript[N, x]) - 1))Z(x)
Z(x)≡(∑_((Subscript[h, 1]), (Subscript[h, 2]), ⋯, 
              (Subscript[h, T]) = 0)^((Subscript[N, h]) - 1))\({∏\+\(t = 1\)\%\(T - 1\)S\_\(h\_t, x\_t, h\_\(t + 1\)\)\%\((+)\)}\)(Subscript[S, (Subscript[h, T]), (Subscript[x, T])]^((hx)))

. Maximising G

1 shall now use the U(x,h) defined in equation () to evaluate ∂(G/∂)(Subscript[λ, i]), using the Gibbs machine training algorithm in equation (). The (Subscript[λ, i]) are the various λ-coefficients that appear in equation (). Firstly define two matrices (Subscript[D, i, j]^((hh))) and (Subscript[D, i, j]^((hx))) as

(Subscript[D, i, j]^((hh)))(t)≡\(\((1 - δ\_\(j, 0\))\) \(δ\_\(h\_t, i\)\) δ\_\(h\_\(t + 1\), j\)\)
(Subscript[D, i, j]^((hx)))(t)≡\(\((1 - δ\_\(j, 0\))\) \(δ\_\(h\_t, i\)\) δ\_\(x\_t, j\)\)

Using these definitions I obtain two types of derivative

((∂G)/(∂(Subscript[λ, i, j]^((hh)))))=(∑_(t = 1)^(T - 1))[(Subscript[⟨(Subscript[D, i, j]^((hh)))(t)⟩, free])-(Subscript[⟨(Subscript[D, i, j]^((hh)))(t)⟩, clamped])]
((∂G)/(∂(Subscript[λ, i, j]^((hx)))))=(∑_(t = 1)^T)[(Subscript[⟨(Subscript[D, i, j]^((hx)))(t)⟩, free])-(Subscript[⟨(Subscript[D, i, j]^((hx)))(t)⟩, clamped])]

where free denotes averaging with respect to Q(x,h), and clamped denotes averaging with respect to Q(h|x)P(x), using the same notation as defined in equation ().

Each of the averages in equation () is a linear combination of average state occupations. These are analogous to histograms or frequencies of occurrence. For instance, ⟨(Subscript[δ, (Subscript[h, t]), 
              i]) (Subscript[δ, (Subscript[x, t]), j])⟩ averages (over Q(x,h)) a quantity that is 1 when (Subscript[h, t]) = i and (Subscript[x, t]) = j, and zero otherwise. This average is the answer to the question "What proportion of the time is ((Subscript[h, t]), (Subscript[x, t])) to be found in state (i, j) when the Gibbs distribution is allowed to explore its states without external influence?". There are analogous interpretations of each of the other averages in equation ().

The free and clamped averages are explicitly

(Subscript[\(⟨D\_\(i, j\)\%\((·)\)⟩\), free])≡(∑_((Subscript[x, 1]), (Subscript[x, 2]), ⋯, 
              (Subscript[x, T]) = 0)^((Subscript[N, x]) - 1))(∑_((Subscript[h, 1]), (Subscript[h, 2]), ⋯, 
                (Subscript[h, T]) = 0)^((Subscript[N, h]) - 1))Q(x,h)(Subscript[D, i, j]^((·)))
(Subscript[\(⟨D\_\(i, j\)\%\((·)\)⟩\), clamped])≡(∑_(x∈training set))(((∑_((Subscript[h, 1]), (Subscript[h, 2]), ⋯, 
                  (Subscript[h, T]) = 0)^((Subscript[N, h]) - 1))Q(x,h)(Subscript[D, i, j]^((·))))/(Q(x)))

The gradients ((∂G)/(∂(Subscript[λ, i, j]^((hh))))) and ((∂G)/(∂(Subscript[λ, i, j]^((hx))))) can thus be evaluated by measuring the various state occupations that appear in equation () to yield

((∂G)/(∂(Subscript[λ, i, j]^((hh)))))=(1/Z)(∑_(τ = 1)^(T - 1))(∑_((Subscript[x, 1]), (Subscript[x, 2]), ⋯, 
                    (Subscript[x, T]) = 0)^((Subscript[N, x]) - 1))(∑_((Subscript[h, 1]), (Subscript[h, 2]), ⋯, 
                      (Subscript[h, T]) = 0)^((Subscript[N, h]) - 1))(Subscript[D, i, j]^((hh)))(τ)\({∏\+\(t = 1\)\%\(τ - 1\)S\_\(h\_t, x\_t, h\_\(t + 1\)\)\%\((+)\)}\)(Subscript[S, (Subscript[h, τ]), (Subscript[x, τ])]^((hx)))\({∏\+\(t = τ\)\%\(T - 1\)S\_\(h\_t, h\_\(t + 1\), x\_\(t + 1\)\)\%\((-)\)}\)-(1/N)(∑_(τ = 1)^(T - 1))(∑_(x∈training set))(1/(Z(x)))(∑_((Subscript[h, 1]), (Subscript[h, 2]), ⋯, 
                        (Subscript[h, T]) = 0)^((Subscript[N, h]) - 1))(Subscript[D, i, j]^((hh)))(τ)\({∏\+\(t = 1\)\%\(τ - 1\)S\_\(h\_t, x\_t, h\_\(t + 1\)\)\%\((+)\)}\)(Subscript[S, (Subscript[h, τ]), (Subscript[x, τ])]^((hx)))\({∏\+\(t = τ\)\%\(T - 1\)S\_\(h\_t, h\_\(t + 1\), x\_\(t + 1\)\)\%\((-)\)}\)
((∂G)/(∂(Subscript[λ, i, j]^((hx)))))=ditto with (Subscript[D, i, j]^((hh)))(τ)⟶(Subscript[D, i, j]^((hx)))(τ) and T-1⟶T

where Z and Z(x) are defined in equation (). The free average is self explanatory. The clamped average depends on a sum over the training set, from which x should be selected at random in order to emulate the probability P(x).

The value of G itself can also be calculated up to an additive constant. Thus

G=-log[Z]+(1/N)(∑_(x∈training set))log[Z(x)]+constant

Equation () can be used to verify that a gradient ascent algorithm based on equation () is indeed increasing the value of GFOOTNOTE.

Algorithmic details

In this section I shall explain how a general purpose hidden Markov model training algorithm can be written, which allows both Gibbs machine and Baum-Welch training algorithms to be run.

. Clamped matrix products

The hidden Markov model training algorithms of interest all reduce to evaluating products of matrices. For instance, Baum-Welch re-estimation involves the recursive computation of so-called α and β matrices [2]. The Gibbs machine training algorithm reduces to equation (), where the "clamped" case is essentially the same as Baum-Welch, and the "free" case involves a complete sum over histories. It would be prudent to write an algorithm that allows products of (S^((±))) matrices to be evaluated whilst clamping arbitrary subsets of indices that would normally be summed over. Everything that I require may be evaluated as one or another special case of this general algorithm.

The matrix products that I need to evaluate are

(Subscript[M, (Subscript[h, t]), (Subscript[x, t])]^((+))) (t)=(Subscript[C, (Subscript[h, t]), (Subscript[x, t])]^((hx)))(t)(Subscript[S, (Subscript[h, t]), (Subscript[x, t])]^((hx)))(∑_((Subscript[h, t]) = 0)^((Subscript[N, h]) - 1))(∑_((Subscript[x, t]) = 0)^((Subscript[N, x]) - 1))(Subscript[C, (Subscript[h, t - 1]), (Subscript[h, t])]^((hh)))(t - 1)(Subscript[S, (Subscript[h, t - 1]), (Subscript[h, t])]^((hh)))(Subscript[M, (Subscript[h, t - 1]), (Subscript[x, t - 1])]^((+)))(
                  t - 1)
(Subscript[M, (Subscript[h, t]), (Subscript[x, t])]^((-))) (t)=(Subscript[C, (Subscript[h, t]), (Subscript[x, t])]^((hx)))(t)(Subscript[S, (Subscript[h, t]), (Subscript[x, t])]^((hx)))(∑_((Subscript[h, t]) = 0)^((Subscript[N, h]) - 1))(∑_((Subscript[x, t]) = 0)^((Subscript[N, x]) - 1))(Subscript[C, (Subscript[h, t]), (Subscript[h, t + 1])]^((hh)))(t)(Subscript[S, (Subscript[h, t]), (Subscript[h, t + 1])]^((hh)))(Subscript[M, (Subscript[h, t + 1]), (Subscript[x, t + 1])]^((-)))(
                  t + 1)

where I use the boundary conditions

(Subscript[M, (Subscript[h, 1]), (Subscript[x, 1])]^((+)))(
        t = 1) = (Subscript[δ, (Subscript[h, 1]), 
            0]) (Subscript[δ, (Subscript[x, 1]), 0])
(Subscript[M, (Subscript[h, 1]), (Subscript[x, 1])]^((-)))(
        t = T) = (Subscript[δ, (Subscript[h, T]), 
            (Subscript[N, h - 1])]) (Subscript[δ, (Subscript[x, 1]), 0])

The matrices (Subscript[C, (Subscript[h, t]), (Subscript[x, t])]^((hx)))(t), 2 ≤ t ≤ T - 1 and (Subscript[C, (Subscript[h, t]), (Subscript[h, t + 1])]^((hh)))(t), 1 ≤ t ≤ T - 1 specify the clamps as follows

(Subscript[C, (Subscript[h, t]), (Subscript[x, t])]^((hx)))(t) or (Subscript[C, (Subscript[h, t]), (Subscript[h, t + 1])]^((hh)))(t)={({{1, permitted by clamp}, {0, forbidden by clamp}})

In practice I encode the boundary conditions of equation () in the clamps. The clamp representation in equation () is both very general and mathematically convenient. However, it is highly redundant for the usual "all or one" type of clamp that I need for my current purposesFOOTNOTE. In such cases an appropriate coding of the matrix should be used to save memory.

. Evaluating the training gradients

I shall now express equation () in terms of the (M^((±))) matrices. ((∂G)/(∂(Subscript[λ, i, j]^((hh))))) becomes

((∂G)/(∂(Subscript[λ, i, j]^((hh)))))=(((1 - (Subscript[δ, j, 0])))/Z)(∑_(τ = 1)^(T - 1))(∑_((Subscript[x, τ]), 
                    (Subscript[x, τ + 1]) = 0)^((Subscript[N, x]) - 1))(Subscript[M, i, (Subscript[x, τ])]^((+)))(τ)(Subscript[S, i, j]^((hh)))(Subscript[M, j, 
                        (Subscript[x, τ + 1])]^((-)))(τ + 1)-(((1 - (Subscript[δ, j, 0])))/N)(∑_(τ = 1)^(T - 1))(∑_(x∈training set))(1/(Z(x)))(Subscript[M, i, (Subscript[x, τ])]^((+)))(τ)(Subscript[S, i, j]^((hh)))(Subscript[M, j, 
                        (Subscript[x, τ + 1])]^((-)))(τ + 
                      1)

Similarly, ((∂G)/(∂(Subscript[λ, i, j]^((hx))))) becomes

((∂G)/(∂(Subscript[λ, i, j]^((hx)))))=(((1 - (Subscript[δ, j, 0])))/Z)(∑_(τ = 1)^T)(Subscript[M, i, j]^((+)))(τ)(1/(Subscript[S, i, j]^((hx))))(Subscript[M, i, 
                      j]^((-)))(τ)-(((1 - (Subscript[δ, j, 0])))/N)(∑_(τ = 1)^T)(∑_(x∈training set))(1/(Z(x)))(Subscript[M, i, j]^((+)))(τ)(1/(Subscript[S, i, j]^((hx))))(Subscript[M, i, 
                        j]^((-)))(τ)

In the second term of equation () and equation () the (M^((±))) matrices should be appropriately clamped by the training set xFOOTNOTE. Z and Z(x) may be expressed as

Z=({{(Subscript[M, (Subscript[N, h]) - 1, 0]^((+)))(T),  , boundaries clamped}})
Z(x)=({{(Subscript[M, (Subscript[N, h]) - 1, 0]^((+)))(T),  , both x and boundaries clamped}})

. Imposing the Hammersley-Clifford constraints

In the discussion following the Harnmersley-Clifford expansion in equation () I showed how the Hammersley-Clifford expansion places the correct number of constraints on the total number of independent parameters. These constraints become the (1 - (Subscript[δ, j, 0])) factors in equation () and equation (). Thus ((∂G)/(∂(Subscript[λ, i, 0]^((hh))))) and ((∂G)/(∂(Subscript[λ, i, 0]^((hx))))) are zero because there are no such λ-coefficients in equation (). Thus the zeroth column of each λ-matrix is constrained. This is aesthetically ugly, because it selects a specific column as the "odd one out"FOOTNOTE.

It proves to be mathematically convenient to simply ignore the (1 - (Subscript[δ, j, 0])) factors in equation () and equation (), and to impose a row normalisation constraint. Thus the (S^((hh))) and (S^((hx))) matrices in equation () should be redefined as

(Subscript[S, i, j]^((hh)))≡(e^(-(Subscript[λ, i, j]^((hh))))), (∑_(j = 0)^((Subscript[N, h]) - 1))(Subscript[S, i, j]^((hh)))=1
(Subscript[S, i, j]^((hx)))≡(e^(-(Subscript[λ, i, j]^((hx))))), (∑_(j = 0)^((Subscript[N, x]) - 1))(Subscript[S, i, j]^((hx)))=1

This row normalisation procedure turns (S^((hh))) and (S^((hx))) into stochastic matrices.

I can now compare the Gibbs machine training algorithm directly with the BaumWelch re-estimation method. Whereas the Gibbs machine training algorithm makes small adjustments to the model based on the gradient of G, the Baum-Welch re-estimation method makes potentially large adjustments to the model to maximise (a lower bound for) G at each stage. Furthermore, given the row normalisation prescription for the Hammersley-Clifford expansion of the hidden Markov model of §4.1, the model is parameterised in the same way for both types of training methodFOOTNOTE. At this point it is clear that for the hidden Markov model structure of §. the Gibbs machine approach offers nothing that Baum-Welch reestimation cannot do better.

My pedagogical approach to comparing the Gibbs machine training algorithm with the Baum-Welch re-estimation method has made clear the intimate connection between the two methods. In the appendix I show how this connection can be explored much further.

Conclusions

I have demonstrated how my Gibbs machine training algorithm [9] is intimately related to the Baum-Welch re-estimation method for hidden Markov models [2]. In the appendix I present a proof of this intimate relationship in the general case, and I point out how one can combine the two methods in a hybrid training scheme.

The principal disadvantage of Gibbs machine training is that it is a gradient ascent method, so it tends to be slower to converge than Baum-Welch re-estimation. Furthermore, the Gibbs machine can require extensive Monte Carlo simulations - although not for the hidden Markov model that I have studied.

The principal advantage of using a Gibbs machine is the Hammersley-Clifford expansion which allows one systematically to parameterise Markov models starting with 1-point terms, progressing through 2-point terms, until finally a complete parameterisation is recovered. This allows one to construct "cheap and cheerful" Markov models which model some aspects of nth order correlations, without actually introducing functions of n variables into one's model.

A secondary advantage of the Gibbs machine training algorithm is its ability to fill in where the Baum-Welch method fails to work - such as when the Markov model has loops. This applies independently of whether one is using a Gibbs machine in its "cheap and cheerful" mode or not.

In summary, my Gibbs machine training algorithm [9] and the Baurn-Welch re-estimation method [2] are two sides of the same coin. One should use both methods in tandem to train a general Markov model - the Baum-Welch re-estimation method trains the tree-like periphery of the model, whilst the Gibbs machine training algorithm deals with the rest (see the appendix for more details).

Finally, this work would not have occurred if I had not been irnmersed in an environment consisting of (lots of ) hidden Markovians and (one or two) Bayesians. In particular, I thank John Bridle for saying "Baum-Welch" often enough that I eventually read [2, 8], Martin Russell for convincing me that my Gibbs machine might have something useful to offer to hidden Markov modellers, and Andrew Varga for bringing [11] to my attention.

Appendix A The Baum-Welch re-estimation method

In this section I shall rederive the Baum-Welch re-estimation method [2] using the minimum possible set of assumptions, and I then deduce the re-estimation formulae for a mixture a Gaussians [8]. I should point out that there is no essential difference between the discrete and the continuum versions of the re-estimation methodFOOTNOTE. I should record here that the EM algorithm [6] is essentially the same as the Baum-Welch method.

A. The fundamental derivation

For compatibility with my other derivations, I shall maximise a relative entropy G, estimated using the whole training set, rather than a likelihood that is measured for a single training set vector [2] - there is no essential difference between these two optimisation criteria. Thus define G(λ) as

G(λ)≡-∫dx P(x)log[((P(x))/(Q(x|λ)))]≤0

where I now explicitly write the λ dependence of the model Q(x|λ)FOOTNOTE. It irnmediately follows that

ΔG(λ',λ)≡G(λ')-G(λ)=∫dx P(x)log[((Q(x|λ'))/(Q(x|λ)))]

where λ' is the sought-after re-estimate of λ that maximally increases the value of G(λ).

The problem of maximising ΔG(λ',λ) is non-trivial because it has an integration (over hidden variables) inside the logarithm function. The solution to this impasse is to find a lower bound for ΔG(λ',λ) that uses Jensen's inequality (for convex functions) to move the integral outside the logarithm function. Thus

log[((Q(x|λ'))/(Q(x|λ)))]=log[((∫dh Q(x,h|λ'))/(Q(x|λ)))]
=log[∫dh ((Q(x,h|λ))/(Q(x|λ)))((Q(x,h|λ'))/(Q(x,h|λ)))]
=log[∫dh Q(h|x,λ)((Q(x,h|λ'))/(Q(x,h|λ)))]
≥∫dh Q(h|x,λ)log[((Q(x,h|λ'))/(Q(x,h|λ)))]

The 4 steps of this derivation have the following interpretation:

• Step 1. Write out the integration over hidden variables h explicitly.

• Steps 2 and 3. You must now make sure that the integration has a measure with unit total mass (i.e. normalised), because Jensen's inequality needs this. You could use either Q(h|x,λ) or Q(h|x,λ'). However, Jensen's inequality will move this measure outside the logarithm, so it is better to use Q(h|x,λ) for which we already have available a means of performing the average (i.e. clamping).

• Step 4. Apply Jensen's inequality to the convex logarithm function. Equality holds if and only if Q(h|x,λ')=Q(h|x,λ) ∀h.

This completes the derivation of the lower bound, so ΔG(λ',λ) now satisfies the inequality

ΔG(λ',λ)≥∫dx P(x)∫dh Q(h|x,λ)log[((Q(x,h|λ'))/(Q(x,h|λ)))]

In equation () there is no integration over h to be found inside the logarithm function, which makes its evaluation relatively simple, as intended.

At this point it is convenient to depart from the derivation in [2], by assuming that Q(x,h|λ') is a Gibbs distribution. This is not actually a restriction, since the Hammersley-Clifford expansion can be used to construct any model that we require. Naturally, the Gibbs distribution formalism includes the hidden Markov model discussed in [2] as a special case. Thus Q(x,h|λ') becomes

Q(x,h|λ')≡(1/(Z(λ')))(e^(-λ'.U(x,h)))
Z(λ')≡∫dx ∫dh (e^(-λ'.U(x,h)))

Insert equation () into equation () to obtainFOOTNOTE

ΔG(λ',λ)≥∫dx P(x)∫dh Q(h|x,λ) [-log[Z(λ')]-λ'.U(x,h)]+constant
=-log[Z(λ')]-λ'.(Subscript[⟨U⟩, clamped])+constant

where (Subscript[\(⟨⋯⟩\), clamped]) denotes an average which uses P(x) to sample from the training set, and then uses Q(h|x,λ) to determine the integral over h.

In equation () the log[Z(λ')] depends on λ' in a complicated fashion (see the definition of Z(λ') in equation ()). This prevents easy computation of the λ' that maximises this lower bound for ΔG(λ',λ).

However, there is a special case where Z(λ') is constant. Thus if Q(x,h|λ') can be expressed entirely as a product of conditional probabilities, then Z(λ') is fixed at 1, by construction. This type of Q(x,h|λ') has a tree-like structure, with a conditional probability factor associated with each branch. The expression for the lower bound in equation () thus simplifies to

ΔG(λ',λ)≥-λ'.(Subscript[⟨U⟩, clamped])+constant

Now regroup the ((Subscript[λ, i])^′) and (Subscript[U, i])(x,h) into matrices, each of which corresponds to one of the conditional probability factors in Q(x,h|λ')

λ'.U(x,h)=(∑_(i, j, b))(Subscript[((λ^′)), i, j]^b)(Subscript[U, i, j]^b)(x,h)
\(∑\+j e\^\(-\(λ\^′\)\_\(i, j\)\%b\) = 1\)

where b sums over branches of the tree, and i and j sum over the "parent" and "child" states (respectively) associated with each branch. I have assumed that the (Subscript[U, i, j]^b)(x,h) are a complete set of non-overlapping sampling functions - Dirac deltas in the continuous case, Kronecker deltas in the discrete caseFOOTNOTE. The row normalisation of the conditional probability factors reveals itself in the normalisation of sums of (e^(-(Subscript[((λ^′)), i, j]^b))) factors. I can now write equation () as

ΔG(λ',λ)≥-(∑_(i, j, 
                b))(Subscript[((λ^′)), i, j]^b) 
                (Subscript[n, i, j]^b)+constant

where (Subscript[n, i, j]^b)≡(Subscript[\(⟨U\_\(i, j\)\%b⟩\), clamped]) is the relative frequency with which state (i, j) on branch b occurs when the Gibbs distribution is clamped.

It is now a simple matter to derive the re-estimate of (Subscript[((λ^′)), i, j]^b) that maximises ΔG(λ',λ). Thus introduce a Lagrange multiplier term for each row-normalisation constraint, and differentiate with respect to (Subscript[((λ^′)), i, j]^b)

\(∂\/∂\(λ\^′\)\_\(i, j\)\%b\)[-(∑_(i, j, 
                  b))(Subscript[((λ^′)), i, j]^b) 
                  (Subscript[n, i, j]^b)-(∑_(i, b))(Subscript[μ, i]^b)(∑_j)\(e\^\(-\(λ\^′\)\_\(i, j\)\%b\)\)]=-(Subscript[n, i, j]^b) + (Subscript[μ, i]^b) 
            (e^(-(Subscript[((λ^′)), i, j]^b)))

This must be zero when \(TraditionalForm\`\(λ\^′\)\_\(i, j\)\%b = λ\_\(i, j\)\%b\). Sum over j and use the normalisation condition to obtain (Subscript[μ, i]^b) = (∑_j) (Subscript[n, i, j]^b), whence

(Subscript[((λ^′)), i, j]^b) = 
      ((Subscript[n, i, j]^b)/((∑_j) (Subscript[n, i, j]^b)))

This is the Baum-Welch re-estimation scheme [2]. 1 have derived it assuming only that the model Q(x,h|λ) is tree-like, which can be expressed as a product of conditional probabilities with a fixed total probability mass. My derivation is similar to, but much (simpler than, that used in [6], where the (now fashionable) EM algorithm for re-estimating model parameters is derivedFOOTNOTE.

A. Application to Gaussian mixtures

A special case of a Gibbs distribution with hidden variables, Q(x,h|λ), is a Gaussian whose parameters are conditioned on the state of a hidden variable. Thus define

Q(x|m(s),C(s))ρ(
            s)≡ρ(s)(1/(Sqrt[det(2π C(s))]))(e^(-(1/2)(((x-m(s)))^T)(C^(-1))(s)(x-m(s))))

where s is a hidden variable, ρ(s) is its probability density function (ρ(s) is normalised as ∫ds ρ(s)=1), m(s) is a "mean vector", C(s) is a syrnrnetric "covariance matrix", and Q(x|m(s),C(s)) is the probability density function of x conditioned on s.

I can relate equation () to equation () as follows

λ'.U=-(1/2)log det(2π C'(s))+log((ρ^′)(s))-(1/2)(((x-m'(s)))^T)((C')^(-1))(s)(x-m'(s))
Z(λ')=1

Z = 1 because Q(x|m(s),C(s))ρ(
            s) is normalised (by construction), so I can use the simplified form of equation () given in equation (). This permits me to use the standard Baum-Welch method to re-estimate ρ(s), m(s) and C(s) [8]FOOTNOTE. In this case the clamped average is computed as follows (using the old unprimed parameters)

(Subscript[\(⟨\([\)\(⋯\)\(]\)⟩\), clamped])≡∫dx ds Q(s|x)P(x) [⋯]

where [⋯] can be replaced by any sampling function that one wishes to average.

In order to re-estimate the model parameters I need to differentiate λ'.U with respect to each parameter of interest, and perform a clamped average using the original parameter values. There are three cases to consider:

• Re-estimate ρ(s). The steps in the derivation are a standard application of Lagrange multiplier methods, as follows:

({{λ'.U, λ'.U, ⟶, log((ρ^′)(s))}, {Average, (Subscript[⟨λ'.U⟩, clamped]), =, log((ρ^′)(s))∫dx Q(s|x)P(x)}, {Lagrange, ℒ, =, -μ∫ds (ρ^′)(s)}, {Differentiate, (∂/(∂(ρ^′)(s)))[(Subscript[⟨λ'.U⟩, clamped])+ℒ], =, (1/((ρ^′)(s)))∫dx Q(s|x)P(x)-μ}, {μ, μ, =, ∫ds∫dx Q(s|x)P(x)=1}, {Re - estimate, (OverHat[ρ])(s), =, ∫dx Q(s|x)P(x)}})

• Re-estimate m'(s) as follows:

({{λ'.U, λ'.U, ⟶, (x^T)((C')^(-1))(s)m'(s)-(1/2)((m'(s))^T)((C')^(-1))(s)m'(s)}, {Differentiate, ((∂λ'.U)/(∂m'(s))), =, ((C')^(-1))(s)[x-m'(s)]}, {Average, (Subscript[⟨((∂λ'.U)/(∂m'(s)))⟩, clamped]), =, ((C')^(-1))(s)∫dx Q(s|x)P(x)x-((C')^(-1))(s)m'(s)∫dx Q(s|x)P(x)}, {Re - estimate, (OverHat[m])(s), =, ((∫dx Q(s|x)P(x)x)/(∫dx Q(s|x)P(x)))}})

• Re-estimate m'(s) as follows:

({{λ'.U, λ'.U, ⟶, (1/2)log det((1/(2  π)) ((C')^(-1))(s))-(1/2)(((x-m'(s)))^T)((C')^(-1))(s)(x-m'(s))}, {Differentiate, ((∂λ'.U)/(∂((C')^(-1))(s))), =, (1/2)C'(s)-(1/2)(x-m'(s))(((x-m'(s)))^T)}, {Average, (Subscript[⟨((∂λ'.U)/(∂((C')^(-1))(s)))⟩, clamped]), =, (1/2)C'(s)∫dx Q(s|x)P(x)-(1/2)∫dx Q(s|x)P(x)(x-m'(s))(((x-m'(s)))^T)}, {Re - estimate, (OverHat[C])(s), =, ((∫dx Q(s|x)P(x)(x-(OverHat[m])(s))(((x-(OverHat[m])(s)))^T))/(∫dx Q(s|x)P(x)))}})

In order to obtain the result for (OverHat[C])(s) I use the identities

log det ≡ trace log
\(TraditionalForm\`\(e\^A\) 
        e\^B ≡ e\^\(A + B + commutator\ terms\)\)
log[A B] ≡ log[A] + 
        log[B] + commutator terms
trace log[A B] ≡ trace(
        log[A] + log[B])

The three re-estimation formulae above can be used to maximise G when the model Q(x,h|λ) has the form shown in equation () [8]. Note that (OverHat[m])(s) must be computed before (OverHat[C])(s), because (OverHat[C])(s) depends on (OverHat[m])(s).

A. Hybrid re-estimation method

Recall the fundamental inequality in equation (), and suppose that Q(x,h|λ') factorises as follows

Q(x,h|λ')=Q((Subscript[x, 2]),(Subscript[h, 2])|(Subscript[x, 1]),(Subscript[h, 1]),(Subscript[λ, 2])')Q((Subscript[x, 1]),(Subscript[h, 1])|(Subscript[λ, 1])')

where x≡((Subscript[x, 1]),(Subscript[x, 2])), h≡((Subscript[h, 1]),(Subscript[h, 2])), and λ≡((Subscript[λ, 1]),(Subscript[λ, 2])). The Q((Subscript[x, 1]),(Subscript[h, 1])|(Subscript[λ, 1])') factor models the properties of some of the variables ((Subscript[x, 1]),(Subscript[h, 1])), whilst the Q((Subscript[x, 2]),(Subscript[h, 2])|(Subscript[x, 1]),(Subscript[h, 1]),(Subscript[λ, 2])') factor models the remaining variables conditioned on ((Subscript[x, 1]),(Subscript[h, 1])).

The essential step is to assume that the Q((Subscript[x, 2]),(Subscript[h, 2])|(Subscript[x, 1]),(Subscript[h, 1]),(Subscript[λ, 2])') factor is tree-like, whereas the Q((Subscript[x, 1]),(Subscript[h, 1])|(Subscript[λ, 1])') factor is potentially fully connectedFOOTNOTE. In a sense, one performs a decomposition of (x,h) into "internal" variables ((Subscript[x, 1]),(Subscript[h, 1])) and "peripheral" variables ((Subscript[x, 2]),(Subscript[h, 2])). ΔG(λ',λ) thus satisfies the inequality

ΔG(λ',λ)≥∫dx P(x)∫dh Q(h|x,λ)[log[Q((Subscript[x, 2]),(Subscript[h, 2])|(Subscript[x, 1]),(Subscript[h, 1]),(Subscript[λ, 2])')]+log[Q((Subscript[x, 1]),(Subscript[h, 1])|(Subscript[λ, 1])')]]+constant

• The Baum-Welch piece. The Q((Subscript[x, 2]),(Subscript[h, 2])|(Subscript[x, 1]),(Subscript[h, 1]),(Subscript[λ, 2])') has Z = 1 (by construction), so this term reduces to the case in equation (). The clamped average must sum over all permitted states, which includes ((Subscript[x, 1]),(Subscript[h, 1])), so the re-estimate of the peripheral piece of the model depends on a probabilistic "boundary condition" supplied by the internal piece of the model.

• The Gibbs distribution piece. The Q((Subscript[x, 1]),(Subscript[h, 1])|(Subscript[λ, 1])') does not have a simple Z factor, so it is of the class shown in equation (), which does not yield to Baum-Welch reestimation. Because this term must be improved by gradient ascent, I shall replace λ' by λ in the differentiations. Introduce a Gibbs distribution model

Q((Subscript[x, 1]),(Subscript[h, 1])|(Subscript[λ, 1]))≡(1/(Z((Subscript[λ, 1]))))(e^(-(Subscript[λ, 1]).U((Subscript[x, 1]),(Subscript[h, 1]))))
Z((Subscript[λ, 1]))≡∫(Subscript[dx, 1]) ∫(Subscript[dh, 1]) (e^(-(Subscript[λ, 1]).U((Subscript[x, 1]),(Subscript[h, 1]))))

Differentiate ΔG(λ',λ) with respect to (Subscript[λ, 1])' to obtain

(Subscript[((∂ΔG(λ',λ))/(∂(Subscript[λ, 1])'))|, λ'=λ])=(Subscript[⟨U⟩, free])-(Subscript[⟨U⟩, clamped])

where.the "free" average is over Q((Subscript[x, 1]),(Subscript[h, 1])|(Subscript[λ, 1])), and the "clamped" average is over Q(h|x,λ)P(x). This "clamped" average and the probabilistic boundary condition used in the Baum-Welch re-estirnation involve the same calculation.

References

[1] Ackley D H, Hinton G E and Sejnowski T J, 1985, Cogn. Sci., 9, 147-169, A learning algorithm for Boltzmann machines.

[2] Baum L E, 1972, Inequalities, 3, 1-8, An inequality and associated maxirnisation technique in statistical estimation of probabilistic functions of a Markov process.

[3] Besag J E, 1974, J. Roy. Statist. Soc. Ser. B, 36, 192-236, Spatial interaction and the statistical analysis of lattice systems.

[4] Bridle J S and Moore R K, 1984, Proc. Inst. Acoust. Autumn Conf., Windermere, 315-322, Boltzmann machines for speech pattern processing.

[5] Cox R P, 1946, Am. J. Phys., 17, 1-13, Probability, frequency and reasonable expectation.

[6] Dempster A P, Laird N M and Rubin D B, 1977, J. Roy. Statist. Soc. Ser. B, 39, 1-37, Maximum likelihood from incomplete data via the EM algorithm.

[7] Jeffries H, 1939, Theory of probability, Oxford University Press.

[8] Liporace L A, 1982, IEEE Trans. Inform. Theory, 28(5), 729-734, Maximum likelihood estimation for multivariate observations of Markov sources.

OPEN [9] Luttrell S P, 1985, RSRE Memo., 3815, The implications of Boltzmann-type machines for SAR data processing: a preliminary survey.

[10] Metropolis N, Rosenbluth A W, Rosenbluth M N, Teller A H and Teller E, 1953, J. Chem. Phys., 21, 1087-1092, Equation of state calculations by fast computing machines.

[11] Zhao Y, Atlas L and Zhuang X, 1988, Proc. ICASSP'88, 1, 28-31, Application of the Gibbs distribution to hidden Markov modelling in isolated word recognition.

Footnotes

1. A Markov random field model is a particular type of "neural network" model.

2. It is unnecessary, and sometimes rnisleading, to interpret this parameterisation of P(x) in terms of statistical thermodynamics.

3. This is the principle advantage of the Hanunersley-Clifford expansion.

4. Even quantum mechanics is a hidden variables theory in this picture - the state vector itself is the set of hidden variables.

5. The conventional Boltzmann machine is of limited use because it uses binary variables and quadratic interactions to emulate the behaviour of non-binary variables and higher order interactions.

6. This intuitive picture has a rigorous information theoretic derivation [9].

7. Recently the phrase Gibbs sampler has crept into popular use. This is the Metropolis algorithm by another name - there is nothing new to be found in the Gibbs sampler.

8. Strictly speaking, the statistical thermodynamic flavour that has surrounded interpretations of such computer simulations of Q(x,h) should not be carried over to interpretations of the theoretical development of expressions involving Q(x,h).

9. Bridle et al suggest that a conventional (i.e. binary) Boltzmann machine should be used to simulate an hidden Markov model, and that a "continuous flow annealing" scheme should be used. There are two main problems here: a binary Boltzmann machine is inappropriate for a non-binary hidden Markov model, and annealing is computationally horrendously expensive.

10. Zhao et al apparently overlook the connection between the Gibbs distribution and the Boltzmann machine. Despite being aware of the Hanunersley-Clifford expansion, inexplicably they fail to parameterise their Gibbs distribution correctly, leading to an impaired performance when compared with the hidden Markov model approach.

11. The languages of discrete and continuous mathematics are only superficially different.

12. It is a trivial task to transform conditional probability into Gibbs distribution language, but the converse is usually not possible in closed form.

13. I perform this sleight of hand to simplify the mathematics.

14. The sum over histories is familiar to those who have used Feynman's path integral formulation of quantum mechanics, for instance.

15. This is a luxury that is not enjoyed by the original Metropolis algorithm based Boltzmann machine [1].

16. In this case either all elements or only a single element of the clamp matrix have value 1.

17. I know that it is sloppy to write (1/(Subscript[S, i, j]^((hx)))) in the second term of equation ()!

18. Blame the original form of the Hammersley-Clifford expansion in equation () for this.

19. This is not a property of Markov models in general.

20. This remark is prompted by comments that I have heard being made about supposed differences between discrete and continuous variables.

21. I have not assumed that Q(x|λ) is a Gibbs distribution.

22. From this point onwards in this appendix I made many corrections on 16 November 1990 to use primed variables where appropriate.

23. Sometimes these are called indicator functions.

24. The EM algorithm [6] is exactly the same as the Baum-Welch re-estimation algorithm [2]. In fact, in the discussion following the main text of [6], Baum himself records some interesting remarks about who thought of the method first!

25. In a sense, there is nothing in [8] that is not already reported in [2].

25. All networks can be decomposed in this fashion.