Viterbi Algorithm - Most likely sequence vs sequence of most likely states

I'm trying to understand why the following pseudo-code function is correct:

# probability == p. Tm: the transition matrix. Em: the emission matrix.
function viterbi(O, S, Π, Tm, Em): best_path
trellis ← matrix(length(S), length(O))  # To hold p. of each state given each observation.
# Determine each hidden state's p. at time 0…
for s in range(length(S)):
trellis[s, 0] ← Π[s] * Em[s, O[0]]
# …and afterwards, assuming each state's most likely prior state, k.
for o in range(1, length(O)):
for s in range(length(S)):
k ← argmax(k in trellis[k, o-1] * Tm[k, s] * Em[s, o])
trellis[s, o] ← trellis[k, o-1] * Tm[k, s] * Em[s, o]
best_path ← list()
for o in range(-1, -(length(O)+1), -1):  # Backtrack from last observation.
k ← argmax(k in trellis[k, o])  # Most likely state at o
best_path.insert(0, S[k])       # is noted for return.
return best_path


To me it seems that the best_path is just the sequence of most likely states and not the most likely sequence. This is taken from this Wikipedia page. There is another pseudo-code function on that Wikipedia page that stores back pointers to determine the best path, which the above function does not do.

Why are the two functions equivalent?

You’re right to question the equivalence—they’re not equal. And it’s not the sequence of most likely states, either. It’s actually worse than that, Tom.

The Viterbi algorithm is split into a forward pass and a backward pass. The forward pass fills in the trellis. The backward pass reconstructs the most likely sequence of states efficiently using the memoization from the forward pass. Each cell $$(s, o)$$ contains the maximum (over all prefixes state sequences) of the joint probability of the observations and states, given that the current state is $$s$$.

But there isn’t enough information here to determine which sequence that is, yet. One way to do this is to track backpointers of which predecessor led to our best score. The other, less common approach is to use automatic differentiation to determine it. In the forward pass, we would use computations that track their gradients (à la Torch). In the backward pass, we would select each predecessor state with highest gradient. This is the one that won in our $$\max$$ from the forward pass.

The algorithm you show in the article doesn’t do either of these things. At each point, it picks the state that has the best overall history. It ignores, though, whether this state is most compatible with its successor!

So far, we’ve established that this isn’t the Viterbi algorithm—the code doesn’t construct the most likely sequence of states. But it’s not the sequence of most likely states either.

Remember that the most likely state at a given point is the $$\arg\max_{S_i} \sum_{S_{1, \ldots, i-1, i+1, \ldots N}} P(S_{1, \ldots, N}, O_{1, \ldots, N})$$. That is—which state is most compatible with both what’s before it and what follows? We marginalize out (i.e., sum over) the other state candidates. Recovering this sequence of states is called posterior decoding. It gives us the states with highest posterior marginal probability. The resulting sequence may not be likely itself, though! You can efficiently compute the sequence of most likely states with the forward-backward algorithm.

But this differs from what Wikipedia shows, too! Instead, it’s finding at each point the state that’s in the most likely sequence up to time $$o$$. It's concerned only with the prefix, not what follows. (What follows is what Viterbi helps with.) In simpler words, it's "What did I think was the best when I'd gotten this far during the forward pass?" I can’t think of a time that people would ever want this.

If you take a look at the Talk section of the Wikipedia article, you can see that there’s a lot of contention about the pseudocode and the Python code, different camps saying it works or doesn’t. Clearly your expertise would be valued to improve the article.