Comparative Analysis of Viterbi Training and Maximum Likelihood Estimation for HMMs

We present an asymptotic analysis of Viterbi Training (VT) and contrast it with a more conventional Maximum Likelihood (ML) approach to parameter estimation in Hidden Markov Models. While ML estimator works by (locally) maximizing the likelihood of the observed data, VT seeks to maximize the probability of the most likely hidden state sequence. We develop an analytical framework based on a generating function formalism and illustrate it on an exactly solvable model of HMM with one unambiguous symbol. For this particular model the ML objective function is continuously degenerate. VT objective, in contrast, is shown to have only finite degeneracy. Furthermore, VT converges faster and results in sparser (simpler) models, thus realizing an automatic Occam's razor for HMM learning. For more general scenario VT can be worse compared to ML but still capable of correctly recovering most of the parameters.

READ FULL TEXT VIEW PDF

page 1

page 2

page 3

page 4

04/11/2018

Maximum likelihood estimation in hidden Markov models with inhomogeneous noise

We consider parameter estimation in hidden finite state space Markov mod...
10/27/2020

Fast Stochastic Quadrature for Approximate Maximum-Likelihood Estimation

Recent stochastic quadrature techniques for undirected graphical models...
06/08/2021

General-order observation-driven models: ergodicity and consistency of the maximum likelihood estimator

The class of observation-driven models (ODMs) includes many models of no...
06/30/2017

Neural Sequence Model Training via α-divergence Minimization

We propose a new neural sequence model training method in which the obje...
12/19/2017

Approximate Profile Maximum Likelihood

We propose an efficient algorithm for approximate computation of the pro...
11/07/2018

Iterative Marginal Maximum Likelihood DOD and DOA Estimation for MIMO Radar in the Presence of SIRP Clutter

The spherically invariant random process (SIRP) clutter model is commonl...
12/26/2018

BlinkML: Efficient Maximum Likelihood Estimation with Probabilistic Guarantees

The rising volume of datasets has made training machine learning (ML) mo...

1 Introduction

Hidden Markov Models (HMM) provide one of the simplest examples of structured data observed through a noisy channel. The inference problems of HMM naturally divide into two classes [20, 9]: i) recovering the hidden sequence of states given the observed sequence, and ii)

estimating the model parameters (transition probabilities of the hidden Markov chain and/or conditional probabilities of observations) from the observed sequence. The first class of problems is usually solved via the maximum a posteriori (MAP) method and its computational implementation known as Viterbi algorithm

[20, 9]

. For the parameter estimation problem, the prevailing method is maximum likelihood (ML) estimation, which finds the parameters by maximizing the likelihood of the observed data. Since global optimization is generally intractable, in practice it is implemented through an expectation–maximization (EM) procedure known as Baum–Welch algorithm 

[20, 9].

An alternative approach to parameter learning is Viterbi Training (VT), also known in the literature as segmental K-means, Baum–Viterbi algorithm, classification EM, hard EM, etc. Instead of maximizing the likelihood of the observed data, VT seeks to maximize the probability of the most likely hidden state sequence. Maximizing VT objective function is hard 

[8], so in practice it is implemented via an EM-style iterations between calculating the MAP sequence and adjusting the model parameters based on the sequence statistics. It is known that VT lacks some of the desired features of ML estimation such as consistency [17], and in fact, can produce biased estimates [9]. However, it has been shown to perform well in practice, which explains its widespread use in applications such as speech recognition [16], unsupervised dependency parsing [24], grammar induction [6], ion channel modeling [19]. It is generally assumed that VT is more robust and faster but usually less accurate, although for certain tasks it outperforms conventional EM [24].

The current understanding of when and under what circumstances one method should be preferred over the other is not well–established. For HMMs with continuos observations, Ref. [18] established an upper bound on the difference between the ML and VT objective functions, and showed that both approaches produce asymptotically similar estimates when the dimensionality of the observation space is very large. Note, however, that this asymptotic limit is not very interesting as it makes the structure imposed by the Markovian process irrelevant. A similar attempt to compare both approaches on discrete models (for stochastic context free grammars) was presented in [23]. However, the established bound was very loose.

