# On Maximum a Posteriori Estimation of Hidden Markov Processes

###### Abstract

We present a theoretical analysis of Maximum a Posteriori (MAP) sequence estimation for binary symmetric hidden Markov processes. We reduce the MAP estimation to the energy minimization of an appropriately defined Ising spin model, and focus on the performance of MAP as characterized by its accuracy and the number of solutions corresponding to a typical observed sequence. It is shown that for a finite range of sufficiently low noise levels, the solution is uniquely related to the observed sequence, while the accuracy degrades linearly with increasing the noise strength. For intermediate noise values, the accuracy is nearly noise-independent, but now there are exponentially many solutions to the estimation problem, which is reflected in non–zero ground–state entropy for the Ising model. Finally, for even larger noise intensities, the number of solutions reduces again, but the accuracy is poor. It is shown that these regimes are different thermodynamic phases of the Ising model that are related to each other via first-order phase transitions.

On Maximum a Posteriori Estimation of Hidden Markov Processes

Armen Allahverdyan Yerevan Physics Institute Yerevan 375036, Armenia aarmen@mail.yerphi.am Aram Galstyan USC Information Sciences Institute Marina del Rey, CA 90292, USA galstyan@isi.edu

## 1 Introduction

Hidden Markov Models (HMM) are used extensively for modeling sequential data in various areas [9, 4]: information theory, signal processing, bioinformatics, mathematical economics, linguistics, etc. One of the main problems underlying many applications of HMMs amounts to inferring the hidden state sequence based on noise-corrupted observation sequence . This is often done through maximum a posteriori (MAP) approach, which finds an estimate by maximizing the posterior probability .

The computational solution to the MAP optimization problem is readily available via the Viterbi algorithm [9]. Despite its extensive use in many applications, however, the properties of MAP estimation, and specifically, the structure of its solution space, have received surprisingly little attention. On the other hand, it is clear that choosing a single state sequence might be insufficient for adequately understanding the structure of the inferred process. To get a more complete picture, one needs to know whether there are other nearly optimal sequences, how many of them, how they compare with the optimal solution, and so on.

Generally, the structure of an inference method can be characterized by the accuracy of the estimation, and the number of solutions that the method can produce in response to a given sequence . In this paper we study the structure of MAP inference for the simplest binary, symmetric HMM. As an accuracy measure we employ the moments of the estimated sequence in comparison of those of the actual sequence , while the number of possible estimates will be characterized by its averaged logarithm . The binary symmetric HMM is studied by reducing it to the Ising model in random fields, a relation well-known both in computer science [6] and statistical physics [8]. In this way, the average cost of MAP and the logarithm of the number of solutions relate, respectively, to the energy and the entropy of the Ising model at the zero temperature.

Our results indicate that even for a simple process such as binary symmetric HMM, MAP yields a very rich and non-trivial solution structure. The main findings can be summarized as follows: For a small, but finite range of noise values the MAP solution is uniquely related to the observed sequence, and the accuracy of the solution degrades linearly with increasing the noise strength. For intermediate values of noise the accuracy is nearly noise-independent, but now there are exponentially many solutions to the estimation problem, which is reflected in non–zero ground–state entropy for the Ising model. Finally, for larger noise intensities the number of solutions is reduced again, but the accuracy is poor. Furthermore, those regimes are the manifestation of different thermodynamic phases of the Ising model, which are related to each other via first-order phase transitions.

The rest of the paper is organized as follows: After some general discussion of MAP scheme in Section 2, we define the model studied here in Section 3. Its solution is given in Sections 4 and 5. The latter also discusses our concrete findings on the structure of MAP for the binary symmetric HMM. We conclude the paper by discussion of our results and future work.

## 2 Maximum a posteriori (MAP) estimation: general description

Let and be realizations of discrete-time random processes and , respectively. We write their probabilities as and . We assume that is the noisy observation of ; the influence of noise is described by the conditional probability . Let us further assume that we are given an observed sequence , and we know the probabilities and . We do not know which specific sequence generated the observation . MAP offers a method for estimating the generating sequence on the ground of : is found by maximizing over the posterior probability . Since does not depend on , we can equally well minimize

