ML Math Derivations (15): Hidden Markov Models

You hear footsteps behind you in the fog. You can’t see the walker, only the sounds. From the rhythm and pitch — short, soft, hurried — can you guess whether they are walking, running, or limping? And if you observed an entire sequence, which gait sequence is most likely? How likely is any sequence of sounds under your model of how walking works?

These are the three problems of HMMs, and the surprise is that all three reduce to one trick: write the joint $P(\mathbf{O}, \mathbf{I})$ as a product of local factors along time, then share sub-computations across time with dynamic programming. Brute force costs $O(N^T)$ . Forward-Backward, Viterbi, and Baum-Welch all cost $O(N^2 T)$ . The exponent collapses because the Markov assumption makes the future conditionally independent of the past given the present.

What You Will Learn#

• The HMM joint distribution and its two conditional-independence assumptions
• Forward / Backward: $P(\mathbf{O}\mid\lambda)$ and posterior smoothing $\gamma_t(i), \xi_t(i,j)$
• Viterbi: MAP decoding via DP, with the single change of sum -> max
• Baum-Welch: EM specialised to HMMs; why the M-step formulae are expected counts
• Numerical pitfalls (underflow, scaling) and modern descendants (CRF, RNN, CTC)

Prerequisites#

• Markov chains and stochastic matrices
• Dynamic programming / memoisation
• The EM algorithm and the ELBO (Part 13 )

The Model#

An HMM is the simplest non-trivial latent-variable model for sequences. Two parallel chains run in time:

• a hidden state chain $\mathbf{I}=(i_1,\dots,i_T)$ with $i_t\in\{1,\dots,N\}$ , and
• an observed emission chain $\mathbf{O}=(o_1,\dots,o_T)$ with $o_t\in\{v_1,\dots,v_M\}$ .

The model is fully described by three parameter blocks $\lambda=(\boldsymbol{\pi}, \mathbf{A}, \mathbf{B})$ :

Two conditional-independence assumptions make everything tractable:

1. First-order Markov on states — $P(i_{t+1}\mid i_{1:t}) = P(i_{t+1}\mid i_t)$ .
2. Observation independence — $o_t$ depends only on $i_t$ , not on other states or observations.

Every algorithm in this article is a clever way to sum or maximise this product without enumerating the $N^T$ hidden paths.

The 3-State Weather Toy#

Throughout the article we lean on a 3-state weather example — Sunny, Rainy, Cloudy — with the transitions shown above. Sunny is sticky ($a_{\text{SS}}=0.70$ ), and Cloudy is the hub that connects everything else. This is small enough to inspect by hand and rich enough to expose every algorithmic subtlety.

Three Problems, One Joint#

Given $\lambda$ , three questions exhaust what we can ask:

A naive approach to (1) sums the joint over all $N^T$ hidden sequences — already $\approx 10^{170}$ for $N{=}50, T{=}100$ . The next three sections make all three problems polynomial.

Forward Algorithm: Evaluation by Left-to-Right Sweep#

It is the joint probability of having generated the observations seen so far and ending up in state $i$ at time $t$ .

The bracketed sum is the only computation that crosses the time boundary; everything else is local. This is dynamic programming in the cleanest possible form.

Cost. Each of $T$ steps has $N$ targets, each requires summing over $N$ predecessors: $O(N^2 T)$ . For $N{=}50, T{=}100$ this is $2.5\times 10^{5}$ operations — about $10^{165}\times$ faster than the brute-force sum.

Underflow. Probabilities multiply geometrically, so $\alpha_t$ eventually underflows. Two standard fixes:

• Log-space with logsumexp: $\log\alpha_t(j) = \mathrm{lse}_i\!\big(\log\alpha_{t-1}(i) + \log a_{ij}\big) + \log b_j(o_t)$ .
• Scaled forward: at each $t$ divide $\alpha_t$ by $c_t = \sum_j \alpha_t(j)$ and accumulate $\log P(\mathbf{O}) = \sum_t \log c_t$ .