Our goal here is to understand, both qualitatively and quantitatively, the difference between the two estimation methods. We develop an analytical approach based on generating functions for examining the asymptotic properties of both approaches. Previously, a similar approach was used for calculating entropy rate of a hidden Markov process [1]. Here we provide a non-trivial extension of the methods that allows to perform comparative asymptotic analysis of ML and VT estimation. It is shown that both estimation methods correspond to certain free-energy minimization problem at different temperatures. Furthermore, we demonstrate the approach on a particular class of HMM with one unambiguous symbol and obtain a closed–form solution to the estimation problem. This class of HMMs is sufficiently rich so as to include models where not all parameters can be determined from the observations, i.e., the model is not identifiable [7, 14, 9].

We find that for the considered model VT is a better option if the ML objective is degenerate (i.e., not all parameters can be obtained from observations). Namely, not only VT recovers the identifiable parameters but it also provides a simple (in the sense that non-identifiable parameters are set to zero) and optimal (in the sense of the MAP performance) solution. Hence, VT realizes an automatic Occam’s razor for the HMM learning. In addition, we show that the VT algorithm for this model converges faster than the conventional EM approach. Whenever the ML objective is not degenerate, VT leads generally to inferior results that, nevertheless, may be partially correct in the sense of recovering certain (not all) parameters.

2 Hidden Markov Process

Let be a discrete-time, stationary, Markov process with conditional probability

(1)

where is an integer. Each realization

of the random variable

takes values . We assume that is mixing: it has a unique stationary distribution , , that is established from any initial probability in the long time limit.

Let random variables , with realizations , be noisy observations of : the (time-invariant) conditional probability of observing given the realization of the Markov process is . Defining , , the joint probability of reads

(2)

where the transfer-matrix with matrix elements is defined as

(3)

is called a hidden Markov process. Generally, it is not Markov, but it inherits stationarity and mixing from [9]. The probabilities for can be represented as follows:

(4)

where is a product of transfer matrices.

3 Parameter Estimation

3.1 Maximum Likelihood Estimation

The unknown parameters of an HMM are the transition probabilities of the Markov process and the observation probabilities ; see (2). They have to be estimated from the observed sequence . This is standardly done via the maximum-likelihood approach: one starts with some trial values , of the parameters and calculates the (log)-likelihood , where means the probability (4) calculated at the trial values of the parameters. Next, one maximizes over and for the given observed sequence (in practice this is done via the Baum-Welch algorithm [20, 9]). The rationale of this approach is as follows. Provided that the length of the observed sequence is long, and recaling that is mixing (due to the analogous feature of

) we get probability-one convergence (law of large numbers)

[9]:

(5)

where the average is taken over the true probability that generated . Since the relative entropy is non-negative, , the global maximum of as a function of and is reached for and . This argument is silent on how unique this global maximum is and how difficult to reach it.

3.2 Viterbi Training

An alternative approach to the parameter learning employs the maximal a posteriori (MAP) estimation and proceeds as follows: Instead of maximizing the likelihood of observed data (5) one tries to maximize the probability of the most likely sequence [20, 9]. Given the joint probability at trial values of parameters, and given the observed sequence , one estimates the generating state-sequence via maximizing the a posteriori probability

(6)

over . Since does not depend on , one can maximize . If the number of observations is sufficiently large , one can substitute by its average over [see (5)] and instead maximize (over model parameters)

(7)

To relate (7) to the free energy concept (see e.g. [2, 4]), we define an auxiliary (Gibbsian) probability

(8)

where is a parameter. As a function of (and for a fixed ), concentrates on those that maximize :

(9)

where is the Kronecker delta, are equivalent outcomes of the maximization, and is the number of such outcomes. Further, define

(10)

Within statistical mechanics Eqs. 8 and 10 refer to, respectively, the Gibbs distribution and free energy of a physical system with Hamiltonian coupled to a thermal bath at inverse temperature [2, 4]. It is then clear that ML and Viterbi Training correspond to minimizing the free energy Eq. 10 at and , respectively. Note that , which yields .

3.3 Local Optimization