(1) |

The advantage of using is that if is ergodic (in the sense of weak law of large numbers) [2], which we assume from now on, then for , will be independent from , if belongs to the typical set [2]. The typical set has the overall probability converging to one: . Since all elements of have (nearly) equal probability, we can employ with probability one the averaged quantity instead of .

If the noise is very weak, (with being the Kronecker delta), we recover the generating sequence almost exactly. For a strong noise the estimation is dominated by the prior distribution , so that the estimation is not informative. When no priors are put, , the MAP estimation reduces to the Maximum Likelihood (ML) estimation scheme. The latter also reproduces the source sequence almost exactly if the noise is weak.

According to the Viterbi algortithm, for a given the mimimization of in (1) produces one single estimate . However, it is possible that there are other sequences for which , though greater than , is almost equal to the latter in the sense of . All such sequences are equivalent for and we list them as possible solutions:

(2) |

If , we repeat the above ergodicity argument and get for the logarithm of the number of solutions corresponding to a typical observed sequence

(3) |

A finite means that there are exponentially many outcomes of minimizing over . We call entropy, since it relates to the entropy of the Ising model; see below.

We can calculate various moments of , which are random variables due to the dependence on , and employ them for characterizing the accuracy of the estimation; see below for examples. For small noise values these moments will be close to those of the original process . Another useful quantity is the average overlap between the estimated sequences , and the observed sequence (definition of overlap is clarified below). A small overlap means that the estimation is not dominated by observations.

## 3 Binary symmetric hidden markov model (HMM).

### 3.1 Definition.

We consider the MAP estimation of a binary, discrete-time Markov stochastic process . Each random variable has only two realizations . The Markov feature implies

(4) |

where is a time-independent transition probability of the Markov process. For the considered binary symmetric situation it is parameterized by a single number , , , and the stationary distribution is . Furthermore, the noise process is assumed to be memory-less, time-independent and unbiased:

(5) |

where , , and is the probability of error. Here memory-less refers to the factorization in (5), time-independence refers to the fact that in (5) does not depend on , while unbiased means that the noise acts symmetrically on both realizations of the Markov process: .

Note that the composite process with realizations is Markov with transition probabilities

(6) |

However, is in general not a Markov process.

### 3.2 Mapping to the Ising model.

Let us represent the transition probabilities as

(7) |

Likewise, we represent the noise model as

(8) |

We combine (1, 4–5) to represent the log–likelihood as

(9) |

where we have omitted an irrelevant additive factor. is the Hamiltonian of a one–dimensional (1d) Ising spin model with external random fields governed by the probability [11]. The factor in (9) is the spin–spin interaction constant, uniquely determined from the transition probability : If , the constant is positive, which refers to the ferromagnetic situation: the spin–spin interaction tends to align the spins. From now on we assume . We note that the main difference between (9) and other random–field Ising models considered in literature [8, 1], is that in our situation the random fields are not uncorrelated random variables, but display non-Markovian correlations.

### 3.3 Implementation of MAP

To minimize over , we introduce a non-zero temperature , and define the following conditional probability

(10) |

where is the partition function. In the terminology of statistical physics, gives the probability distribution of states for a system with Hamiltonian interacting with a thermal bath at temperature , and with frozen (i.e., fixed for each site) random fields [7]. For , and a given , the function is strongly picked at those [ground states], which minimize . If, however, the limit is taken after the limit , we get

(11) |

where and were defined in (2). From now on we understand the limit in this sense.

The average of [average energy] in the limit will be equal to the minimized over :

(12) |

where we have used the fact that all ground state configurations have the same energy, , for any .

The average logarithm of the number of MAP solutions is equal to the zero-temperature entropy

Let us introduce the the free energy:

(13) |

defined with the Ising Hamiltonian (9). The entropy is expressed via the free energy as [see (3, 12)]:

(14) |

Furthermore, we define the following relevant characteristics of MAP:

(15) |

Here accounts for the correlations between neighbouring spins in the estimated sequence, while measures the overlap between the estimated and the observed sequences (the average Hamming distance between the two is simply ). In the limiting case of very weak noise, when the magnitude of the random fields is large [see (8)], we have (observation-dominance), while is equal to the corresponding value of the Markov process :

(16) |

And for very strong noise (the probability of error is close to , which means ), nullifies, while goes to the corresponding values calculated over the prior distribution : .

## 4 Recursion relation

Let us return to the partition function (10)

We apply to to the following transformations [1]:

where , , and where

(17) | |||

(18) |

Thus, once the first spin is excluded, the field acting on the second spin changes from to . Note the zero-temperature () limits ()

(19) | |||

(20) |

where for and for .

Repeating the above steps we express the partition function as follows:

(21) |

where is obtained from the recursion relation

(22) |

This is a random recursion relation, since are random quantities governed by the probability . Depending on the value of , can take values or . Even when assumes a finite number of values, from (22) can in principle assume an infinite number of values. Fortunately, for , due to the special form (19, 20) of and , the number of values assumed by is finite (though it can be large). It is checked by inspection that the values taken by are parametrized as , where is a positive or negative integer, while can assume only three values . It can also be seen that the states are not recurrent: once takes a value with (note that there is a finite probability for that), it shall never return to the states . In the limit , we can completely disregard the states .

Now recall that the process with probabilities is not Markov. To make it Markov we should enlarge it by adding the random variable ; see (6). Here we write the realizations of this auxiliary Markov process as , so as not to mix them with those of original process . ( and have identical statistical characteristics, but these are different processes: is employed merely for making the composite process Markov.) Likewise, we make the process with realizations Markov by enlarging it to . Let us denote this composite Markov process by . Its conditional probabilities read

(23) |

where and refer to the Markov process and the noise, while takes two values and , depending on whether the corresponding transition is allowed or not by recursion (22). Now the task is to find all possible values of , and then to determine . Before turning to this task, we relate the characteristics of the studied MAP estimation to the stationary probabilities of the composite Markov process . First we get from the stationary probabilities . Next we return to (21) and to the definition of free energy (13). Since the composite Markov process will be seen to be ergodic, the free energy can be written as [1]

(24) |

where the summation is taken over all possible [for a given range of ] values of . Once is found, we can apply (3.3, 15).

As for entropy (14) we get from (18, 21)

(25) |

In this expression we should now select the terms which survive and :

where is the Kronecker symbol. In the limit , should – with probability one, i.e., for the elements of the typical set – converge to , provided that the composite Markov process is ergodic. We thus get [1]

(26) |

The physical meaning of this formula is that the zero-temperature entropy can be extensive only when the external field acting on the spin has the same energy as the spin–spin coupling constant ; see (7). If this is the case, then a macroscopic amount of spins is frustrated, i.e., the factors influencing those spins compensate each other, so that their sign is not predetermined even at the zero temperature.

### 4.1 Stationary states of the recursion

For given and define an integer as

(27) |

Note that the case (and there is no upper limit on ) corresponds to . One can check that for each integer the recurrent states assumed by the recursion (22) can be parametrized as

(28) | |||

(29) |

Note the symmetry and . The transitions between these states—which via the binary function determine the transition matrix in (23)—are illustrated in Figs. 1 and 2 for and , respectively. The reader can easily generalize the latter graph to an arbitrary .

We are now prepared to write down from (23, 28, 29) and Fig. 1 the following transition matrix for the composite Markov process with

(30) |

This is a block matrix composed of matrices (hence the actual size of in (30) is ): means the matrix with all its elements equal to , and

(31) | |||

(32) |

where . Note that is equal to the transition matrix of the Markov process . Once the stationary probability vector of is found from , we get , , , and .

For a general the following matrices serve as building blocks for the matrix

