Derivation of the Baum-Welch algorithm for HMMs

Consider our formulation of HMMs as before, i.e. $\{v_1, \cdots, v_m\}$ is the alphabet, $O_1, \cdots, O_T$ is the output sequence, s_t is the state at time t, A is the transition probability a_ij = P(s_t=j|s_t-1=i) and B is the state-symbol probability for observing a given output symbol b_jk = P(O_t = v_k|s_t=s).

Let Z = f(A,B) be a random variable such that x_ijk = f(A,B) denotes the joint event consisting of the transition from state i to j and the observation of the v_k from state j. So the observed output is completely specified by a sequence $X = x_{i_1j_1k_1} \cdots x_{i_Tj_Tk_T}$. Let p_ijk = P(x_ijk).

Then the parameter of interest, $\theta$, is exactly Z since it represents the totality of transition and state-symbol probabilities.

We wish to find a $\theta'$ such that

$$\sum_X P(X\vert O) \log P'(X,O) > \sum_X P(X\vert O) \log P(X,O)$$

The expression on the left is called Baum's auxiliary function. We have already seen that if the above inequality holds, then P'(O) > P(O). So maximizing the LHS wrt. $\theta$ is equivalent to maximizing P(O).

To do this, we first find a stationary point of the LHS subject to the constraints $\forall i: \sum_j \sum_k p_{ijk}' - 1 = 0$. Using Lagrange multipliers, the function to maximize with respect to the variables p_ijk' is thus $w = \sum_X P(X\vert O) \log P'(X,O) - \sum_N \lambda_i \left(\sum_j \sum_k p_{ijk}'-1\right)$. Then $\frac{\partial w}{\partial p_{ijk}'} =$

 

$$\frac{\partial }{\partial p_{ijk}'} \left[ \sum_X P(X\vert O) \log P'(X,O) - \sum_N \lambda_i \left(\sum_j \sum_k p_{ijk}'-1\right) \right]$$

$$\sum_X P(X\vert O) \frac{\partial/\partial p_{ijk}' P'(X,O)} {P'(X,O)}-\lambda_i$$ = (16)

Since X completely specifies O (recall that x_ijk also specifies that symbol v_k is observed in state j),

$$\begin{eqnarray*}P'(X,O) = P'(X) &=& \prod_{t=0}^T p_{i_tj_tk_t}'\\ \frac{\part... ...rac{\partial }{\partial p_{ijk}'}\prod_{t=1}^T p_{i_tj_tk_t}'\\ \end{eqnarray*}$$

If c_ijk denotes the number of times x_ijk occurs in X, then

$$\begin{eqnarray*}\frac{\partial }{\partial p_{ijk}'}P'(X,O) &=& \frac{\partial ... ..._tj_tk_t}}}{p_{ijk}'}\\ &=& \frac{c_{ijk} P'(X,O)}{p_{ijk}'}\\ \end{eqnarray*}$$

Hence,

$$\frac{\partial/\partial p_{ijk}' P'(X,O)}{P'(X,O)} = \frac{c_{ijk}}{p_{ijk}'}$$

We can also satisfy ourselves that

$$\partial^2 P'(X,O)/\partial (p_{ijk}')^2 = -c_{ijk}/(p_{ijk}')^2 <0$$

and so the stationary point obtained by setting the first derivative to zero is indeed a maximum.

Substituting in (17) and setting to zero,

$$\begin{eqnarray*}\sum_X P(X\vert O) \frac{c_{ijk}}{p_{ijk}'} -\lambda_i = 0 \qqu... ...}\\ p_{ijk}' = \frac{1}{\lambda_i}\sum_X P(X\vert O) c_{ijk}\\ \end{eqnarray*}$$

If we now write $P(X\vert O) = \frac{P(X,O)}{P(O)}$, then

$$\begin{eqnarray*}p_{ijk}' &=& \frac{1}{\lambda_i P(O)}\sum_X c_{ijk}P(X,O) = \f... ...i} \sum_X c_{ijk} P(X,O)\\ K_ip_{ijk}' &=& \sum_X c_{ijk}P(X,O) \end{eqnarray*}$$

where K_i can be interpreted as a normalizing constant. The RHS can now be summed over each of the times $1, \cdots, T$ by grouping the sequences such that in each group x_itjtkt = x_ijk.

$$\begin{eqnarray*}&=& \sum_{t=1}^T P(O_1,\cdots,O_{t-1}, S_{t-1}=i) p_{ijk} P(O_... ...\frac{\sum_{t=1}^T \alpha_t(i)a_{ij}b_j(v_k)\beta_{t+1}(j)}{K_i} \end{eqnarray*}$$

The RHS above is just the reestimated probability of x_ijk because prior to normalization, the sum yields the expected frequency of the x_ijk computed using the current values of the p_ijk.

We can now proceed to calculate each p_ijk' as above. The set of all such p_ijk's constitutes our $\theta'$ which we can use in place of $\theta$ during the next iteration.

The EM theorem tells us that if we continue in this fashion, then after each iteration

1.

either the reestimated $\theta'$ is more likely than the original $\theta$ in the sense that $P(O\vert\theta') > P(O\vert\theta)$ or

2.

we have reached a stationary point of the likelihood function w at which $\theta' = \theta$

This is an algorithm, because it is iterative. There is an expectation step (aka estimation step) and a maximization step (aka modification step) in each iteration.

Note that the algorithm is not guaranteed to converge at the global maximum. But it has been found to yield good results in practice.