As we mentioned, global maximization of neither objective is feasible in the general case. Instead, in practice this maximization is locally implemented via an EM-type algorithm [20, 9]: for a given observed sequence , and for some initial values of the parameters, one calculates the expected value of the objective function under the trial parameters (E-step), and adjusts the parameters to maximize this expectation (M-step). The resulting estimates of the parameters are now employed as new trial parameters and the previous step is repeated. This recursion continues till convergence.

For our purposes, this procedure can be understood as calculating certain statistics of the hidden sequence averaged over the Gibbs distribution Eqs. 8. Indeed, let us introduce and define

(11)

Then, for instance, the (iterative) Viterbi estimate of the transition probabilities are given as follows:

(12)

Conditional probabilities for observations are calculated similarly via a different indicator function.

4 Generating Function

Note from (4) that both and are obtained as matrix-products. For a large number of multipliers the behavior of such products is governed by the multiplicative law of large numbers. We now recall its formulation from [10]: for and generated by the mixing process there is a probability-one convergence:

(13)

where is a matrix norm in the linear space of matrices, and

is the maximal eigenvalue of

. Note that (13) does not depend on the specific norm chosen, because all norms in the finite-dimensional linear space are equivalent; they differ by a multiplicative factor that disappears for [10]. Eqs. (4, 13) also imply . Altogether, we calculate (5) via its probability-one limit

(14)

Note that the multiplicative law of large numbers is normally formulated for the maximal singular value. Its reformulation in terms of the maximal eigenvalue needs additional arguments

[1].

Introducing the generating function

(15)

where is a non-negative number, and where means in degree of , one represents (14) as

(16)

where we took into account , as follows from (15).

The behavior of is better understood after expressing it via the zeta-function [1]

(17)

where is given by (15). Since for a large , [this is the content of (13)], the zeta-function has a zero at :

(18)

Indeed for close (but smaller than) , the series almost diverges and one has . Recalling that and taking in , we get from (16)

(19)

For calculating in (10) we have instead of (19)

(20)

where employs instead of in (19).

Though in this paper we restricted ourselves to the limit , we stress that the knowledge of the generating function allows to analyze the learning algorithms for any finite .

5 Hidden Markov Model with One Unambiguous Symbol

5.1 Definition

Figure 1: The hidden Markov process (2122) for . Gray circles and arrows indicate on the realization and transitions of the internal Markov process; see (21). The light circles and black arrows indicate on the realizations of the observed process.

Given a -state Markov process , the observed process has two states and ; see Fig. 1. All internal states besides one are observed as , while the internal state produces, respectively, and with probabilities and . For we obtain from (1) , , . Hence is unambiguous: if it is observed, the unobserved process was certainly in ; see Fig. 1. The simplest example of such HMM exists already for ; see [12] for analytical features of entropy for this case. We, however, describe in detail the situation, since this case will be seen to be generic (in contrast to ) and it allows straightforward generalizations to . The transition matrix (1) of a general Markov process reads

(21)

where, e.g., is the transition probability ; see Fig. 1. The corresponding transfer matrices read from (3)

(22)

Eq. (22) makes straightforward the reconstruction of the transfer-matrices for . It should also be obvious that for all only the first row of consists of non-zero elements.

For we get from (22) the simplest example of an aggregated HMM, where several Markov states are mapped into one observed state. This model plays a special role for the HMM theory, since it was employed in the pioneering study of the non-identifiability problem [7].

5.2 Solution of the Model

For this model can be calculated exactly, because has only one non-zero row. Using the method outlined in the supplementary material (see also [1, 3]) we get

(23)

where and are the largest eigenvalues of and , respectively

(24)
(25)

Here and are, respectively right and left eigenvalues of , while () are the eigenvalues of . Eqs. (24, 25) obviously extend to hatted quantities.

We get from (23, 19):

(26)
(27)

Note that for , are return probabilities to the state of the -state Markov process. For this interpretation does not hold, but still has a meaning of probability as .

Turning to equations (19, 27) for the free energy, we note that as a function of trial values it depends on the following parameters:

(28)

As a function of the true values, the free energy depends on the same parameters (28) [without hats], though concrete dependencies are different. For the studied class of HMM there are at most unknown parameters: transition probabilities of the unobserved Markov chain, and one parameter coming from observations. We checked numerically that the Jacobian of the transformation from the unknown parameters to the parameters (28) has rank . Any parameters among (28) can be taken as independent ones.

