Computing Sequence Likelihoods with Forward Messages

From Hidden Paths to Forward Messages

For a hidden Markov model, the likelihood of an observed sequence is the sum of the probabilities of all hidden-state paths that could have produced it. That sum grows rapidly with the sequence length. Forward messages organize the same calculation by keeping, at each time, one joint probability for each possible current state. A transition update combines those probabilities with the transition matrix; an emission update accounts for the observation at that time. The final likelihood is the sum across states.

Work on a probability space \((\Omega,\mathcal{F},\mathbb{P})\). Let \((X_n)_{n\ge1}\) be a hidden Markov chain with finite state space \(S=\{1,\ldots,K\}\), initial distribution \(\pi\) as a row vector, and transition matrix \(P=(p_{ij})\). Thus \(p_{ij}=\mathbb{P}(X_{n+1}=j\mid X_n=i)\), and each row of \(P\) sums to 1. Let \((Y_n)_{n\ge1}\) take values in a finite observation alphabet \(\mathcal{Y}\). Conditional on the hidden states, observations are independent, and the emission probabilities are \(b_i(y)=\mathbb{P}(Y_n=y\mid X_n=i)\); for each state \(i\), \(\sum_{y\in\mathcal{Y}}b_i(y)=1\). These specifications define the usual HMM joint probability model.

The time-one message uses the initial distribution and the first emission. For each later time, first predict the hidden state using \(P\), then multiply each destination-state probability by the probability of the observation there. Writing \(\odot\) for componentwise multiplication and \(b(y)=(b_1(y),\ldots,b_K(y))\), the updates are

To justify the update, for any state \(j\), sum over the possible state \(i\) at time \(t-1\). The joint probability of the observation prefix through \(t-1\), \(X_{t-1}=i\), and \(X_t=j\) is \(\alpha_{t-1}(i)p_{ij}\), by the Markov transition rule. Given \(X_t=j\), the probability of the new observation is \(b_j(y_t)\), so multiplying by that emission probability gives the new joint probability. This explains both the order of operations and why the messages retain the likelihood information.

Why the Recursion Gives the Likelihood

Proof. At time \(1\), the joint probability of \(Y_1=y_1\) and \(X_1=j\) is \(\pi_jb_j(y_1)\), which is the \(j\)th entry of \(\pi\odot b(y_1)\). Now suppose the identity holds at time \(t-1\). The events \(X_{t-1}=i\) partition the possible preceding states, so for each \(j\), \[ \begin{aligned} \mathbb{P}(Y_1=y_1,\ldots,Y_t=y_t,X_t=j) &=\sum_{i\in S}\mathbb{P}(Y_1=y_1,\ldots,Y_{t-1}=y_{t-1},X_{t-1}=i)\\ &\qquad\cdot p_{ij}b_j(y_t)\\ &=\left(\sum_{i\in S}\alpha_{t-1}(i)p_{ij}\right)b_j(y_t). \end{aligned} \] The last expression is exactly the \(j\)th entry of \((\alpha_{t-1}P)\odot b(y_t)\). This proves the identity at time \(t\), and induction proves it through \(T\). Finally, the events \(X_T=j\), \(j\in S\), are disjoint and cover all possibilities, so summing the joint probabilities over \(j\) gives the sequence likelihood. \(\square\)

The recursion is a method of evaluating the path sum, not an approximation to it. In practice, keep entries as exact fractions when the model probabilities are rational. Check transition rows and emission distributions before calculating; a single invalid row or emission distribution does not define the stated HMM.

Scaling Messages Without Losing the Likelihood

Forward entries may become very small for long sequences, even when calculated exactly. A useful alternative is to normalize the message after each observation and keep track of the removed scale. This changes the size of the intermediate numbers, not the likelihood. Define \(c_1=\sum_j\alpha_1(j)\), \(f_1=\alpha_1/c_1\), and, for \(t\ge2\), define \(u_t=(f_{t-1}P)\odot b(y_t)\), \(c_t=\sum_j u_t(j)\), and \(f_t=u_t/c_t\). These normalizations require each \(c_t>0\).

Proof. At time 1, \(\alpha_1=c_1f_1\). Suppose that at time \(t-1\), \(\alpha_{t-1}=(\prod_{r=1}^{t-1}c_r)f_{t-1}\). The unnormalized forward update gives \[ \alpha_t=(\alpha_{t-1}P)\odot b(y_t) =\left(\prod_{r=1}^{t-1}c_r\right)\left((f_{t-1}P)\odot b(y_t)\right) =\left(\prod_{r=1}^{t-1}c_r\right)u_t =\left(\prod_{r=1}^{t}c_r\right)f_t. \] This proves the relation by induction. Since the entries of \(f_T\) sum to 1, summing the entries of \(\alpha_T\) gives \(\prod_{r=1}^T c_r\). The sum of \(\alpha_T\) is the sequence likelihood by the forward-message identity, proving the proposition. \(\square\)

If some \(c_t=0\), the observation prefix through time \(t\) has probability zero. Then \(f_t=u_t/c_t\) is undefined, and subsequent normalized updates cannot be used; the unnormalized recursion still gives a zero message. This is a model or data compatibility issue, not a reason to divide by zero or replace the zero by an arbitrary positive number.

Reading the Calculation Correctly

A forward message is not the probability of a hidden state after seeing the observations; it is the joint probability of that state and the observation prefix. In particular, its entries need not sum to 1. Dividing by their sum produces a normalized state distribution when that sum is positive, but the removed sum is part of the likelihood and must be retained if the goal is to compute that likelihood.

Common mistake: In the \(B,A\) example, after predicting \((\tfrac7{32},\tfrac9{32})\), one might add these entries and report \(\tfrac12\) as the likelihood. That sum accounts for the first observation \(B\) and the transition, but it omits the final emission \(A\). The correct update multiplies the entries by the state-specific \(A\)-probabilities, giving \((\tfrac{21}{128},\tfrac9{128})\); only then do we sum, obtaining \(\tfrac{15}{64}\). A reliable check at every time is to confirm that the emission for the current observation has been applied exactly once.

For a short sequence, exact fractions make each update and the final sum auditable. For a long sequence, the same recursion applies; use stepwise scaling when its positive-scale condition holds, or an appropriate log-probability implementation. The essential bookkeeping remains the same: transition into a state, apply that state's emission probability, and sum over the possible final states.