where all not indicated elements are equal to zero. The transition matrix for a general reads [see Figs. 1, 2]

(33) | |||

The left matrices of each tensor product is a block matrix; each block consists of one matrix. The right matrices of each tensor product are also block matrices; now each block consists of one matrix. The zero matrix is written as . The overall size of is , since each state in (28) is augmented by two realizations of the hidden Markov process.

Note that going from one value of to another amounts to changing the dimension of the matrices , , and . Since these matrices are sparse, efficient numerical algorithms of treating them are available, even for larger values of .

## 5 MAP inference

Let us indicate how the quantities of interest are expressed via the stationary probability of the Markov process (obtained from (33)). Recall that since the estimated process is unbiased, we are interested in the second moment , overlap and entropy . The former two quantities have to be obtained via the free energy. To this end, we trace out the redundant variables in the stationary probability of to obtain the following probabilities ():

(34) |

where the equalities in (34) are due to the symmetry of the unbiased situation. We add a lower index to relevant quantities (e.g., to ’s) to indicate the specific value of . Recall that, e.g., and are in general different quantities, since they belong to different Markov processes and , respectively.

Due to (34), we shall need only the probabilities and that normalize to one-half:

(35) |

The free energy then reads (see (20, 24))

(36) |

Now we make use of the fact that free energy is a continuous function of its parameters ^{1}^{1}1Outside phase
transitions free energy is smooth, while at the phase-transition
points it has to be at least continuous, since, besides
being the generating function for calculating various averages, free
energy is also a measure of dynamic stability, and at the phase-transition points both
phases are equally stable by definition (see [7] for more details)., which in our case implies

(37) |

This leads from (36) to

(38) |

One can confirm (38) from (43, 44, 45). Note that (38) will hold for all values of and , since it does not depend on and/or (the formalism holds without requiring any specific relation between and ). Combining (36) with (15) from Section 3.3, we obtain for the second moment of the estimated sequence and the overlap

(39) |

As seen from (28, 29), if the relations (27) hold [recall that they are strict inequalities], there are only two realizations and , which, according to (26), contribute into the entropy. Recalling also (34), we get ()

(40) |

This relation holds for , if we assume .

At the transition points between the various regimes (27), there are more states that contribute into the entropy. The reader can verify that

(41) |

This equation is written down assuming that the value of at does not depend on whether the latter point is reached as or as . This assumptions leads from (41) to:

(42) |

This relation has the same origin as the continuity of the free energy.

### 5.1 The regime or .

We deduce for the stationary probabilities from (33)

(43) | |||

(44) |

This implies from (39) , . This is in fact the Maximum Likelihood (ML) regime: the noise is so small (or is so large) that the estimation is completely governed by the observations: . The second moment of the estimated sequence in this regime is given by the original value (see (16)) times the squared error probability . The entropy in this regime is zero (see (40)): . In this sense the estimation is uniquely determined by observations. We stress that the ML and MAP schemes agree with each other not only for very small, but also for finite noises.

### 5.2 The regimes and .

For more compact presentation of the probabilities, let us introduce separate notations for the noise strength and the Markov correlator , , where and ; see (16). The probabilities obtained from (33) are written as

(45) |

We skip a tedious analytical expressions for .

The values of and deduced from (39, 45) are shown in Figs. 3 and 3. We compare those values with the results obtained by actually finding the MAP estimate via the Viterbi algorithm, and calculating those quantities directly. It is seen that at the regime change points and , and experience sudden jumps, or first-order phase transitions. Remarkably, those features are perfectly reproduced in the simulations, as shown in Figure 3 and 3. For instance, in the ML regime (), the overlap indicating that the estimation is governed solely by observations. At it jumps sharply, and then monotonically decreases in the regime . More generally, decays, both monotonously and via jumps, towards the prior-dominated value .

