Viterbi Algorithm: Dynamic Programming for Optimal Sequence Decoding

When labeling a sequence of tokens, you face a fundamental decision: should you choose each label independently, or should you consider the entire sequence at once? Greedy decoding picks the best label at each position, ignoring how earlier choices affect later ones. But sequence labeling isn't about finding locally optimal tags; it's about finding the globally optimal tag sequence.

The Viterbi algorithm solves this problem elegantly. Given a sequence of observations and a model that defines transition and emission probabilities, Viterbi finds the single best sequence of hidden states. For part-of-speech tagging, it finds the best tag sequence. For named entity recognition, it finds the best BIO label sequence. The algorithm runs in time linear in the sequence length, making it practical for real applications.

This chapter develops the Viterbi algorithm from first principles. You'll understand why exhaustive search is intractable, how dynamic programming makes the problem solvable, and how to implement efficient decoding in practice. We also explore log-space computation to avoid numerical underflow and show how Viterbi forms the foundation for beam search decoding in neural sequence models.

Consider a hidden Markov model for part-of-speech tagging. You observe a sequence of words $w_1, w_2, \ldots, w_n$ and want to find the sequence of tags $t_1, t_2, \ldots, t_n$ that maximizes the joint probability of seeing both the tags and the words together. The HMM factors this joint probability into a product of simpler terms:

where:

• $t_1, \ldots, t_n$: the sequence of hidden states (POS tags) we want to infer
• $w_1, \ldots, w_n$: the observed sequence of words
• $P(t_1)$: the initial state probability, capturing how likely the sentence starts with tag $t_1$
• $P(t_i | t_{i-1})$: the transition probability from tag $t_{i-1}$ to tag $t_i$, encoding grammatical patterns like "determiners often precede nouns"
• $P(w_i | t_i)$: the emission probability of word $w_i$ given tag $t_i$, capturing which words are typical for each tag

The first product captures transition probabilities: how likely is each tag given the previous one? The second product captures emission probabilities: how likely is each word given its tag? This factorization embodies the Markov assumption: each tag depends only on its immediate predecessor, not on the full history.

A greedy approach picks the best tag at each position independently, considering only the previous tag. At each position $i$, greedy decoding selects:

where:

• $t_i^*$: the greedily chosen tag at position $i$
• $\arg\max_{t_i}$: selects the tag $t_i$ that maximizes the following expression
• $P(w_i | t_i)$: the emission probability of the current word given candidate tag $t_i$
• $P(t_i | t_{i-1}^*)$: the transition probability from the previously chosen tag $t_{i-1}^*$ to candidate tag $t_i$

This seems reasonable but can produce suboptimal sequences. The problem is that a locally good choice might force bad choices later. Consider tagging "The can can hold water." A greedy tagger might assign MODAL to the first "can" because modals are common after determiners. But this makes the second "can" difficult: modals rarely follow modals. The globally optimal solution assigns NOUN to the first "can" and VERB to the second, even though NOUN isn't the highest-probability tag for "can" in isolation.

This example demonstrates why greedy decoding fails. The locally optimal choice of MD for the first "can" yields score 0.180, beating NN's 0.125. But this decision creates a problem: the second "can" must follow a modal, and modals rarely follow modals (P(MD|MD) = 0.05). The alternative path through NN enables a much higher score at the next position (0.060 vs 0.030), potentially yielding a better overall sequence probability.

Let's trace the cumulative scores for both paths:

The greedy path starts stronger (0.180 vs 0.125) because MD has higher emission probability for "can." But this local advantage creates a problem: the MD→VB transition at position 3 yields only 0.0216, while NN→VB would have enabled a score of 0.0075. The globally optimal path requires looking ahead, something greedy decoding cannot do.

If greedy decoding doesn't work, why not try all possible tag sequences? With $k$ possible tags and a sequence of length $n$, there are $k^n$ possible sequences. For a modest vocabulary of 45 tags and a 20-word sentence, that's $45^{20} \approx 10^{33}$ sequences, far more than the number of atoms in the observable universe.

Even with a modest 45-tag vocabulary, a 20-word sentence produces more possible sequences than atoms in the observable universe. Checking each sequence individually would take longer than the age of the universe on any computer. The key insight is that we don't need to enumerate every sequence. Many sequences share common prefixes, and we can avoid recomputing scores for these shared subsequences. This is where dynamic programming enters.

