Viterbi algorithm

Why you need this

A designer hands you a level built by hand. You want to know which column was a Run, which a Gap, which a Climb, so you can count how sections follow each other and seed your generator with real numbers. Nobody labeled the columns. Viterbi labels them for you: given the tiles and a hidden Markov model, it returns the one hidden path most likely to have produced them.

The idea

The decoding question is: which hidden sequence $X_1, \dots, X_T$ maximizes $P(X_1..X_T, O_1..O_T)$ for the observations you have? There are $N^T$ candidates, and as with the forward algorithm, dynamic programming collapses them into a table of $T \times N$ cells.

Define the Viterbi variable $\delta_t(j)$ as the best log-probability of any hidden path that ends in state $j$ at time $t$ and emits $O_1..O_t$ along the way. Logs, because a path's probability is a product of $2T$ factors that would underflow, and because in log space products become sums (see logarithms and underflow). Since $\log$ is increasing, the path with the largest log-probability is the path with the largest probability.

Initialization.

$$\delta_1(j) = \log \pi_j + \log B_{j O_1}.$$

Recursion. The best path ending in $j$ at $t+1$ came from the best path ending in some $i$ at $t$, then took the step $i \to j$ and emitted $O_{t+1}$:

$$\delta_{t+1}(j) = \max_{i} \big[ \delta_t(i) + \log A_{ij} \big] + \log B_{j O_{t+1}}.$$

Compare this with the forward recursion, $\alpha_{t+1}(j) = [\sum_i \alpha_t(i) A_{ij}] B_{j O_{t+1}}$. Same shape, same table, with $\max$ where the sum was and logs where the products were. Forward asks "how much probability flows into $j$ from all predecessors"; Viterbi asks "which single predecessor sends the most".

Backpointers. Record which $i$ won each max:

$$\psi_{t+1}(j) = \arg\max_i \big[ \delta_t(i) + \log A_{ij} \big].$$

Termination. The best final state is $X_T^* = \arg\max_j \delta_T(j)$, and the best path's log-probability is $\delta_T(X_T^*)$.

Backtracking. Walk the pointers from the end: $X_t^* = \psi_{t+1}(X_{t+1}^*)$ for $t = T-1$ down to $1$.

In log space with signs flipped, $-\log A_{ij}$ and $-\log B_{jk}$ are non-negative costs, and Viterbi is a shortest-path search through a graph with $T$ layers of $N$ nodes: the same relaxation you would write for a layered DAG, and the reason it is exact rather than a heuristic.

Worked example

Two hidden states, 1 Run and 2 Gap; two tiles, 1 flat and 2 pit; the same model as the forward page:

$$\pi = (0.8, 0.2), \qquad A = \begin{pmatrix} 0.7 & 0.3 \\ 0.6 & 0.4 \end{pmatrix}, \qquad B = \begin{pmatrix} 0.9 & 0.1 \\ 0.2 & 0.8 \end{pmatrix}.$$

Observations: $O = (\text{flat}, \text{pit}, \text{pit}, \text{flat})$, so $T = 4$.

The code works in logs. For hand checking it is easier to carry the raw product and take the log at the end; both columns are shown. Each cell lists the two candidates $\delta_{t-1}(i) \cdot A_{ij}$ (in raw form), the winner times the emission, and the backpointer.

	$t = 1$ (flat)	$t = 2$ (pit)	$t = 3$ (pit)	$t = 4$ (flat)
Run candidates		from Run $0.72 \cdot 0.7 = 0.504$; from Gap $0.04 \cdot 0.6 = 0.024$	from Run $0.0504 \cdot 0.7 = 0.03528$; from Gap $0.1728 \cdot 0.6 = 0.10368$	from Run $0.010368 \cdot 0.7 = 0.0072576$; from Gap $0.055296 \cdot 0.6 = 0.0331776$
$\delta_t(\text{Run})$ raw	$0.8 \cdot 0.9 = 0.72$	$0.504 \cdot 0.1 = 0.0504$	$0.10368 \cdot 0.1 = 0.010368$	$0.0331776 \cdot 0.9 = 0.0298598$
$\delta_t(\text{Run})$ log	$-0.329$	$-2.988$	$-4.569$	$-3.511$
$\psi_t(\text{Run})$		Run	Gap	Gap
Gap candidates		from Run $0.72 \cdot 0.3 = 0.216$; from Gap $0.04 \cdot 0.4 = 0.016$	from Run $0.0504 \cdot 0.3 = 0.01512$; from Gap $0.1728 \cdot 0.4 = 0.06912$	from Run $0.010368 \cdot 0.3 = 0.0031104$; from Gap $0.055296 \cdot 0.4 = 0.0221184$
$\delta_t(\text{Gap})$ raw	$0.2 \cdot 0.2 = 0.04$	$0.216 \cdot 0.8 = 0.1728$	$0.06912 \cdot 0.8 = 0.055296$	$0.0221184 \cdot 0.2 = 0.0044237$
$\delta_t(\text{Gap})$ log	$-3.219$	$-1.756$	$-2.895$	$-5.421$
$\psi_t(\text{Gap})$		Run	Gap	Gap

