Stable HMM Computation with Scaling and Log Probabilities

Why Likelihood Computation Needs Numerical Stabilization

A hidden Markov model assigns a probability to an observed sequence by summing over all possible hidden-state sequences. Even when each individual probability is represented accurately, the sum can become extremely small as the sequence grows. Directly multiplying and adding such numbers can produce numerical underflow: a computer rounds a positive likelihood to zero, and its logarithm then cannot be computed. Scaling and log probabilities address this problem in two complementary ways.

Let \((\Omega,\mathcal{F},\mathbb{P})\) be a probability space. Consider a discrete-time hidden Markov model with latent process \((Z_t)_{t\ge1}\) on a finite state space \(S=\{1,\ldots,K\}\) and observations \(Y_1,\ldots,Y_T\). The latent chain has initial distribution \(\pi\), written as a row vector, and transition matrix \(A=(a_{ij})\), where \(a_{ij}=\mathbb{P}(Z_{t+1}=j\mid Z_t=i)\). Each row of \(A\) sums to one. Given the latent states, the observations are conditionally independent, and the observation at time \(t\), given \(Z_t=j\), has probability mass function or density \(b_j(y_t)\). For continuous observations, the sequence likelihood below is a density evaluated at the observed values.

Summing the joint probabilities or densities over the hidden paths gives the sequence likelihood:

This formula has \(K^T\) terms. Grouping paths according to their final state gives a practical recursion. Define \(q_t(j)\) as the joint probability or density of \(Y_1=y_1,\ldots,Y_t=y_t\) and \(Z_t=j\). Then \(q_1(j)=\pi_jb_j(y_1)\), and for \(t\ge2\),

The recursion follows by summing the contribution from each possible previous state \(i\). It is the arithmetic being stabilized here; the emphasis is on preserving the likelihood accurately, not on deriving the recursion in full.

Scaling at Every Time Step

Instead of allowing the unnormalized quantities \(q_t(j)\) to shrink, normalize the state weights after each observation. At time \(1\), calculate \(u_1(j)=\pi_jb_j(y_1)\), set \(c_1=\sum_j u_1(j)\), and define \(\alpha_1(j)=u_1(j)/c_1\). At time \(t\ge2\), calculate

Each \(\alpha_t\) is a row vector of nonnegative entries summing to one, provided \(c_t>0\). The factors \(c_t\) retain the probability mass removed by normalization. In particular, they need not be less than one when observations have densities, rather than probabilities.

Proof. Let \(C_0=1\) and \(C_t=\prod_{r=1}^{t}c_r\). We prove by induction that \(q_t(j)=C_t\alpha_t(j)\) for every state \(j\). At \(t=1\), \(q_1(j)=u_1(j)=c_1\alpha_1(j)=C_1\alpha_1(j)\). Suppose the identity holds at \(t-1\). The recursion for \(q_t\) gives \(q_t(j)=b_j(y_t)\sum_i C_{t-1}\alpha_{t-1}(i)a_{ij}=C_{t-1}u_t(j)=C_{t-1}c_t\alpha_t(j)=C_t\alpha_t(j)\). Summing over \(j\) yields \(L_T=\sum_jq_T(j)=C_T\sum_j\alpha_T(j)=C_T\), since the normalized weights sum to one. Taking logarithms of the positive product gives the log-likelihood identity. \(\square\)

If \(c_t=0\) at any step, then every \(u_t(j)=0\), so the observed sequence has zero likelihood under the model. The normalized vector \(\alpha_t\) is then undefined; do not divide by zero. Report log likelihood \(-\infty\) and handle the zero-likelihood case explicitly.

Log Probabilities and Log-Sum-Exp

Scaling keeps the weights at a manageable magnitude. A second strategy is to store logarithms and replace products by sums. A sum of positive quantities still requires care: adding their logarithms is not valid. For real inputs \(x_1,\ldots,x_n\), the log-sum-exp function is

Directly exponentiating very negative inputs can underflow. The stable form subtracts the largest input before exponentiating.

Proof. Since \(M\) is finite, \(e^{x_r}=e^M e^{x_r-M}\) for every \(r\). Therefore \(\sum_{r=1}^{n}e^{x_r}=e^M\sum_{r=1}^{n}e^{x_r-M}\). Taking the logarithm and using \(\log(ab)=\log a+\log b\) for positive \(a,b\) gives the stated identity. Each shifted exponent has \(x_r-M\le0\), so none exceeds one. \(\square\)

In a log-domain HMM calculation, set \(\log 0=-\infty\). Products of probabilities become sums of their logarithms, while sums over states become log-sum-exp operations. For example, the update for a state \(j\) uses \(\log b_j(y_t)+\operatorname{LSE}_{i\in S}(\log\alpha_{t-1}(i)+\log a_{ij})\). After computing these unnormalized log weights, use log-sum-exp across \(j\) to obtain \(\log c_t\), then subtract \(\log c_t\) from each finite log weight to normalize. If all terms in a log-sum-exp are \(-\infty\), its value is \(-\infty\); that update has zero total weight and cannot be normalized.

Common Mistakes and Practical Checks

A frequent mistake with log-sum-exp is to subtract the maximum but forget to add it back. In the preceding example, this error reports \(\log(\tfrac13+1)=\log(\tfrac43)\) instead of \(-1000\). The correction is to use \(M+\log\sum_r e^{x_r-M}\), including the \(M\) term. Another mistake is to multiply scaling factors but then also multiply by the sum of the final normalized weights as though it contained additional likelihood information. That sum is exactly one, so the likelihood is already \(\prod_t c_t\).

For a reliable implementation, check that transition rows sum to one, each normalization factor is nonnegative, and every positive normalization produces weights summing to one. Compare the scaled likelihood with direct path enumeration on a short sequence, as in the first example. When working in logs, zero probabilities must remain \(-\infty\); replacing them by an arbitrary finite number changes the model. If observations have densities, a density value and its scaling factor can exceed one. The normalization and product identities still apply, but the factors should not automatically be interpreted as probabilities bounded by one.