For the number of effective independent parameters that affect the free energy is smaller than the number of parameters. So if the number of unknown parameters is larger than , neither of them can be found explicitly. One can only determine the values of the effective parameters.

6 The Simplest Non-Trivial Scenario

The following example allows the full analytical treatment, but is generic in the sense that it contains all the key features of the more general situation given above (21). Assume that and

(29)

Note the following explicit expressions

(30)
(31)
(32)

Eqs. (3032) with obvious changes for every symbol hold for , and . Note a consequence of :

(33)

6.1 Optimization of

Eqs. (27) and (3032) imply ,

(34)
(35)

The free energy depends on three independent parameters [recall (33)]. Hence, minimizing we get (), but we do not obtain a definite solution for the unknown parameters: any four numbers , , , satisfying three equations , , , minimize .

6.2 Optimization of

In deriving (35) we used no particular feature of , , . Hence, as seen from (20), the free energy at is recovered from (35) by equating its LHS to and by taking in its RHS: , , , . The zero-temperature free energy reads from (35)

(36)

We now minimize over the trial parameters . This is not what is done by the VT algorithm; see the discussion after (12). But at any rate both procedures aim to minimize the same target. VT recursion for this models will be studied in section 6.3 — it leads to the same result. Minimizing over the trial parameters produces four distinct solutions:

(37)

For each of these four solutions: () and . The easiest way to get these results is to minimize under conditions (for ), obtain and then to conclude that due to the inequality the conditional minimization led to the global minimization. The logics of (37) is that the unambiguous state tends to get detached from the ambiguous ones, since the probabilities nullifying in (37) refer to transitions from or to the unambiguous state.

Note that although minimizing either and produces correct values of the independent variables , , , in the present situation minimizing is preferable, because it leads to the four-fold degenerate set of solutions (37) instead of the continuously degenerate set. For instance, if the solution with is chosen we get for other parameters

(38)

Furthermore, a more elaborate analysis reveals that for each fixed set of correct parameters only one among the four solutions Eq. 37 provides the best value for the quality of the MAP reconstruction, i.e. for the overlap between the original and MAP-decoded sequences.

Finally, we note that minimizing allows one to get the correct values of the independent variables , and only if their number is less than the number of unknown parameters. This is not a drawback, since once the number of unknown parameters is sufficiently small [less than four for the present case (29)] their exact values are obtained by minimizing . Even then, the minimization of can provide partially correct answers. Assume in (36) that the parameter is known, . Now has three local minima given by , and ; cf. with (37). The minimum with is the global one and it allows to obtain the exact values of the two effective parameters: and . These effective parameters are recovered, because they do not depend on the known parameter . Two other minima have greater values of , and they allow to recover only one effective parameter: . If in addition to also is known, the two local minimia of ( and ) allow to recover only. In contrast, if (or ) is known exactly, there are three local minima again—, , —but now none of effective parameters is equal to its true value: ().

6.3 Viterbi EM

Recall the description of the VT algorithm given after (12). For calculating via (11, 12) we modify the transfer matrix element in (15, 17) as , which produces from (11, 12) for the MAP-estimates of the transition probabilities

(39)
(40)
(41)

where , . The limit of and is obvious: each of them is equal to or depending on the ratios and . The EM approach amounts to starting with some trial values , , , and using , , , as new trial parameters (and so on). We see from (3941) that the algorithm converges just in one step: (3941) are equal to the parameters given by one of four solutions (37)—which one among the solutions (37) is selected depends on the on initial trial parameters in (3941)—recovering the correct effective parameters (3032); e.g. cf. (38) with (39, 41) under . Hence, VT converges in one step in contrast to the Baum-Welch algorithm (that uses EM to locally minimize ) which, for the present model, obviously does not converge in one step. There is possibly a deeper point in the one-step convergence that can explain why in practice VT converges faster than the Baum-Welch algorithm [9, 21]: recall that, e.g. the Newton method for local optimization works precisely in one step for quadratic functions, but generally there is a class of functions, where it performs faster than (say) the steepest descent method. Further research should show whether our situation is similar: the VT works just in one step for this exactly solvable HMM model that belongs to a class of models, where VT generally performs faster than ML.