Backward Algorithm and Posterior Smoothing#

The forward sweep alone evaluates $P(\mathbf{O})$ . To compute the posterior over hidden states at every $t$ , we also need a right-to-left sweep.

Note the conditioning: $\beta$ is a conditional probability of the future, while $\alpha$ is a joint probability of the past. They mirror each other.

Boundary. Past time $T$ there are no observations, so $\beta_T(i) = 1$ for all $i$ .

Two Posteriors That Drive Learning#

Once both $\alpha$ and $\beta$ are tabulated we can read off, in $O(1)$ extra work per cell:

These are the two statistics Baum-Welch needs in its E-step.

Viterbi: From Sum to Max#

(The denominator $P(\mathbf{O})$ does not depend on $\mathbf{I}$ , so it can be dropped.)

After the forward pass, terminate at $i_T^* = \arg\max_i \delta_T(i)$ and backtrack: $i_t^* = \psi_{t+1}(i_{t+1}^*)$ .

Why does swapping sum for max work? Both operators distribute over the time-factorised joint — “max-times” forms a commutative semiring, just like “sum-times”. The dynamic-programming bookkeeping is identical; only the local reduction changes.

Now everything is an addition, so underflow is impossible.

Baum-Welch: EM for HMMs#

When $\lambda$ is unknown, learn it by maximising $\log P(\mathbf{O}\mid\lambda)$ . The hidden states $\mathbf{I}$ are latent, so apply EM. The result — the Baum-Welch algorithm — predates the general EM framework by several years (Baum & Petrie, 1966) and is one of the prettiest examples of EM in the wild.

E-step: Posterior Statistics#

Given the current $\lambda^{(k)}$ , run forward-backward to obtain $\gamma_t(i)$ and $\xi_t(i,j)$ . These are the expected values, under the current model, of the indicator variables $\mathbb{1}[i_t=i]$ and $\mathbb{1}[i_t=i, i_{t+1}=j]$ .

The three terms are decoupled in $\boldsymbol{\pi}, \mathbf{A}, \mathbf{B}$ — the M-step solves three independent constrained maximisations.

M-step: Three Closed-Form Updates#

Read each ratio as expected-transitions-from-$i$ -to-$j$ over expected-departures-from-$i$ , and similarly for emissions. Identical in spirit to MLE on observed counts; only “observed” is replaced by “expected under the posterior”.

Convergence#

EM guarantees $P(\mathbf{O}\mid\lambda^{(k+1)}) \geq P(\mathbf{O}\mid\lambda^{(k)})$ — the curve in the figure is provably monotonic. The catch: the surface is multi-modal, so Baum-Welch finds a local maximum. Standard remedies are random restarts, k-means / segmental k-means initialisation, or informative priors.

Application: POS Tagging#

Part-of-speech (POS) tagging is the textbook HMM application. Hidden states are tags (PRON, VERB, ADJ, NOUN, …) and observations are words. Transitions encode grammar (a determiner is usually followed by an adjective or noun); emissions encode the lexicon (the word love is most often a verb, sometimes a noun).

For “I love natural language processing” the Viterbi path lights up as PRON / VERB / ADJ / NOUN / NOUN — correct, despite processing being lexically ambiguous, because the ADJ -> NOUN transition is highly probable.

The same engine appears in speech recognition (states = phonemes, observations = MFCC frames), bioinformatics profile HMMs (states = match/insert/delete columns), and gesture recognition — whenever a discrete latent process generates a noisy observable stream.

Numerical scaling: why naive forward-backward fails#

The forward variables $\alpha_t(i) = P(\mathbf{o}_{1:t}, q_t = s_i \mid \lambda)$ shrink geometrically: each step multiplies by transition and emission probabilities, both bounded by 1. After $T = 200$ tokens with average factor $0.1$ , $\alpha_T \sim 10^{-200}$ , well below the float32 underflow threshold $\approx 10^{-38}$ and below float64’s $\approx 10^{-308}$ for any non-trivial sequence length.