Consider the second moment of the estimated sequence shown in Figure 3. We see that is nearly a constant for each given . This is the main virtue of MAP scheme as compared to the ML scheme: While the latter predicts a that quickly decays with the noise as (the dotted line in the plot), the proper MAP value of is not far from its noise-free value , and is nearly a constant for a finite range of noise strength . This advantage of MAP over ML is due to supporting the estimation process by the priors. Indeed, the values of the overlap indicate that the estimated sequence is not completely driven by the observations, though it is still not very far from them. Upon increasing towards its maximal value , experiences jumps during each regime change. For larger these jumps are smaller and more frequent, leading to its prior-dominated value .

0.07308 | 0.06587 | 0.05925 | 0.05349 |

Now let us focus on the entropy: It naturally nullifies in the ML regime (), while in the regime the entropy is monotonously increasing with , as shown in Figure 3. At (the phase transition point, where changes from 1 to 2) experiences a jump, which is again usual for first-order phase-transitions. maximizes at an intermediate value of , and then decays to zero for , see Table 1; at this point the present approach reduces to a ferromagnetic 1d Ising model without magnetic fields. This model has a trivial ground-state structure and hence zero entropy. We also note that right at the transition points the values of is different; see Table 2. The largest value is attained for .

0.4051 | |||||

0.1629 | 0.1462 | 0.1220 | 0.0992 | 0.0831 |

Finally, we would like to note that the second moment of the estimated sequence, , is an indirect measure of accuracy. In practice, one is restricted to use such indirect measures as information about the true sequence might not available. In Figure 4 we present the average error rate for the MAP estimation, which is given by the normalized Hamming distance between the true and Viterbi–decoded sequences, plotted against the noise intensity. Also shown is the average error rate of ML estimation, which is simply . For vanishing noise, both MAP and ML yield the same average error. Upon increasing the noise intensity, the MAP estimation error behaves differently depending on the parameter : For small values of , MAP is always superior to ML for a wide range of noise intensities. For larger values of , however, the situation is more complicated: Although both methods perform similarly, there are some differences and crossovers between the two at intermediate noise intensities, as shown in Figure 4.

## 6 Discussion

We theoretically examined Maximum a Posteriori (MAP) estimation for hidden Markov sequences, and found that MAP yields a non–trivial solution structure even for the simple binary and unbiased hidden Markov process considered here. We demonstrated that for a finite range of noise intensities, there is no difference between MAP and Maximum Likelihood (ML) estimations, as the solution is observation–dominated. While it was expected that the two methods agree for a vanishing noise, the fact of their exact agreement for a finite range of the noise is non-trivial. Furthermore, upon increasing the noise intensity the MAP solution switches between different operational regimes that are separated by first–order phase transitions. In particular, a first-order phase-transition separates the regime where MAP and ML agree exactly. At this transition point the influence of the prior information becomes comparable to the influence of observations. In the vicinity of the first-order phase-transitions the performance of MAP (e.g., characteristics of the estimated sequence) changes abruptly. This means that a small change in the noise intensity may lead to a large change in the performance. In other words, the phase-transition points should be avoided in applications.

For practical applications of HMM (e.g., in speech recognition, or machine translation) it is not enough to know the single solution that provides the largest posteriori probability [5]. At the very least, one should also know how many sequences have a posteriori probabilities sufficiently close to the optimal one. Motivated by this fact, we studied the number of MAP-solutions that have (for long sequences ) almost equal logarithms of the posterior probability. A finite means that there is an exponential number of solutions with posterior probabilities slightly less than the optimal. We found that is finite whenever MAP differs from ML. We believe that this theoretical result might have practical implications as well. For instance, in applications such as statistical machine translation, one usually considers top solutions to the inference problem, and then chooses one according to some heuristics. Our result suggests that one needs to be careful with this practice whenever is non–zero, as one might discard a large number of nearly optimal solutions if is not chosen sufficiently large.