We conclude this section by noting that the solvable case (29) is generic: its key results extend to the general situation defined above (21). We checked this fact numerically for several values of . In particular, the minimization of nullifies as many trial parameters as necessary to express the remaining parameters via independent effective parameters . Hence for and two such trial parameters are nullified; cf. with discussion around (28). If the true error probability , the trial value is among the nullified parameters. Again, there is a discrete degeneracy in solutions provided by minimizing .

7 Summary

We presented a method for analyzing two basic techniques for parameter estimation in HMMs, and illustrated it on a specific class of HMMs with one unambiguous symbol. The virtue of this class of models is that it is exactly solvable, hence the sought quantities can be obtained in a closed form via generating functions. This is a rare occasion, because characteristics of HMM such as likelihood or entropy are notoriously difficult to calculate explicitly [1]. An important feature of the example considered here is that the set of unknown parameters is not completely identifiable in the maximum likelihood sense [7, 14]. This corresponds to the zero eigenvalue of the Hessian for the ML (maximum-likelihood) objective function. In practice, one can have weaker degeneracy of the objective function resulting in very small values for the Hessian eigenvalues. This scenario occurs often in various models of physics and computational biology [11]. Hence, it is a drawback that the theory of HMM learning was developed assuming complete identifiably [5].

One of our main result is that in contrast to the ML approach that produces continuously degenerate solutions, VT results in finitely degenerate solution that is sparse, i.e., some [non-identifiable] parameters are set to zero, and, furthermore, converges faster. Note that sparsity might be a desired feature in many practical applications. For instance, imposing sparsity on conventional EM-type learning has been shown to produce better results part of speech tagging applications [25]. Whereas  [25] had to impose sparsity via an additional penalty term in the objective function, in our case sparsity is a natural outcome of maximizing the likelihood of the best sequence. While our results were obtained on a class of exactly-solvable model, it is plausible that they hold more generally.

The fact that VT provides simpler and more definite solutions—among all choices of the parameters compatible with the observed data—can be viewed as a type of the Occam’s razor for the parameter learning. Note finally that statistical mechanics intuition behind these results is that the aposteriori likelihood is (negative) zero-temperature free energy of a certain physical system. Minimizing this free energy makes physical sense: this is the premise of the second law of thermodynamics that ensures relaxation towards a more equilibrium state. In that zero-temperature equilibrium state certain types of motion are frozen, which means nullifying the corresponding transition probabilities. In that way the second law relates to the Occam’s razor. Other connections of this type are discussed in [15].

Acknowledgments

This research was supported in part by the US ARO MURI grant No. W911NF0610094 and US DTRA grant HDTRA1-10-1-0086.