The backward pass uses the same $c_t$ for normalisation, which keeps the smoothed posterior $\gamma_t(i) = \hat\alpha_t(i)\hat\beta_t(i)$ correct without further bookkeeping.

Log-space costs an extra exp and log per step but cleanly handles emission probabilities that are already log-densities (e.g. continuous Gaussian emissions). Most production HMM libraries — hmmlearn, pomegranate — use log-space by default. Viterbi already lives in log-space because $\max$ commutes with monotone transforms, so there is no scaling subtlety there.

Cost, complexity, and when to reach for HMMs#

The three core algorithms have identical asymptotic cost $O(N^2 T)$ where $N$ is the state count and $T$ the sequence length. The quadratic in $N$ comes from the transition sum; the linear in $T$ from the sweep. EM (Baum-Welch) is one forward-backward per iteration, so a full training run is $O(N^2 T \cdot I)$ for $I$ iterations and a single sequence, or $O(N^2 \sum_s T_s \cdot I)$ for a corpus.

When does this become a problem? On a 200-token sequence with $N = 45$ POS tags, $N^2 T = 4 \times 10^5$ — instant. On speech recognition with $N \approx 5000$ HMM states (one per context-dependent triphone), $N^2 = 2.5 \times 10^7$ per token. This is exactly why production speech systems used beam search on Viterbi (prune all but the top-$K$ states per step) and why the field eventually moved to RNN/Transformer encoders that avoid the $N^2$ entirely.

When should you still pick HMM over an RNN or Transformer in 2026?

1. Tiny labelled data, strong structural prior. When you have 50 sentences and 10 hidden states, the inductive bias of “discrete state at each step, Markov transition” regularises the fit far better than a neural net.
2. Exact posterior matters. Viterbi gives the global MAP path; beam search on a Transformer does not. Bioinformatics (gene finding, protein secondary structure) still uses HMMs for this reason.
3. Streaming with bounded memory. Forward filtering is $O(N)$ memory regardless of $T$ . Transformer attention is $O(T)$ .
4. Interpretability. You can read off $A_{ij}$ as “probability of moving from state $i$ to $j$ ” and validate it against domain knowledge. A Transformer’s attention head means much less.

For everything else — long sequences, large labelled corpora, complex emission distributions — modern sequence models dominate. HMMs are the right teaching example because they cleanly separate what to compute (filter, smooth, decode, learn) from how to compute it.

FAQ#

Forward vs. Viterbi — why does swapping operators matter?#

Forward returns the marginal $P(\mathbf{O}\mid\lambda) = \sum_{\mathbf{I}} P(\mathbf{O},\mathbf{I})$ ; Viterbi returns $\max_{\mathbf{I}} P(\mathbf{O},\mathbf{I})$ . Same DP skeleton, different semiring (sum-product vs. max-product). Forward answers “how plausible is this evidence?”; Viterbi answers “what story best explains it?”

Why does Viterbi maximise the joint instead of the posterior?#

Because $P(\mathbf{I}\mid\mathbf{O}) = P(\mathbf{O},\mathbf{I})/P(\mathbf{O})$ and the denominator is a constant in $\mathbf{I}$ . Maximising the joint is therefore equivalent and avoids one normalisation.

When does Baum-Welch fail?#

Three classic failure modes: (a) bad initialisation — it lands in a flat or trivial local optimum; (b) label switching — states are only identified up to permutation; (c) observation collapse — a state’s emission concentrates on observed symbols only, leaving zero probability for unseen ones. Smooth with a small Dirichlet prior to fix (c).

Why CRFs over HMMs for sequence labelling?#