We also note that our work is directly related to the notion of trackability, which can be intuitively defined as one’s ability to (accurately) track certain stochastic processes [3, 10]. In fact, a similar phase–transition in the number of solutions was reported by Crespi et. al. [3] for so called weak models, where the entries in the HMP transition and emission matrices are either or . For more general stochastic processes, an information-theoretical characterization of trackability was suggested in [10]. Within this approach, the accuracy is characterized by the probability of the estimated sequence not being equal to the actual one, while the structure of the solution space is described via the number of elements in the (conditional) typical set of sequences given an observed sequence (complexity). Both these quantities relate to the conditional entropy . We note that whereas the accuracy and the complexity measures of [10] deteriorate even for a small (but generic) noise intensity, our approach of defining trackability in terms of the zero–temperature entropy of the Ising Hamiltonian (Equation 9) suggests that a process can be trackable in the MAP sense even in the presence of moderate noise.

Finally, we would like to note that another interesting feature of the MAP estimation is that its characteristics ( and ) change only slightly in between the phase-transition points. In contrast to ML estimation, which deteriorates (at least linearly) when increasing the noise, MAP estimation is stable for a finite range of noise intensities. Thus, although MAP estimation may be less accurate compared to ML, it can be still useful as far as its stability is concerned, provided that its range of application is selected carefully.

There are several directions for further developments. First, we intend to obtain analytical results for the average error rate to complement the empirical analysis presented in Figure 4. Furthermore, one can think of a semi–supervised MAP estimation, where one has (possibly noisy) knowledge about the states of the hidden process at particular times. Remarkably, the framework presented here allows a natural generalization to this case. Indeed, one simply needs to modify the Ising energy function by adding quenched fields at the corresponding locations in the chain. Finally, it will be interesting to generalize the analysis presented here beyond the binary hidden Markov processes considered here. In this case, the MAP optimization problem can be mapped to a Potts model. We would like to note that the behavior observed in the simple binary model can be explained by the emergence of a finite fraction of “frustrated” spins, where the frustration can be attributed to two competing tendencies – accommodating observations on one hand, and the hidden (Markovian) dynamical model on the other. Since this mechanism is rather general, we believe that most features of the MAP scheme uncovered here via an exact analysis of the simplest binary model will survive in more general situations.

#### Acknowledgements

Armen Allahverdyan thanks USC Information Sciences Institute for hospitality and support. This research was supported by the U.S. ARO MURI grant W911NF–06–1–0094.

## References

- [1] U. Behn and V.A. Zagrebnov, One-dimensional Markovian-field Ising model: physical properties and characteristics of discrete stochastic mapping, J. Phys. A 21, 2151 (1988).
- [2] T. M. Cover and J. A. Thomas, Elements of Information Theory, (Wiley, New York, 1991).
- [3] V. Crespi, G. Cybenko, G., and G. Jiang, The theory of trackability with applications to sensor networks, ACM Trans. Sen. Netw. 4, 3 (2008).
- [4] Y. Ephraim and N. Merhav, Hidden Markov processes, IEEE Trans. Inf. Th., 48, 1518-1569, (2002).
- [5] L. A. Foreman, Generalization of the Viterbi algorithm, IMA J. Math. Appl. Bus. Ind., 4, 351 (1993).
- [6] S. Geman and D. Geman, Stochastic relaxation, Gibbs distribution and the Bayesian restoration of images, IEEE Trans. Pattern Analysis Machine Intell. 6, 721 (1984).
- [7] L.D. Landau and E.M. Lifshitz, Statistical Physics, I, (Pergamon Press, Oxford, 1978).
- [8] J.M. Pryce and A.D. Bryce, Statistical mechanics of image restoration, J. Phys. A 28, 511 (1995).
- [9] L. R. Rabiner, A tutorial on hidden Markov models and selected applications in speech recognition, Proc. IEEE, 77, 257-286, (1989).
- [10] Y. Sheng, G. Cybenko, V. Crespi, and G. Jiang, Trackability analysis of multiple processes using multi-distributed agents, International Conference on Integration of Knowledge Intensive Multi-Agent Systems, (2005).
- [11] O. Zuk, I. Kanter and E. Domany, The Entropy of a Binary Hidden Markov Process, J. Stat. Phys. 121, 343 (2005).