HMM Algorithms: Forward, Backward, Viterbi, Baum-Welch

Source: Unit 3 §13

#NumbersThe canonical exam HMM, used on every algorithm below

• States: S1, S2. Initial: $\pi_1 = 0.2$, $\pi_2 = 0.8$.
• Transitions: $a_{11} = 0.4$, $a_{12} = 0.6$, $a_{21} = 0.3$, $a_{22} = 0.7$.
• Emissions: $b_1(V_1) = 0.1$, $b_1(V_2) = 0.4$, $b_1(V_3) = 0.5$; $b_2(V_1) = 0.3$, $b_2(V_2) = 0.5$, $b_2(V_3) = 0.2$.
• Observation sequence: O = {V1, V3, V2}.

Forward algorithm: problem 1, likelihood

The forward probability $\alpha_i(t) = P(O_1 O_2 \dots O_t, X_t = s_i \mid \lambda)$ is the probability of the observations up to time $t$ and being in state $i$ at time $t$.

⚠GotchaWhy not brute force

Enumerating all state sequences is exponential in the length of the sequence. The forward algorithm uses dynamic programming - it solves time $t$ using the result already computed for $t-1$ - which makes it linear in sequence length.

$$\alpha_1(i) = \pi_i \cdot b_i(O_1)$$ $$\alpha_{t+1}(j) = \left[ \sum_i \alpha_t(i) \cdot a_{ij} \right] \cdot b_j(O_{t+1})$$ $$P(O \mid \lambda) = \sum_i \alpha_T(i)$$

The three lines are initialisation, recursion and termination. Termination is just "sum the last column".

→StepsThe alpha table, computed

1. Column 1, V1 observed. $\alpha_1(S_1) = \pi_1 b_1(V_1) = 0.2 \times 0.1 = \mathbf{0.02}$ and $\alpha_1(S_2) = \pi_2 b_2(V_1) = 0.8 \times 0.3 = \mathbf{0.24}$.
2. Column 2, V3 observed. $\alpha_2(S_1) = (0.02 \times 0.4 + 0.24 \times 0.3) \times 0.5 = 0.08 \times 0.5 = \mathbf{0.04}$.
3. $\alpha_2(S_2) = (0.02 \times 0.6 + 0.24 \times 0.7) \times 0.2 = 0.18 \times 0.2 = \mathbf{0.036}$.
4. Column 3, V2 observed. $\alpha_3(S_1) = (0.04 \times 0.4 + 0.036 \times 0.3) \times 0.4 = 0.0268 \times 0.4 = \mathbf{0.01072}$.
5. $\alpha_3(S_2) = (0.04 \times 0.6 + 0.036 \times 0.7) \times 0.5 = 0.0492 \times 0.5 = \mathbf{0.0246}$.
6. Terminate: $P(O \mid \lambda) = 0.01072 + 0.0246 = \mathbf{0.03532}$.

Alpha	t=1 (V1)	t=2 (V3)	t=3 (V2)
S1	0.02	0.04	0.01072
S2	0.24	0.036	0.0246

Every node holds the total probability of all paths that reach it, so each column is computed once and reused. That reuse is what turns an exponential enumeration into a linear scan.

Backward algorithm

The backward probability $\beta_i(t)$ is the probability of seeing the observations from time $t+1$ to the end, given that we are in state $i$ at time $t$.

$$\beta_T(i) = 1 \quad \text{for all } i$$ $$\beta_t(i) = \sum_j a_{ij} \cdot b_j(O_{t+1}) \cdot \beta_{t+1}(j)$$

⚠GotchaBeta fills right to left, and its emission index is t+1

The last column is all 1s, not all 0s - there are no remaining observations to explain, and an empty product is 1. Note also that $\beta$ uses $b_j(O_{t+1})$, the next observation, whereas $\alpha$ uses the current one. Getting that index wrong is the usual source of a wrong beta table.

→StepsThe beta table, computed right to left

1. Column 3, initialisation. $\beta_3(S_1) = 1$ and $\beta_3(S_2) = 1$.
2. Column 2, with V2 as $O_3$. $\beta_2(S_1) = 0.4 \times 0.4 + 0.6 \times 0.5 = 0.16 + 0.30 = \mathbf{0.46}$.
3. $\beta_2(S_2) = 0.3 \times 0.4 + 0.7 \times 0.5 = 0.12 + 0.35 = \mathbf{0.47}$.
4. Column 1, with V3 as $O_2$. $\beta_1(S_1) = 0.4 \times 0.5 \times 0.46 + 0.6 \times 0.2 \times 0.47 = 0.092 + 0.0564 = \mathbf{0.1484}$.
5. $\beta_1(S_2) = 0.3 \times 0.5 \times 0.46 + 0.7 \times 0.2 \times 0.47 = 0.069 + 0.0658 = \mathbf{0.1348}$.

Beta	t=1 (V1)	t=2 (V3)	t=3 (V2)
S1	0.1484	0.46	1
S2	0.1348	0.47	1

Combining alpha and beta

$P(O \mid \lambda)$ can be recovered from any column as the dot product of that column's $\alpha$ and $\beta$ values, and every column gives the same answer. That is the standard way to check your arithmetic.

→StepsThree columns, one answer

1. Column 1: $0.02 \times 0.1484 + 0.24 \times 0.1348 = 0.03532$.
2. Column 2: $0.04 \times 0.46 + 0.036 \times 0.47 = 0.03532$.
3. Column 3: $0.01072 \times 1 + 0.0246 \times 1 = 0.03532$.

The occupation or gamma probability is the probability of being in state $i$ at time $t$ given all the observations, past and future:

$$\gamma_t(i) = \frac{\alpha_t(i) \cdot \beta_t(i)}{P(O \mid \lambda)}$$

→StepsThe gamma table, dividing through by 0.03532

1. $\gamma_1(S_1) = 0.02 \times 0.1484 / 0.03532 = \mathbf{0.08403}$ and $\gamma_1(S_2) = 0.24 \times 0.1348 / 0.03532 = \mathbf{0.91596}$.
2. $\gamma_2(S_1) = 0.04 \times 0.46 / 0.03532 = \mathbf{0.52095}$ and $\gamma_2(S_2) = 0.036 \times 0.47 / 0.03532 = \mathbf{0.47904}$.
3. $\gamma_3(S_1) = 0.01072 \times 1 / 0.03532 = \mathbf{0.30351}$ and $\gamma_3(S_2) = 0.0246 \times 1 / 0.03532 = \mathbf{0.69648}$.

Gamma	t=1	t=2	t=3
S1	0.08403	0.52095	0.30351
S2	0.91596	0.47904	0.69648

◎Exam cueEach gamma column sums to 1

It has to: at any time $t$ the chain is in exactly one state, so the probabilities over states must total 1. Picking the highest-gamma state at each time is called posterior decoding, and it is a different answer from Viterbi.

Viterbi algorithm: problem 2, decoding

Decoding means finding the single most likely hidden state sequence given the observations.

Viterbi is almost identical to the forward algorithm, with the sum replaced by a max:

$$v_t(j) = \max_i \left[ v_{t-1}(i) \cdot a_{ij} \right] \cdot b_j(O_t)$$

⇄CompareThe one-operator difference

Forwardαₜ(j) = [ Σᵢ αₜ₋₁(i)·aᵢⱼ ] · bⱼ(Oₜ). Sums over all paths into the node, so each node holds a TOTAL probability. Answers "how likely is this observation sequence".

Viterbivₜ(j) = [ maxᵢ vₜ₋₁(i)·aᵢⱼ ] · bⱼ(Oₜ). Keeps only the most likely path into the node, and records a BACKPOINTER to the max-giving predecessor. Answers "which states produced it".

Viterbi also records the back-pointer - the predecessor that gave the max - at every node, and backtracks from the best final node at the end to recover the optimal state sequence.

→StepsViterbi on the canonical HMM

1. t = 1, V1. $v_1(S_1) = 0.2 \times 0.1 = \mathbf{0.02}$, $v_1(S_2) = 0.8 \times 0.3 = \mathbf{0.24}$. No backpointers yet.
2. t = 2 into S1, V3. Candidates $0.02 \times 0.4 = 0.008$ and $0.24 \times 0.3 = 0.072$. Max is 0.072 from S2, so $v_2(S_1) = 0.072 \times 0.5 = \mathbf{0.036}$.
3. t = 2 into S2. Candidates $0.02 \times 0.6 = 0.012$ and $0.24 \times 0.7 = 0.168$. Max is 0.168 from S2, so $v_2(S_2) = 0.168 \times 0.2 = \mathbf{0.0336}$.
4. t = 3 into S1, V2. Candidates $0.036 \times 0.4 = 0.0144$ and $0.0336 \times 0.3 = 0.01008$. Max is 0.0144 from S1, so $v_3(S_1) = 0.0144 \times 0.4 = \mathbf{0.00576}$.
5. t = 3 into S2. Candidates $0.036 \times 0.6 = 0.0216$ and $0.0336 \times 0.7 = 0.02352$. Max is 0.02352 from S2, so $v_3(S_2) = 0.02352 \times 0.5 = \mathbf{0.01176}$.
6. Terminate and backtrack. The larger final value is $0.01176$ at S2; its backpointer is S2 at t=2, whose backpointer is S2 at t=1. The best sequence is S2 → S2 → S2.

Same trellis, one operator changed. Only the winning edge into each node survives, and following the bold edges back from the best final node recovers the optimal state sequence S2, S2, S2.

⇄CompareLikelihood versus Viterbi, the dishonest-casino intuition

Likelihood (γ)At each time, sum over sequences to get the state distribution at time t. It fluctuates a lot from one step to the next - a narrow, local view.

ViterbiThe single most likely path through the whole trellis. Smoother, and it reflects the transition probabilities properly - a global view.

Baum-Welch algorithm: problem 3, learning

Baum-Welch learns the HMM parameters $\lambda = (A, B, \pi)$ from observations alone. It is EM applied to HMMs.

≡FactsThe chicken and egg, and the extra quantity it needs

• If we knew the state-occupation probability $\gamma$, we could derive the emission and transition probabilities. If we knew those, we could derive $\gamma$. Neither is available first.
• $\xi$ (xi) is the probability of transitioning from state $i$ to state $j$ after time $t$, given all the observations. Like $\gamma$, it is computed from $\alpha$ and $\beta$.

→StepsThe EM loop for an HMM

1. E-step: with the current model, run forward-backward to compute $\gamma$ and $\xi$ - that is, establish $p(\gamma, \xi \mid x, a, b)$.
2. M-step: re-estimate the best transition $a$ and emission $b$ as point estimates that maximise the objective.
3. Repeat: fix one set of parameters to improve the other, until convergence.

Baum-Welch is EM wearing HMM notation. Forward-backward supplies the soft counts, and the soft counts supply the next model.

≡FactsApplications of HMMs

• Speech recognition.
• DNA and base-pair analysis.
• Part-of-speech tagging.
• The dishonest casino - inferring which die, fair or biased, is being used from the observed rolls.