Termination: at $t = 4$, Run ($-3.511$) beats Gap ($-5.421$), so $X_4^* = \text{Run}$.

Backtracking: $\psi_4(\text{Run}) = \text{Gap}$, so $X_3^* = \text{Gap}$. $\psi_3(\text{Gap}) = \text{Gap}$, so $X_2^* = \text{Gap}$. $\psi_2(\text{Gap}) = \text{Run}$, so $X_1^* = \text{Run}$.

The decoded path is Run, Gap, Gap, Run, with probability $0.0298598$ (log $-3.511$). Read it against the tiles: flat under Run, two pits under two Gaps, flat under Run. Note that at $t = 2$ the forward algorithm's filtered belief and Viterbi's choice agree (Gap), but they need not in general; Viterbi commits to a whole path, and the filtered column does not.

Note also that $\delta_t(\text{Run})$ at $t = 4$ is larger than at $t = 3$. That is fine: the emission at $t = 4$ is flat under Run (0.9), much larger than pit under Run (0.1) at $t = 3$, and the switch of predecessor from Run to Gap brings a bigger $\delta_3$ along.

JavaScript

// pi[i], A[i][j], B[j][k]: probabilities. O: symbol indices (0-based).
// Returns the most likely state path and its log-probability.
function viterbi(pi, A, B, O) {
  const N = pi.length, T = O.length;
  const lg = (x) => (x > 0 ? Math.log(x) : -Infinity);
  const delta = [], psi = [];
  delta[0] = pi.map((p, j) => lg(p) + lg(B[j][O[0]]));
  psi[0] = new Array(N).fill(-1);
  for (let t = 1; t < T; t++) {
    delta[t] = new Array(N); psi[t] = new Array(N);
    for (let j = 0; j < N; j++) {
      let best = -Infinity, arg = 0;
      for (let i = 0; i < N; i++) {
        const v = delta[t - 1][i] + lg(A[i][j]);
        if (v > best) { best = v; arg = i; }
      }
      delta[t][j] = best + lg(B[j][O[t]]);
      psi[t][j] = arg;
    }
  }
  let last = 0;
  for (let j = 1; j < N; j++) if (delta[T - 1][j] > delta[T - 1][last]) last = j;
  const path = new Array(T);
  path[T - 1] = last;
  for (let t = T - 1; t > 0; t--) path[t - 1] = psi[t][path[t]];
  return { path, logP: delta[T - 1][last] };
}
const pi = [0.8, 0.2], A = [[0.7, 0.3], [0.6, 0.4]], B = [[0.9, 0.1], [0.2, 0.8]];
console.log(viterbi(pi, A, B, [0, 1, 1, 0]));   // path [0, 1, 1, 0], logP -3.511...

In a map generator

• Labeling approved levels. Take a hand-authored or approved level, turn each column into a tile symbol, and decode. Each column now carries a section label.
• Counting to initialize. With labels in hand, count how often each label follows each label and normalize the rows: that is an estimate of $A$. Count how often each tile appears under each label: that is an estimate of $B$. Add 1 to every count first so nothing is zero. This is cheap and gives Baum-Welch a starting point near a sensible answer instead of a random one.
• Segmenting for the editor. Show the decoded labels as a colored strip above the level in the design tool. Designers can correct labels, and the corrections become training data.
• Explaining a rejected level. When the forward score flags a generated level as unlikely, the Viterbi path tells you where it went wrong: the run of six Hazards at columns 40 to 45.

Common mistakes

• Using probabilities instead of logs. Fine at $T = 4$; at $T = 200$ every cell is 0 and every backpointer is 0. Log from the first line.
• $\log 0$. A zero in $A$ or $B$ must become $-\infty$, not NaN. Math.log(0) is $-\infty$ in JavaScript, which is fine, but a floored value of 1e-300 that then gets multiplied is not. Keep the guard explicit.
• Storing only the last backpointer column. Backtracking needs all $T$ columns of $\psi$. Memory is $O(T N)$ and that is unavoidable for the path.
• Reading the per-column argmax as the path. The state that is individually most likely at each $t$ can string together into a path that is impossible (a transition with $A_{ij} = 0$). Viterbi is the only thing that returns a coherent path.
• Ties broken inconsistently. Two candidates with equal log-probability must resolve the same way every run or the labels flicker between builds. Use strict > and iterate in a fixed order.

Cost

Time is $O(T \cdot N^2)$: $T$ observations, $N$ cells per column, $N$ candidates per cell. Memory is $O(T \cdot N)$ for the backpointers (the $\delta$ table itself can be two columns). For $N = 6$ and a 200-column level that is 7,200 comparisons and 1,200 stored bytes, far below anything that matters. It hurts only when $N$ grows into the hundreds, which happens if you encode a second-order model as pairs of states.

Going further

• Baum-Welch, which replaces hard Viterbi labels with soft expected counts.
• The forward-backward algorithm and posterior decoding, for the most likely state at each time rather than the most likely whole path.
• Beam search, an approximation that keeps only the top few $\delta_t$ entries when $N$ is large.
• Dijkstra's algorithm, to see Viterbi as shortest path on a layered graph.