References

  • [1] A. E. Allahverdyan, Entropy of Hidden Markov Processes via Cycle Expansion, J. Stat. Phys. 133, 535 (2008).
  • [2] A.E. Allahverdyan and A. Galstyan,

    On Maximum a Posteriori Estimation of Hidden Markov Processes

    , Proc. of UAI, (2009).
  • [3] R. Artuso. E. Aurell and P. Cvitanovic, Recycling of strange sets, Nonlinearity 3, 325 (1990).
  • [4] P. Baldi and S. Brunak, Bioinformatics, MIT Press, Cambridge, USA (2001).
  • [5] L. E. Baum and T. Petrie, Statistical inference for probabilistic functions of finite state Markov chains, Ann. Math. Stat. 37, 1554 (1966).
  • [6] J.M. Benedi, J.A. Sanchez, Estimation of stochastic context-free grammars and their use as language models, Comp. Speech and Lang. 19, pp. 249-274 (2005).
  • [7] D. Blackwell and L. Koopmans, On the identifiability problem for functions of finite Markov chains, Ann. Math. Statist. 28, 1011 (1957).
  • [8] S. B. Cohen and N. A. Smith, Viterbi training for PCFGs: Hardness results and competitiveness of uniform initialization, Procs. of ACL (2010).
  • [9] Y. Ephraim and N. Merhav, Hidden Markov processes, IEEE Trans. Inf. Th., 48, 1518-1569, (2002).
  • [10] L.Y. Goldsheid and G.A. Margulis, Lyapunov indices of a product of random matrices, Russ. Math. Surveys 44, 11 (1989).
  • [11] R. N. Gutenkunst et al., Universally Sloppy Parameter Sensitivities in Systems Biology Models, PLoS Computational Biology, 3, 1871 (2007).
  • [12] G. Han and B. Marcus, Analyticity of entropy rate of hidden Markov chains, IEEE Trans. Inf. Th., 52, 5251 (2006).
  • [13] R. A. Horn and C. R. Johnson, Matrix Analysis (Cambridge University Press, New Jersey, USA, 1985).
  • [14] H. Ito, S. Amari, and K. Kobayashi, Identifiability of Hidden Markov Information Sources, IEEE Trans. Inf. Th., 38, 324 (1992).
  • [15] D. Janzing, On causally asymmetric versions of Occam’s Razor and their relation to thermodynamics, arXiv:0708.3411 (2007).
  • [16] B. H. Juang and L. R. Rabiner, The segmental k-means algorithm for estimating parameters of hidden Markov models, IEEE Trans. Acoustics, Speech, and Signal Processing, ASSP-38, no.9, pp.1639-1641, (1990).
  • [17] B. G. Leroux, Maximum-Likelihood Estimation for Hidden Markov Models, Stochastic Processes and Their Applications, 40, 127 (1992).
  • [18] N. Merhav and Y. Ephraim, Maximum likelihood hidden Markov modeling using a dominant sequence of states, IEEE Transactions on Signal Processing, vol.39, no.9, pp.2111-2115 (1991).
  • [19] F. Qin, Restoration of single-channel currents using the segmental k-means method based on hidden Markov modeling, Biophys J 86, 1488–1501 (2004).
  • [20] L. R. Rabiner, A tutorial on hidden Markov models and selected applications in speech recognition, Proc. IEEE, 77, 257 (1989).
  • [21] L. J. Rodriguez and I. Torres, Comparative Study of the Baum-Welch and Viterbi Training Algorithms

    , Pattern Recognition and Image Analysis, Lecture Notes in Computer Science,

    2652/2003, 847 (2003).
  • [22] D. Ruelle, Statistical Mechanics, Thermodynamic Formalism, (Reading, MA: Addison-Wesley, 1978).
  • [23] J. Sanchez, J. Benedi, F. Casacuberta, Comparison between the inside-outside algorithm and the Viterbi algorithm for stochastic context-free grammars, in Adv. in Struct. and Synt. Pattern Recognition (1996).
  • [24] V. I. Spitkovsky, H. Alshawi, D. Jurafsky, and C. D. Manning, Viterbi Training Improves Unsupervised Dependency Parsing, in Proc. of the 14th Conference on Computational Natural Language Learning (2010).
  • [25] A. Vaswani, A. Pauls, and D. Chiang, Efficient optimization of an MDL-inspired objective function for unsupervised part-of-speech tagging, in Proc. ACL (2010).

Supplementary Material

Here we recall how to calculate the moment-generating function

via zeta-function [22] and periodic orbits [3, 1]. Let be the maximal eigenvalue of matrix with non-negative elements [13]. Since and have identical eigenvalues, we get , ( is an integer).

Recall the content of section 4. Eqs. (15, 3, 4) lead to

(42)
(43)

where we have introduced a notation for better readability. We obtain

(44)

where and are arbitrary sequences of symbols . One can prove for [22]:

where are the indices referring to realizations of the HMM, and where means that the summation goes over all that divide , e.g., for . Here contains sequences selected according to the following rules: i) turns to itself after successive cyclic permutations, but does not turn to itself after any smaller (than ) number of successive cyclic permutations; ii) if is in , then contains none of those sequences obtained from under successive cyclic permutations.

Starting from (Supplementary Material) and introducing notations , , we transform as

The summation over , , yields

(45)

where , (all the notations introduced generalize—via introducing a hat—to functions with trial values of the parameters, e.g., ). are obtained from (45). We write them down assuming that (two realizations of the observed process)