We've established that greedy decoding fails because local choices can force globally suboptimal sequences, and that exhaustive search is computationally impossible. How can we find the optimal path without examining every possibility? The answer lies in a beautiful insight from dynamic programming: we can build the optimal solution incrementally from smaller optimal solutions.

Imagine you're at position 3 in a sequence, and you want to know the best path that ends at the tag NOUN. Here's the crucial observation: regardless of how that path reached position 2, the best continuation to NOUN at position 3 depends only on:

1. What state we're in at position 2 (not how we got there)
2. The transition probability to NOUN
3. The emission probability of the observed word from NOUN

This property is called optimal substructure. If the best path through the entire sequence passes through state $s'$ at position 2 and then reaches NOUN at position 3, then the portion of that path from the start to $s'$ must itself be optimal. Why? If there were a better path to $s'$, we could substitute it and get a better overall path, contradicting our assumption that we had the best path.

This principle transforms our approach entirely. Instead of tracking complete paths (exponentially many), we only need to track the best score for reaching each state at each position. At any position, we have just $k$ values to maintain, one per possible state, rather than an exponentially growing set of paths.

Let's formalize this insight. We define $V_i(s)$ as the probability of the most likely path that:

• Ends at state $s$ at position $i$
• Generates the observations $w_1, \ldots, w_i$ seen so far

How do we compute $V_i(s)$? Consider any path ending at state $s$ at position $i$. That path must have come from some predecessor state $s'$ at position $i-1$. The probability of this path factors into three terms:

1. The probability of the best path to $s'$: This is $V_{i-1}(s')$, which we've already computed
2. The transition probability from $s'$ to $s$: This is $P(s | s')$
3. The probability of emitting the current observation from $s$: This is $P(w_i | s)$

Since we want the best path to $s$, we take the maximum over all possible predecessor states. This gives us the Viterbi recurrence:

where:

• $V_i(s)$: the probability of the best path ending at state $s$ at position $i$
• $\max_{s'}$: takes the maximum over all possible predecessor states $s'$
• $V_{i-1}(s')$: the best score for reaching state $s'$ at the previous position (already computed)
• $P(s | s')$: the transition probability from predecessor state $s'$ to current state $s$
• $P(w_i | s)$: the emission probability of observing word $w_i$ from state $s$

Notice that $P(w_i | s)$ appears inside the maximization but doesn't depend on the predecessor $s'$. This is because once we decide to be in state $s$ at position $i$, the emission probability is fixed regardless of how we arrived. We could factor it out:

This form is computationally equivalent but makes clear that we're choosing the best predecessor based solely on the accumulated path score and transition, then multiplying by the emission.

The recurrence needs a starting point. At position 1, there is no predecessor to consider. The probability of the best path ending at state $s$ after seeing just the first observation is simply:

where:

• $P(s)$: the initial state probability (how likely the sequence starts in state $s$)
• $P(w_1 | s)$: the emission probability of the first observation from state $s$

This captures both how likely we are to start in state $s$ and how well state $s$ explains the first observed word.

Let's count the operations. At each of the $n$ positions (after the first), we compute $V_i(s)$ for each of the $k$ states. Computing each $V_i(s)$ requires examining $k$ predecessors. Total: $O(n \cdot k^2)$ operations.

Compare this to the $k^n$ paths we would examine exhaustively:

The exponential has become polynomial. This is the power of dynamic programming: by recognizing that we can solve the problem incrementally and that we only need to remember the best score per state (not every path), we collapse an intractable problem into an efficient algorithm.

It helps to visualize the Viterbi computation as filling a grid called a trellis. The horizontal axis represents positions in the sequence (1 through $n$), and the vertical axis represents possible states (the $k$ tags or labels). Each cell $(i, s)$ holds the score $V_i(s)$, the probability of the best path ending at state $s$ at position $i$.

Arrows connect cells at adjacent positions. Each arrow from $(i-1, s')$ to $(i, s)$ represents one term in the maximization: "if I came from state $s'$, what would my score be?" The algorithm fills this grid column by column, left to right, computing each cell from the cells in the previous column.

In the figure above, gray arrows show all possible transitions between adjacent positions. The red path highlights the optimal sequence discovered by Viterbi: DT → NN → VB → NN for the sentence "The can hold water." Notice that:

• At position 1, DT (determiner) wins because "The" is overwhelmingly a determiner
• At position 2, NN (noun) wins, not because "can" is most often a noun, but because this choice enables better continuations
• The algorithm found this globally optimal path without explicitly enumerating all 81 possible sequences (3 states × 4 positions)

Now we understand how Viterbi computes scores efficiently. But there's a missing piece: at the end of the algorithm, we know the score of the best path, but we haven't recorded the path itself. How do we reconstruct which states we actually passed through?

The Viterbi recurrence gives us the score of the best path to each state, but we've lost something in the process: which path achieved that score? At each step, we took the maximum over predecessors, keeping only the winning score. To reconstruct the actual optimal sequence, we need to remember which predecessor won at each cell.

Think of it this way: when we compute $V_3(\text{NOUN})$, we're answering "what's the best path probability ending at NOUN at position 3?" But we're not recording "and that path came through VERB at position 2." Without this information, when we finish computing the entire trellis, we'll know that the best final state is (say) NOUN with score 0.0073, but we won't know the sequence of states that led there.

The solution is elegantly simple: alongside each score, store the predecessor that achieved it. We call this a backpointer because it points backward to the previous state in the optimal path.

The backpointer $B_i(s)$ records which predecessor state $s'$ maximized the score when computing $V_i(s)$:

where:

• $B_i(s)$: the backpointer for state $s$ at position $i$, storing which predecessor led to the best score
• $\arg\max_{s'}$: returns the state $s'$ (not the value) that maximizes the expression
• The expression inside is identical to the Viterbi recurrence

Notice the difference between $\max$ and $\arg\max$:

• $\max$ returns the value: "the best score is 0.042"
• $\arg\max$ returns the argument: "the best predecessor is VERB"

We compute both simultaneously: the same computation that finds the maximum score also identifies which predecessor achieved it.

Once we've filled the entire trellis, we have scores and backpointers for every (position, state) pair. Reconstruction proceeds in three steps:

Step 1: Find the best final state. The optimal path ends at whichever state has the highest score at the final position:

where:

• $s_n^*$: the optimal state at the final position $n$
• $V_n(s)$: the Viterbi score at the final position for each state $s$

Step 2: Follow backpointers backward. Starting from $s_n^*$, we trace back through the stored pointers:

where:

• $s_i^*$: the optimal state at position $i$
• $B_{i+1}(s_{i+1}^*)$: the backpointer stored at position $i+1$ for state $s_{i+1}^*$, which tells us which state at position $i$ led to the best path through $s_{i+1}^*$

Each backpointer tells us: "to reach this state optimally, I came from that state." By following this chain, we recover the complete path.

Step 3: Reverse the sequence. Since we traced from position $n$ back to position 1, the recovered sequence is in reverse order. A simple reversal gives us the forward path.

The beauty of this approach is that backpointers require only $O(n \cdot k)$ additional storage, one pointer per cell in the trellis, and add no complexity to the forward pass (we're already examining all predecessors to compute the maximum). The backward pass is a simple $O(n)$ traversal.

Let's walk through how this code implements the formulas we developed:

1. Initialization (lines 17-24): We create two data structures: V for Viterbi scores and backpointer for tracking predecessors. Each is a list of dictionaries, one per position.

2. Base case (lines 26-30): For the first observation, we apply $V_1(s) = P(s) \cdot P(w_1 | s)$. There's no predecessor, so backpointers are None.

3. Recursion (lines 32-50): The nested loops implement the recurrence. For each position i and each state s, we examine all predecessors s_prev, compute the score $V_{i-1}(s') \cdot P(s|s') \cdot P(w_i|s)$, and keep the maximum. We store both the score (V[i][s]) and which predecessor achieved it (backpointer[i][s]).

4. Termination and backtracking (lines 52-63): We find the state with the highest final score (this is $s_n^* = \arg\max_s V_n(s)$), then follow backpointers backward to reconstruct the path.

The implementation uses a small epsilon (1e-10) instead of zero for missing probabilities to avoid numerical issues with exact zeros.

Now let's trace through this algorithm on a concrete problem. We'll use a simple HMM for weather inference: given a sequence of observed activities (walking, shopping, cleaning), infer the most likely sequence of weather states (sunny or rainy).

The HMM parameters encode two types of knowledge about the weather-activity relationship:

Transition Probabilities $P(s_t | s_{t-1})$: How likely is each weather state given the previous day's weather?

Emission Probabilities $P(w | s)$: How likely is each activity given the current weather?

These tables reveal the model's structure. Weather tends to persist: sunny days have 70% chance of staying sunny, rainy days 60% chance of staying rainy. For activities, walking is a strong indicator of sunny weather, while cleaning suggests rain. Shopping is the most ambiguous activity, slightly favoring rainy days.

Let's interpret this step-by-step computation:

Position 0 (observing "walk"): We apply the base case formula $V_1(s) = P(s) \cdot P(w_1|s)$. Sunny scores higher (0.360 vs 0.040) for two reinforcing reasons: the prior favors sunny days (0.6 vs 0.4), and walking is much more likely on sunny days (0.6 vs 0.1).

Position 1 (observing "shop"): Now the recurrence kicks in. For each current state, we consider both possible predecessors:

• To reach Sunny: the best predecessor is Sunny (score 0.0756), since we benefit from both the high previous score and a favorable transition (sunny days tend to stay sunny)
• To reach Rainy: again the best predecessor is Sunny (score 0.0432), despite the transition penalty, because the previous Sunny score was so much higher than Rainy

The result: The optimal path Sunny → Sunny → Rainy aligns with intuition: walking strongly suggests good weather, shopping is ambiguous (slightly favoring rain), and cleaning strongly suggests rain. The algorithm found this without examining all 8 possible paths. It computed only 12 values (2 states × 3 positions × 2 predecessor checks).

Multiplying many small probabilities causes numerical underflow. For a sequence of 100 tokens with probabilities around 0.01, the path probability would be approximately $0.01^{100} = 10^{-200}$, far below the smallest representable floating-point number (approximately $10^{-308}$ for 64-bit floats).

The solution is to work in log space. The logarithm function has a useful property: it converts products into sums. Since $\log(ab) = \log(a) + \log(b)$, we can rewrite the Viterbi recurrence entirely in terms of log probabilities:

where:

• $\log V_i(s)$: the log probability of the best path ending at state $s$ at position $i$
• $\log V_{i-1}(s')$: the log probability of the best path to predecessor state $s'$ (already computed)
• $\log P(s | s')$: the log transition probability from $s'$ to $s$
• $\log P(w_i | s)$: the log emission probability of word $w_i$ from state $s$

Log probabilities are negative since probabilities are between 0 and 1, and $\log(x) < 0$ when $x < 1$. The best path has the highest (least negative) log probability. The $\max$ operation is unchanged because log is monotonic: if $a > b$, then $\log(a) > \log(b)$. This means finding the maximum log probability gives the same answer as finding the maximum probability.

Both methods find the identical optimal path, confirming that log-space computation preserves correctness. The probability-space result (best_prob) and the exponentiated log-space result match to numerical precision. For the 100-token demonstration, probability space underflows to exactly zero while log space produces a meaningful value of -230.26. This is why all production Viterbi implementations work in log space.

The visualization makes the underflow problem concrete. With a typical per-step probability of 0.1, probability-space values hit the float64 minimum around sequence length 308 and become indistinguishable from zero. Log-space values, in contrast, decrease linearly and remain perfectly representable regardless of sequence length. Even a 1000-word document produces a log probability around -2300, large but entirely computable.

Log-space computation is essential for practical implementations. The algorithm produces identical paths but with stable numerics.

The Viterbi algorithm's efficiency comes from its dynamic programming structure. Let's analyze the time and space requirements.

For a sequence of length $n$ with $k$ possible states:

1. Initialization: $O(k)$ operations to set up the first column
2. Recursion: For each of the $n-1$ remaining positions, we compute $k$ cells. Each cell considers $k$ predecessors. Total: $O(n \cdot k^2)$
3. Backtracking: $O(n)$ to trace back through the path

The overall time complexity is $O(n \cdot k^2)$, which is linear in sequence length and quadratic in the number of states.

The benchmarks confirm the theoretical complexity. Doubling sequence length roughly doubles runtime (linear scaling), while doubling the number of states roughly quadruples runtime (quadratic scaling). Even with 20 states and 200-word sequences, Viterbi completes in tens of milliseconds, enabling real-time tagging in production systems.

The algorithm stores:

1. Viterbi table: $O(n \cdot k)$ scores
2. Backpointers: $O(n \cdot k)$ references

Total space is $O(n \cdot k)$. However, if we only need the final path and not the full table, we can reduce space to $O(k)$ by keeping only two columns at a time, though we'd need to run a second pass to recover the path.

A production implementation incorporates several optimizations beyond the basic algorithm. Let's build a more efficient version.

This implementation uses NumPy's vectorized operations to compute scores for all states simultaneously, avoiding Python loops over states. The key insight is that the inner loop over predecessors becomes a matrix operation.

The NumPy implementation produces the same results as the dictionary-based version but runs significantly faster for large state spaces. The log probability is negative (as expected for log probabilities) and represents the joint probability of the path and observations. The vectorized approach eliminates Python loops over states, leveraging NumPy's optimized C implementations.

Many HMMs have sparse transition matrices where most transitions have zero probability. We can exploit this by only considering valid transitions:

For BIO tagging where I-tags can only follow B-tags or I-tags of the same type, sparse transitions significantly reduce computation. With 9 entity types, you have 19 tags: one O tag plus B and I tags for each entity type (9 × 2 + 1 = 19). A dense transition matrix has $19^2 = 361$ entries, but valid transitions are far fewer since I-LOC cannot follow B-PER.

While Viterbi finds the single best path, many applications need the $k$ best paths. Neural sequence models often use beam search, which maintains the top $k$ partial hypotheses at each step. Beam search is essentially a pruned version of Viterbi.

The connection to Viterbi is direct: if we set the beam width to equal the number of states, beam search becomes equivalent to Viterbi. Smaller beam widths trade optimality for speed.

The comparison reveals the tradeoff between beam width and optimality. Greedy search (k=1) may find a suboptimal path because it prunes too aggressively. As beam width increases, we retain more hypotheses and are more likely to find the optimal solution. With beam width equal to the number of states, beam search becomes exact Viterbi. In practice, small beam widths (3-10) often find near-optimal solutions while running much faster than exact Viterbi.

Beam search with $k=1$ is greedy decoding. As $k$ increases, we approach Viterbi's optimal solution. The tradeoff is computation time: beam search is $O(n \cdot k \cdot |S|)$ where $k$ is the beam width, compared to Viterbi's $O(n \cdot |S|^2)$.

In sequence labeling, certain tag transitions are invalid. An I-LOC tag cannot follow a B-PER tag. Rather than learning these constraints from data, we can enforce them during decoding by setting invalid transition probabilities to zero (or negative infinity in log space).

The constraints follow BIO semantics: I-tags can only continue entities of the same type, while B-tags and O can follow any tag. Invalid transitions like O → I-PER (starting an entity with a continuation tag) or B-PER → I-LOC (continuing a PER entity as LOC) are explicitly forbidden. By setting these transition probabilities to negative infinity during Viterbi decoding, we guarantee the output is always a well-formed BIO sequence.

The heatmap reveals the structure of BIO constraints at a glance. The first row and column (O tag) show that O can follow anything and anything can precede O. The B-tags (columns 2, 4, 6) are fully green because they can follow any tag since they start new entities. The distinctive pattern appears in the I-tag columns: each I-tag has exactly two green cells, corresponding to its matching B and I tags. This sparsity (25 valid out of 49 possible transitions) enables significant optimization during Viterbi decoding.

Constrained decoding ensures the output is always a valid BIO sequence, even if the model's raw scores would suggest otherwise. This is particularly important for neural models that don't explicitly model transition constraints.

The Viterbi algorithm is optimal for first-order Markov models where the current state depends only on the previous state. Several limitations arise in practice:

• Higher-order dependencies exist in natural language. The optimal tag for a word might depend on tags two or three positions back. Extending Viterbi to second-order Markov models is possible but increases complexity to $O(n \cdot k^3)$ since we must track pairs of previous states.

• Non-local features are common in neural sequence models. A BiLSTM encodes the entire sentence before making predictions, so the score at each position depends on the full context, not just the previous tag. CRF layers on top of neural models use Viterbi for decoding, but the emission scores come from the neural network rather than simple probability tables.

• Approximate inference becomes necessary with large state spaces. In machine translation, the "states" are partial translations with exponentially many possibilities. Beam search provides a practical approximation, trading optimality for tractability.

• Label dependencies beyond transitions sometimes matter. In some tagging schemes, the entity type of a word might depend on entity types earlier in the sentence, not just the immediately preceding tag. Such long-range dependencies require more sophisticated inference methods.