CRFs are discriminative: they model $P(\mathbf{I}\mid\mathbf{O})$ directly and can exploit overlapping global features of $\mathbf{O}$ (capitalisation, suffix templates, surrounding words) without the conditional-independence straitjacket. HMMs are still preferred when you need to generate sequences or when training data is scarce.

Are HMMs obsolete in the deep-learning era?#

As stand-alone end-to-end models, mostly yes — RNN/Transformer encoders dominate. But the inference algorithms live on. CTC decoding in modern speech systems is essentially Forward over an alignment lattice; sequence-level distillation uses Viterbi; structured-output Transformers borrow from CRFs which borrow from HMMs.

How do I choose $N$ ?#

Start small, then use information criteria ($\text{AIC} = -2\log L + 2|\lambda|$ , $\text{BIC} = -2\log L + |\lambda|\log T$ ) or held-out likelihood. Bayesian non-parametrics (the iHMM with a hierarchical Dirichlet process prior) place a prior over $N$ and let the data decide.

How do I handle continuous observations?#

Replace the emission matrix with a density: a single Gaussian, a Gaussian mixture (the classical GMM-HMM in speech), or a neural density (Mixture Density Network -> “neural HMM”). The forward-backward recursions are identical; only $b_j(o_t)$ becomes a likelihood evaluation.

Exercises#

E1. With $N=2$ , $\boldsymbol{\pi}=(0.6, 0.4)$ , $\mathbf{B} = \begin{pmatrix}0.5 & 0.5 \\ 0.4 & 0.6\end{pmatrix}$ , and $o_1 = v_1$ , compute $\alpha_1$ .

E2. For $N{=}50, T{=}100$ , compare Forward to brute-force enumeration.

E3. Interpret $\hat a_{ij} = \frac{\sum_t \xi_t(i,j)}{\sum_t \gamma_t(i)}$ .

E4. Show that running Forward then Backward gives the same $P(\mathbf{O}\mid\lambda)$ as Forward alone.

E5. Explain why Viterbi requires only $O(N^2 T)$ time but $O(NT)$ memory for back-pointers.

What’s next#

HMMs are generative — they model $P(x, z)$ and invert to $P(z\mid x)$ . On supervised sequence labelling (named-entity recognition, POS tagging) this wastes effort: I only care about $P(z\mid x)$ , but I am forced to model $P(x)$ as well, often under feature-independence and Markov assumptions that conflict with real data. The next chapter is the discriminative answer: conditional random fields (CRFs).

CRFs model $P(z\mid x)$ directly, make no distributional assumption on $x$ , and therefore allow arbitrary feature functions (windowed context, word shape, gazetteers). The training objective is the conditional log-likelihood, convex; inference is still forward-backward and Viterbi, almost identical dynamic programming to HMMs. The HMM-to-CRF transition is the same “generative → discriminative” shift the supervised line went through, replayed on sequences — and it is the key to understanding what the CRF head is actually doing inside BiLSTM-CRF, Transformer-CRF, and similar deep architectures.

References#

1. Rabiner, L. R. (1989). A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition. Proc. IEEE 77(2), 257-286.
2. Baum, L. E., & Petrie, T. (1966). Statistical Inference for Probabilistic Functions of Finite State Markov Chains. Annals of Math. Statist. 37(6).
3. Viterbi, A. J. (1967). Error Bounds for Convolutional Codes and an Asymptotically Optimum Decoding Algorithm. IEEE Trans. IT 13(2).
4. Murphy, K. P. (2012). Machine Learning: A Probabilistic Perspective, Ch. 17. MIT Press.
5. Bishop, C. M. (2006). Pattern Recognition and Machine Learning, Ch. 13. Springer.
6. Jelinek, F. (1997). Statistical Methods for Speech Recognition. MIT Press.
7. Lafferty, J., McCallum, A., & Pereira, F. (2001). Conditional Random Fields. ICML.
8. Eddy, S. R. (1998). Profile Hidden Markov Models. Bioinformatics 14(9).