Forward–backward algorithm

Forward probabilities

The following description will use matrices of probability values instead of probability distributions. However, the forward-backward algorithm can generally be applied to both continuous and discrete probability models.

We transform the probability distributions related to a given hidden Markov model into matrix notation as follows. The transition probabilities ${\displaystyle \mathbf{P}(X_t\mid X_{t-1})}$ of a given random variable ${\displaystyle X_t}$ representing all possible states in the hidden Markov model will be represented by the matrix ${\displaystyle \mathbf{T}}$ where the column index ${\displaystyle j}$ will represent the target state and the row index ${\displaystyle i}$ represents the start state. A transition from row-vector state ${\displaystyle {\pi }_{\mathbf{t}}}$ to the incremental row-vector state ${\displaystyle {\pi }_{\mathbf{t}\mathbf{+}\mathbf{1}}}$ is written as ${\displaystyle {\pi }_{\mathbf{t}\mathbf{+}\mathbf{1}}={\pi }_{\mathbf{t}}\mathbf{T}}$. The example below represents a system where the probability of staying in the same state after each step is 70% and the probability of transitioning to the other state is 30%. The transition matrix is then:

${\displaystyle \mathbf{T}=\left(\begin{array}{cc}0.7& 0.3\\ 0.3& 0.7\end{array}\right)}$

In a typical Markov model, we would multiply a state vector by this matrix to obtain the probabilities for the subsequent state. In a hidden Markov model the state is unknown, and we instead observe events associated with the possible states. An event matrix of the form:

${\displaystyle \mathbf{B}=\left(\begin{array}{cc}0.9& 0.1\\ 0.2& 0.8\end{array}\right)}$

provides the probabilities for observing events given a particular state. In the above example, event 1 will be observed 90% of the time if we are in state 1 while event 2 has a 10% probability of occurring in this state. In contrast, event 1 will only be observed 20% of the time if we are in state 2 and event 2 has an 80% chance of occurring. Given an arbitrary row-vector describing the state of the system (${\displaystyle \pi }$), the probability of observing event j is then:

${\displaystyle \mathbf{P}(O=j)=\sum _i{\pi }_i{B}_{i,j}}$

The probability of a given state leading to the observed event j can be represented in matrix form by multiplying the state row-vector (${\displaystyle \pi }$) with an observation matrix (${\displaystyle {\mathbf{O}}_{\mathbf{j}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(B_{\ast ,o_j})}$) containing only diagonal entries. Continuing the above example, the observation matrix for event 1 would be:

${\displaystyle {\mathbf{O}}_{\mathbf{1}}=\left(\begin{array}{cc}0.9& 0.0\\ 0.0& 0.2\end{array}\right)}$

This allows us to calculate the new unnormalized probabilities state vector ${\displaystyle {\pi }^{\prime }}$ through Bayes rule, weighting by the likelihood that each element of ${\displaystyle \pi }$ generated event 1 as:

${\displaystyle {\pi }^{\prime }=\pi {\mathbf{O}}_{\mathbf{1}}}$

We can now make this general procedure specific to our series of observations. Assuming an initial state vector ${\displaystyle {\pi }_0}$, (which can be optimized as a parameter through repetitions of the forward-backward procedure), we begin with ${\displaystyle {\mathbf{f}}_{\mathbf{0}\mathbf{:}\mathbf{0}}={\pi }_0}$, then updating the state distribution and weighting by the likelihood of the first observation:

${\displaystyle {\mathbf{f}}_{\mathbf{0}\mathbf{:}\mathbf{1}}={\pi }_0\mathbf{T}{\mathbf{O}}_{{\mathbf{o}}_{\mathbf{1}}}}$

This process can be carried forward with additional observations using:

${\displaystyle {\mathbf{f}}_{\mathbf{0}\mathbf{:}\mathbf{t}}={\mathbf{f}}_{\mathbf{0}\mathbf{:}\mathbf{t}\mathbf{-}\mathbf{1}}\mathbf{T}{\mathbf{O}}_{{\mathbf{o}}_{\mathbf{t}}}}$

This value is the forward unnormalized probability vector. The i'th entry of this vector provides:

${\displaystyle {\mathbf{f}}_{\mathbf{0}\mathbf{:}\mathbf{t}}(i)=\mathbf{P}(o_1,o_2,\dots ,o_t,X_t=x_i|{\pi }_0)}$

Typically, we will normalize the probability vector at each step so that its entries sum to 1. A scaling factor is thus introduced at each step such that:

${\displaystyle {\hat{\mathbf{f}}}_{\mathbf{0}\mathbf{:}\mathbf{t}}=c_t^{-1}{\hat{\mathbf{f}}}_{\mathbf{0}\mathbf{:}\mathbf{t}\mathbf{-}\mathbf{1}}\mathbf{T}{\mathbf{O}}_{{\mathbf{o}}_{\mathbf{t}}}}$

where ${\displaystyle {\hat{\mathbf{f}}}_{\mathbf{0}\mathbf{:}\mathbf{t}\mathbf{-}\mathbf{1}}}$ represents the scaled vector from the previous step and ${\displaystyle c_t}$ represents the scaling factor that causes the resulting vector's entries to sum to 1. The product of the scaling factors is the total probability for observing the given events irrespective of the final states:

${\displaystyle \mathbf{P}(o_1,o_2,\dots ,o_t|{\pi }_0)=\prod _{s=1}^t{c}_s}$

This allows us to interpret the scaled probability vector as:

${\displaystyle {\hat{\mathbf{f}}}_{\mathbf{0}\mathbf{:}\mathbf{t}}(i)=\frac{{\mathbf{f}}_{\mathbf{0}\mathbf{:}\mathbf{t}}(i)}{\prod _{s=1}^t{c}_s}=\frac{\mathbf{P}(o_1,o_2,\dots ,o_t,X_t=x_i|{\pi }_0)}{\mathbf{P}(o_1,o_2,\dots ,o_t|{\pi }_0)}=\mathbf{P}(X_t=x_i|o_1,o_2,\dots ,o_t,{\pi }_0)}$

We thus find that the product of the scaling factors provides us with the total probability for observing the given sequence up to time t and that the scaled probability vector provides us with the probability of being in each state at this time.

Backward probabilities

A similar procedure can be constructed to find backward probabilities. These intend to provide the probabilities:

${\displaystyle {\mathbf{b}}_{\mathbf{t}\mathbf{:}\mathbf{T}}(i)=\mathbf{P}(o_{t+1},o_{t+2},\dots ,o_T|X_t=x_i)}$

That is, we now want to assume that we start in a particular state (${\displaystyle X_t=x_i}$), and we are now interested in the probability of observing all future events from this state. Since the initial state is assumed as given (i.e. the prior probability of this state = 100%), we begin with:

${\displaystyle {\mathbf{b}}_{\mathbf{T}\mathbf{:}\mathbf{T}}=[111\dots {]}^T}$

Notice that we are now using a column vector while the forward probabilities used row vectors. We can then work backwards using:

${\displaystyle {\mathbf{b}}_{\mathbf{t}\mathbf{-}\mathbf{1}\mathbf{:}\mathbf{T}}=\mathbf{T}{\mathbf{O}}_{\mathbf{t}}{\mathbf{b}}_{\mathbf{t}\mathbf{:}\mathbf{T}}}$

While we could normalize this vector as well so that its entries sum to one, this is not usually done. Noting that each entry contains the probability of the future event sequence given a particular initial state, normalizing this vector would be equivalent to applying Bayes' theorem to find the likelihood of each initial state given the future events (assuming uniform priors for the final state vector). However, it is more common to scale this vector using the same ${\displaystyle c_t}$ constants used in the forward probability calculations. ${\displaystyle {\mathbf{b}}_{\mathbf{T}\mathbf{:}\mathbf{T}}}$ is not scaled, but subsequent operations use:

${\displaystyle {\hat{\mathbf{b}}}_{\mathbf{t}\mathbf{-}\mathbf{1}\mathbf{:}\mathbf{T}}=c_t^{-1}\mathbf{T}{\mathbf{O}}_{\mathbf{t}}{\hat{\mathbf{b}}}_{\mathbf{t}\mathbf{:}\mathbf{T}}}$

where ${\displaystyle {\hat{\mathbf{b}}}_{\mathbf{t}\mathbf{:}\mathbf{T}}}$ represents the previous, scaled vector. This result is that the scaled probability vector is related to the backward probabilities by:

${\displaystyle {\hat{\mathbf{b}}}_{\mathbf{t}\mathbf{:}\mathbf{T}}(i)=\frac{{\mathbf{b}}_{\mathbf{t}\mathbf{:}\mathbf{T}}(i)}{\prod _{s=t+1}^T{c}_s}}$

This is useful because it allows us to find the total probability of being in each state at a given time, t, by multiplying these values:

${\displaystyle {\gamma }_{\mathbf{t}}(i)=\mathbf{P}(X_t=x_i|o_1,o_2,\dots ,o_T,{\pi }_0)=\frac{\mathbf{P}(o_1,o_2,\dots ,o_T,X_t=x_i|{\pi }_0)}{\mathbf{P}(o_1,o_2,\dots ,o_T|{\pi }_0)}=\frac{{\mathbf{f}}_{\mathbf{0}\mathbf{:}\mathbf{t}}(i)\cdot {\mathbf{b}}_{\mathbf{t}\mathbf{:}\mathbf{T}}(i)}{\prod _{s=1}^T{c}_s}={\hat{\mathbf{f}}}_{\mathbf{0}\mathbf{:}\mathbf{t}}(i)\cdot {\hat{\mathbf{b}}}_{\mathbf{t}\mathbf{:}\mathbf{T}}(i)}$

To understand this, we note that ${\displaystyle {\mathbf{f}}_{\mathbf{0}\mathbf{:}\mathbf{t}}(i)\cdot {\mathbf{b}}_{\mathbf{t}\mathbf{:}\mathbf{T}}(i)}$ provides the probability for observing the given events in a way that passes through state ${\displaystyle x_i}$ at time t. This probability includes the forward probabilities covering all events up to time t as well as the backward probabilities which include all future events. This is the numerator we are looking for in our equation, and we divide by the total probability of the observation sequence to normalize this value and extract only the probability that ${\displaystyle X_t=x_i}$. These values are sometimes called the "smoothed values" as they combine the forward and backward probabilities to compute a final probability.

The values ${\displaystyle {\gamma }_{\mathbf{t}}(i)}$ thus provide the probability of being in each state at time t. As such, they are useful for determining the most probable state at any time. The term "most probable state" is somewhat ambiguous. While the most probable state is the most likely to be correct at a given point, the sequence of individually probable states is not likely to be the most probable sequence. This is because the probabilities for each point are calculated independently of each other. They do not take into account the transition probabilities between states, and it is thus possible to get states at two moments (t and t+1) that are both most probable at those time points but which have very little probability of occurring together, i.e. ${\displaystyle \mathbf{P}(X_t=x_i,X_{t+1}=x_j)\ne \mathbf{P}(X_t=x_i)\mathbf{P}(X_{t+1}=x_j)}$. The most probable sequence of states that produced an observation sequence can be found using the Viterbi algorithm.

Example

This example takes as its basis the umbrella world in Russell & Norvig 2010 Chapter 15 pp. 567 in which we would like to infer the weather given observation of another person either carrying or not carrying an umbrella. We assume two possible states for the weather: state 1 = rain, state 2 = no rain. We assume that the weather has a 70% chance of staying the same each day and a 30% chance of changing. The transition probabilities are then:

${\displaystyle \mathbf{T}=\left(\begin{array}{cc}0.7& 0.3\\ 0.3& 0.7\end{array}\right)}$

We also assume each state generates one of two possible events: event 1 = umbrella, event 2 = no umbrella. The conditional probabilities for these occurring in each state are given by the probability matrix:

${\displaystyle \mathbf{B}=\left(\begin{array}{cc}0.9& 0.1\\ 0.2& 0.8\end{array}\right)}$

We then observe the following sequence of events: {umbrella, umbrella, no umbrella, umbrella, umbrella} which we will represent in our calculations as:

${\displaystyle {\mathbf{O}}_{\mathbf{1}}=\left(\begin{array}{cc}0.9& 0.0\\ 0.0& 0.2\end{array}\right){\mathbf{O}}_{\mathbf{2}}=\left(\begin{array}{cc}0.9& 0.0\\ 0.0& 0.2\end{array}\right){\mathbf{O}}_{\mathbf{3}}=\left(\begin{array}{cc}0.1& 0.0\\ 0.0& 0.8\end{array}\right){\mathbf{O}}_{\mathbf{4}}=\left(\begin{array}{cc}0.9& 0.0\\ 0.0& 0.2\end{array}\right){\mathbf{O}}_{\mathbf{5}}=\left(\begin{array}{cc}0.9& 0.0\\ 0.0& 0.2\end{array}\right)}$

Note that ${\displaystyle {\mathbf{O}}_{\mathbf{3}}}$ differs from the others because of the "no umbrella" observation.

In computing the forward probabilities we begin with:

${\displaystyle {\mathbf{f}}_{\mathbf{0}\mathbf{:}\mathbf{0}}=\left(\begin{array}{cc}0.5& 0.5\end{array}\right)}$

which is our prior state vector indicating that we don't know which state the weather is in before our observations. While a state vector should be given as a row vector, we will use the transpose of the matrix so that the calculations below are easier to read. Our calculations are then written in the form:

${\displaystyle ({\hat{\mathbf{f}}}_{\mathbf{0}\mathbf{:}\mathbf{t}}{)}^T=c_t^{-1}{\mathbf{O}}_{\mathbf{t}}(\mathbf{T}{)}^T({\hat{\mathbf{f}}}_{\mathbf{0}\mathbf{:}\mathbf{t}\mathbf{-}\mathbf{1}}{)}^T}$

instead of:

${\displaystyle {\hat{\mathbf{f}}}_{\mathbf{0}\mathbf{:}\mathbf{t}}=c_t^{-1}{\hat{\mathbf{f}}}_{\mathbf{0}\mathbf{:}\mathbf{t}\mathbf{-}\mathbf{1}}\mathbf{T}{\mathbf{O}}_{\mathbf{t}}}$

Notice that the transformation matrix is also transposed, but in our example the transpose is equal to the original matrix. Performing these calculations and normalizing the results provides:

${\displaystyle ({\hat{\mathbf{f}}}_{\mathbf{0}\mathbf{:}\mathbf{1}}{)}^T=c_1^{-1}\left(\begin{array}{cc}0.9& 0.0\\ 0.0& 0.2\end{array}\right)\left(\begin{array}{cc}0.7& 0.3\\ 0.3& 0.7\end{array}\right)\left(\begin{array}{c}0.5000\\ 0.5000\end{array}\right)=c_1^{-1}\left(\begin{array}{c}0.4500\\ 0.1000\end{array}\right)=\left(\begin{array}{c}0.8182\\ 0.1818\end{array}\right)}$

${\displaystyle ({\hat{\mathbf{f}}}_{\mathbf{0}\mathbf{:}\mathbf{2}}{)}^T=c_2^{-1}\left(\begin{array}{cc}0.9& 0.0\\ 0.0& 0.2\end{array}\right)\left(\begin{array}{cc}0.7& 0.3\\ 0.3& 0.7\end{array}\right)\left(\begin{array}{c}0.8182\\ 0.1818\end{array}\right)=c_2^{-1}\left(\begin{array}{c}0.5645\\ 0.0745\end{array}\right)=\left(\begin{array}{c}0.8834\\ 0.1166\end{array}\right)}$

${\displaystyle ({\hat{\mathbf{f}}}_{\mathbf{0}\mathbf{:}\mathbf{3}}{)}^T=c_3^{-1}\left(\begin{array}{cc}0.1& 0.0\\ 0.0& 0.8\end{array}\right)\left(\begin{array}{cc}0.7& 0.3\\ 0.3& 0.7\end{array}\right)\left(\begin{array}{c}0.8834\\ 0.1166\end{array}\right)=c_3^{-1}\left(\begin{array}{c}0.0653\\ 0.2772\end{array}\right)=\left(\begin{array}{c}0.1907\\ 0.8093\end{array}\right)}$

${\displaystyle ({\hat{\mathbf{f}}}_{\mathbf{0}\mathbf{:}\mathbf{4}}{)}^T=c_4^{-1}\left(\begin{array}{cc}0.9& 0.0\\ 0.0& 0.2\end{array}\right)\left(\begin{array}{cc}0.7& 0.3\\ 0.3& 0.7\end{array}\right)\left(\begin{array}{c}0.1907\\ 0.8093\end{array}\right)=c_4^{-1}\left(\begin{array}{c}0.3386\\ 0.1247\end{array}\right)=\left(\begin{array}{c}0.7308\\ 0.2692\end{array}\right)}$

${\displaystyle ({\hat{\mathbf{f}}}_{\mathbf{0}\mathbf{:}\mathbf{5}}{)}^T=c_5^{-1}\left(\begin{array}{cc}0.9& 0.0\\ 0.0& 0.2\end{array}\right)\left(\begin{array}{cc}0.7& 0.3\\ 0.3& 0.7\end{array}\right)\left(\begin{array}{c}0.7308\\ 0.2692\end{array}\right)=c_5^{-1}\left(\begin{array}{c}0.5331\\ 0.0815\end{array}\right)=\left(\begin{array}{c}0.8673\\ 0.1327\end{array}\right)}$

For the backward probabilities, we start with:

${\displaystyle {\mathbf{b}}_{\mathbf{5}\mathbf{:}\mathbf{5}}=\left(\begin{array}{c}1.0\\ 1.0\end{array}\right)}$

We are then able to compute (using the observations in reverse order and normalizing with different constants):

${\displaystyle {\hat{\mathbf{b}}}_{\mathbf{4}\mathbf{:}\mathbf{5}}=\alpha \left(\begin{array}{cc}0.7& 0.3\\ 0.3& 0.7\end{array}\right)\left(\begin{array}{cc}0.9& 0.0\\ 0.0& 0.2\end{array}\right)\left(\begin{array}{c}1.0000\\ 1.0000\end{array}\right)=\alpha \left(\begin{array}{c}0.6900\\ 0.4100\end{array}\right)=\left(\begin{array}{c}0.6273\\ 0.3727\end{array}\right)}$

${\displaystyle {\hat{\mathbf{b}}}_{\mathbf{3}\mathbf{:}\mathbf{5}}=\alpha \left(\begin{array}{cc}0.7& 0.3\\ 0.3& 0.7\end{array}\right)\left(\begin{array}{cc}0.9& 0.0\\ 0.0& 0.2\end{array}\right)\left(\begin{array}{c}0.6273\\ 0.3727\end{array}\right)=\alpha \left(\begin{array}{c}0.4175\\ 0.2215\end{array}\right)=\left(\begin{array}{c}0.6533\\ 0.3467\end{array}\right)}$

${\displaystyle {\hat{\mathbf{b}}}_{\mathbf{2}\mathbf{:}\mathbf{5}}=\alpha \left(\begin{array}{cc}0.7& 0.3\\ 0.3& 0.7\end{array}\right)\left(\begin{array}{cc}0.1& 0.0\\ 0.0& 0.8\end{array}\right)\left(\begin{array}{c}0.6533\\ 0.3467\end{array}\right)=\alpha \left(\begin{array}{c}0.1289\\ 0.2138\end{array}\right)=\left(\begin{array}{c}0.3763\\ 0.6237\end{array}\right)}$

${\displaystyle {\hat{\mathbf{b}}}_{\mathbf{1}\mathbf{:}\mathbf{5}}=\alpha \left(\begin{array}{cc}0.7& 0.3\\ 0.3& 0.7\end{array}\right)\left(\begin{array}{cc}0.9& 0.0\\ 0.0& 0.2\end{array}\right)\left(\begin{array}{c}0.3763\\ 0.6237\end{array}\right)=\alpha \left(\begin{array}{c}0.2745\\ 0.1889\end{array}\right)=\left(\begin{array}{c}0.5923\\ 0.4077\end{array}\right)}$

${\displaystyle {\hat{\mathbf{b}}}_{\mathbf{0}\mathbf{:}\mathbf{5}}=\alpha \left(\begin{array}{cc}0.7& 0.3\\ 0.3& 0.7\end{array}\right)\left(\begin{array}{cc}0.9& 0.0\\ 0.0& 0.2\end{array}\right)\left(\begin{array}{c}0.5923\\ 0.4077\end{array}\right)=\alpha \left(\begin{array}{c}0.3976\\ 0.2170\end{array}\right)=\left(\begin{array}{c}0.6469\\ 0.3531\end{array}\right)}$

Finally, we will compute the smoothed probability values. These results must also be scaled so that its entries sum to 1 because we did not scale the backward probabilities with the ${\displaystyle c_t}$'s found earlier. The backward probability vectors above thus actually represent the likelihood of each state at time t given the future observations. Because these vectors are proportional to the actual backward probabilities, the result has to be scaled an additional time.

${\displaystyle ({\gamma }_{\mathbf{0}}{)}^T=\alpha \left(\begin{array}{c}0.5000\\ 0.5000\end{array}\right)\circ \left(\begin{array}{c}0.6469\\ 0.3531\end{array}\right)=\alpha \left(\begin{array}{c}0.3235\\ 0.1765\end{array}\right)=\left(\begin{array}{c}0.6469\\ 0.3531\end{array}\right)}$

${\displaystyle ({\gamma }_{\mathbf{1}}{)}^T=\alpha \left(\begin{array}{c}0.8182\\ 0.1818\end{array}\right)\circ \left(\begin{array}{c}0.5923\\ 0.4077\end{array}\right)=\alpha \left(\begin{array}{c}0.4846\\ 0.0741\end{array}\right)=\left(\begin{array}{c}0.8673\\ 0.1327\end{array}\right)}$

${\displaystyle ({\gamma }_{\mathbf{2}}{)}^T=\alpha \left(\begin{array}{c}0.8834\\ 0.1166\end{array}\right)\circ \left(\begin{array}{c}0.3763\\ 0.6237\end{array}\right)=\alpha \left(\begin{array}{c}0.3324\\ 0.0728\end{array}\right)=\left(\begin{array}{c}0.8204\\ 0.1796\end{array}\right)}$

${\displaystyle ({\gamma }_{\mathbf{3}}{)}^T=\alpha \left(\begin{array}{c}0.1907\\ 0.8093\end{array}\right)\circ \left(\begin{array}{c}0.6533\\ 0.3467\end{array}\right)=\alpha \left(\begin{array}{c}0.1246\\ 0.2806\end{array}\right)=\left(\begin{array}{c}0.3075\\ 0.6925\end{array}\right)}$

${\displaystyle ({\gamma }_{\mathbf{4}}{)}^T=\alpha \left(\begin{array}{c}0.7308\\ 0.2692\end{array}\right)\circ \left(\begin{array}{c}0.6273\\ 0.3727\end{array}\right)=\alpha \left(\begin{array}{c}0.4584\\ 0.1003\end{array}\right)=\left(\begin{array}{c}0.8204\\ 0.1796\end{array}\right)}$

${\displaystyle ({\gamma }_{\mathbf{5}}{)}^T=\alpha \left(\begin{array}{c}0.8673\\ 0.1327\end{array}\right)\circ \left(\begin{array}{c}1.0000\\ 1.0000\end{array}\right)=\alpha \left(\begin{array}{c}0.8673\\ 0.1327\end{array}\right)=\left(\begin{array}{c}0.8673\\ 0.1327\end{array}\right)}$

Notice that the value of ${\displaystyle {\gamma }_{\mathbf{0}}}$ is equal to ${\displaystyle {\hat{\mathbf{b}}}_{\mathbf{0}\mathbf{:}\mathbf{5}}}$ and that ${\displaystyle {\gamma }_{\mathbf{5}}}$ is equal to ${\displaystyle {\hat{\mathbf{f}}}_{\mathbf{0}\mathbf{:}\mathbf{5}}}$. This follows naturally because both ${\displaystyle {\hat{\mathbf{f}}}_{\mathbf{0}\mathbf{:}\mathbf{5}}}$ and ${\displaystyle {\hat{\mathbf{b}}}_{\mathbf{0}\mathbf{:}\mathbf{5}}}$ begin with uniform priors over the initial and final state vectors (respectively) and take into account all of the observations. However, ${\displaystyle {\gamma }_{\mathbf{0}}}$ will only be equal to ${\displaystyle {\hat{\mathbf{b}}}_{\mathbf{0}\mathbf{:}\mathbf{5}}}$ when our initial state vector represents a uniform prior (i.e. all entries are equal). When this is not the case ${\displaystyle {\hat{\mathbf{b}}}_{\mathbf{0}\mathbf{:}\mathbf{5}}}$ needs to be combined with the initial state vector to find the most likely initial state. We thus find that the forward probabilities by themselves are sufficient to calculate the most likely final state. Similarly, the backward probabilities can be combined with the initial state vector to provide the most probable initial state given the observations. The forward and backward probabilities need only be combined to infer the most probable states between the initial and final points.

The calculations above reveal that the most probable weather state on every day except for the third one was "rain". They tell us more than this, however, as they now provide a way to quantify the probabilities of each state at different times. Perhaps most importantly, our value at ${\displaystyle {\gamma }_{\mathbf{5}}}$ quantifies our knowledge of the state vector at the end of the observation sequence. We can then use this to predict the probability of the various weather states tomorrow as well as the probability of observing an umbrella.

Performance

The forward–backward algorithm runs with time complexity ${\displaystyle O(S^{2}T)}$ in space ${\displaystyle O(ST)}$, where ${\displaystyle T}$ is the length of the time sequence and ${\displaystyle S}$ is the number of symbols in the state alphabet. The algorithm can also run in constant space with time complexity ${\displaystyle O(S^2{T}^2)}$ by recomputing values at each step. For comparison, a brute-force procedure would generate all possible ${\displaystyle S^T}$ state sequences and calculate the joint probability of each state sequence with the observed series of events, which would have time complexity ${\displaystyle O(T\cdot S^T)}$. Brute force is intractable for realistic problems, as the number of possible hidden node sequences typically is extremely high.

An enhancement to the general forward-backward algorithm, called the Island algorithm, trades smaller memory usage for longer running time, taking ${\displaystyle O(S^{2}T\log T)}$ time and ${\displaystyle O(S\log T)}$ memory. Furthermore, it is possible to invert the process model to obtain an ${\displaystyle O(S)}$ space, ${\displaystyle O(S^{2}T)}$ time algorithm, although the inverted process may not exist or be ill-conditioned.

In addition, algorithms have been developed to compute ${\displaystyle {\mathbf{f}}_{\mathbf{0}\mathbf{:}\mathbf{t}\mathbf{+}\mathbf{1}}}$ efficiently through online smoothing such as the fixed-lag smoothing (FLS) algorithm.

Pseudocode

<b id="mwASE">algorithm</b> forward_backward <b id="mwASI">is</b> <b id="mwASM">input:</b> guessState int <i id="mwASQ">sequenceIndex</i> <b id="mwASU">output:</b> <i id="mwASY">result</i> <b id="mwASc">if</b> <i id="mwASg">sequenceIndex</i> is past the end of the sequence <b id="mwASk">then</b> <b id="mwASo">return</b> 1 <b id="mwASs">if</b> (<i id="mwASw">guessState</i>, <i id="mwAS0">sequenceIndex</i>) has been seen before <b id="mwAS4">then</b> <b id="mwAS8">return</b> saved result <i id="mwATA">result</i> := 0 <b id="mwATE">for each</b> neighboring state n: <i id="mwATI">result</i> := result + (transition probability from <i id="mwATM">guessState</i> to n given observation element at <i id="mwATQ">sequenceIndex</i>) × Backward(n, <i id="mwATU">sequenceIndex</i> + 1) save result for (<i id="mwATY">guessState</i>, <i id="mwATc">sequenceIndex</i>) <b id="mwATg">return</b> <i id="mwATk">result</i>

Python example

Given HMM (just like in Viterbi algorithm) represented in the Python programming language:

<span id="mwAUA"></span><span class="n" id="mwAUE">states</span> <span class="o" id="mwAUI">=</span> <span class="p" id="mwAUM">(</span><span class="s2" id="mwAUQ">"Healthy"</span><span class="p" id="mwAUU">,</span> <span class="s2" id="mwAUY">"Fever"</span><span class="p" id="mwAUc">)</span> <span class="n" id="mwAUg">end_state</span> <span class="o" id="mwAUk">=</span> <span class="s2" id="mwAUo">"E"</span> <span class="n" id="mwAUs">observations</span> <span class="o" id="mwAUw">=</span> <span class="p" id="mwAU0">(</span><span class="s2" id="mwAU4">"normal"</span><span class="p" id="mwAU8">,</span> <span class="s2" id="mwAVA">"cold"</span><span class="p" id="mwAVE">,</span> <span class="s2" id="mwAVI">"dizzy"</span><span class="p" id="mwAVM">)</span> <span class="n" id="mwAVQ">start_probability</span> <span class="o" id="mwAVU">=</span> <span class="p" id="mwAVY">{</span><span class="s2" id="mwAVc">"Healthy"</span><span class="p" id="mwAVg">:</span> <span class="mf" id="mwAVk">0.6</span><span class="p" id="mwAVo">,</span> <span class="s2" id="mwAVs">"Fever"</span><span class="p" id="mwAVw">:</span> <span class="mf" id="mwAV0">0.4</span><span class="p" id="mwAV4">}</span> <span class="n" id="mwAV8">transition_probability</span> <span class="o" id="mwAWA">=</span> <span class="p" id="mwAWE">{</span> <span class="s2" id="mwAWI">"Healthy"</span><span class="p" id="mwAWM">:</span> <span class="p" id="mwAWQ">{</span><span class="s2" id="mwAWU">"Healthy"</span><span class="p" id="mwAWY">:</span> <span class="mf" id="mwAWc">0.69</span><span class="p" id="mwAWg">,</span> <span class="s2" id="mwAWk">"Fever"</span><span class="p" id="mwAWo">:</span> <span class="mf" id="mwAWs">0.3</span><span class="p" id="mwAWw">,</span> <span class="s2" id="mwAW0">"E"</span><span class="p" id="mwAW4">:</span> <span class="mf" id="mwAW8">0.01</span><span class="p" id="mwAXA">},</span> <span class="s2" id="mwAXE">"Fever"</span><span class="p" id="mwAXI">:</span> <span class="p" id="mwAXM">{</span><span class="s2" id="mwAXQ">"Healthy"</span><span class="p" id="mwAXU">:</span> <span class="mf" id="mwAXY">0.4</span><span class="p" id="mwAXc">,</span> <span class="s2" id="mwAXg">"Fever"</span><span class="p" id="mwAXk">:</span> <span class="mf" id="mwAXo">0.59</span><span class="p" id="mwAXs">,</span> <span class="s2" id="mwAXw">"E"</span><span class="p" id="mwAX0">:</span> <span class="mf" id="mwAX4">0.01</span><span class="p" id="mwAX8">},</span> <span class="p" id="mwAYA">}</span> <span class="n" id="mwAYE">emission_probability</span> <span class="o" id="mwAYI">=</span> <span class="p" id="mwAYM">{</span> <span class="s2" id="mwAYQ">"Healthy"</span><span class="p" id="mwAYU">:</span> <span class="p" id="mwAYY">{</span><span class="s2" id="mwAYc">"normal"</span><span class="p" id="mwAYg">:</span> <span class="mf" id="mwAYk">0.5</span><span class="p" id="mwAYo">,</span> <span class="s2" id="mwAYs">"cold"</span><span class="p" id="mwAYw">:</span> <span class="mf" id="mwAY0">0.4</span><span class="p" id="mwAY4">,</span> <span class="s2" id="mwAY8">"dizzy"</span><span class="p" id="mwAZA">:</span> <span class="mf" id="mwAZE">0.1</span><span class="p" id="mwAZI">},</span> <span class="s2" id="mwAZM">"Fever"</span><span class="p" id="mwAZQ">:</span> <span class="p" id="mwAZU">{</span><span class="s2" id="mwAZY">"normal"</span><span class="p" id="mwAZc">:</span> <span class="mf" id="mwAZg">0.1</span><span class="p" id="mwAZk">,</span> <span class="s2" id="mwAZo">"cold"</span><span class="p" id="mwAZs">:</span> <span class="mf" id="mwAZw">0.3</span><span class="p" id="mwAZ0">,</span> <span class="s2" id="mwAZ4">"dizzy"</span><span class="p" id="mwAZ8">:</span> <span class="mf" id="mwAaA">0.6</span><span class="p" id="mwAaE">},</span> <span class="p" id="mwAaI">}</span>

We can write the implementation of the forward-backward algorithm like this:

<span id="mwAaY"></span><span class="k" id="mwAac">def</span><span class="w" id="mwAag"> </span><span class="nf" id="mwAak">fwd_bkw</span><span class="p" id="mwAao">(</span><span class="n" id="mwAas">observations</span><span class="p" id="mwAaw">,</span> <span class="n" id="mwAa0">states</span><span class="p" id="mwAa4">,</span> <span class="n" id="mwAa8">start_prob</span><span class="p" id="mwAbA">,</span> <span class="n" id="mwAbE">trans_prob</span><span class="p" id="mwAbI">,</span> <span class="n" id="mwAbM">emm_prob</span><span class="p" id="mwAbQ">,</span> <span class="n" id="mwAbU">end_st</span><span class="p" id="mwAbY">):</span> <span class="w" id="mwAbc"> </span><span class="sd" id="mwAbg">"""Forward–backward algorithm."""</span> <span class="c1" id="mwAbk"># Forward part of the algorithm</span> <span class="n" id="mwAbo">fwd</span> <span class="o" id="mwAbs">=</span> <span class="p" id="mwAbw">[]</span> <span class="k" id="mwAb0">for</span> <span class="n" id="mwAb4">i</span><span class="p" id="mwAb8">,</span> <span class="n" id="mwAcA">observation_i</span> <span class="ow" id="mwAcE">in</span> <span class="nb" id="mwAcI">enumerate</span><span class="p" id="mwAcM">(</span><span class="n" id="mwAcQ">observations</span><span class="p" id="mwAcU">):</span> <span class="n" id="mwAcY">f_curr</span> <span class="o" id="mwAcc">=</span> <span class="p" id="mwAcg">{}</span> <span class="k" id="mwAck">for</span> <span class="n" id="mwAco">st</span> <span class="ow" id="mwAcs">in</span> <span class="n" id="mwAcw">states</span><span class="p" id="mwAc0">:</span> <span class="k" id="mwAc4">if</span> <span class="n" id="mwAc8">i</span> <span class="o" id="mwAdA">==</span> <span class="mi" id="mwAdE">0</span><span class="p" id="mwAdI">:</span> <span class="c1" id="mwAdM"># base case for the forward part</span> <span class="n" id="mwAdQ">prev_f_sum</span> <span class="o" id="mwAdU">=</span> <span class="n" id="mwAdY">start_prob</span><span class="p" id="mwAdc">[</span><span class="n" id="mwAdg">st</span><span class="p" id="mwAdk">]</span> <span class="k" id="mwAdo">else</span><span class="p" id="mwAds">:</span> <span class="n" id="mwAdw">prev_f_sum</span> <span class="o" id="mwAd0">=</span> <span class="nb" id="mwAd4">sum</span><span class="p" id="mwAd8">(</span><span class="n" id="mwAeA">f_prev</span><span class="p" id="mwAeE">[</span><span class="n" id="mwAeI">k</span><span class="p" id="mwAeM">]</span> <span class="o" id="mwAeQ">*</span> <span class="n" id="mwAeU">trans_prob</span><span class="p" id="mwAeY">[</span><span class="n" id="mwAec">k</span><span class="p" id="mwAeg">][</span><span class="n" id="mwAek">st</span><span class="p" id="mwAeo">]</span> <span class="k" id="mwAes">for</span> <span class="n" id="mwAew">k</span> <span class="ow" id="mwAe0">in</span> <span class="n" id="mwAe4">states</span><span class="p" id="mwAe8">)</span> <span class="n" id="mwAfA">f_curr</span><span class="p" id="mwAfE">[</span><span class="n" id="mwAfI">st</span><span class="p" id="mwAfM">]</span> <span class="o" id="mwAfQ">=</span> <span class="n" id="mwAfU">emm_prob</span><span class="p" id="mwAfY">[</span><span class="n" id="mwAfc">st</span><span class="p" id="mwAfg">][</span><span class="n" id="mwAfk">observation_i</span><span class="p" id="mwAfo">]</span> <span class="o" id="mwAfs">*</span> <span class="n" id="mwAfw">prev_f_sum</span> <span class="n" id="mwAf0">fwd</span><span class="o" id="mwAf4">.</span><span class="n" id="mwAf8">append</span><span class="p" id="mwAgA">(</span><span class="n" id="mwAgE">f_curr</span><span class="p" id="mwAgI">)</span> <span class="n" id="mwAgM">f_prev</span> <span class="o" id="mwAgQ">=</span> <span class="n" id="mwAgU">f_curr</span> <span class="n" id="mwAgY">p_fwd</span> <span class="o" id="mwAgc">=</span> <span class="nb" id="mwAgg">sum</span><span class="p" id="mwAgk">(</span><span class="n" id="mwAgo">f_curr</span><span class="p" id="mwAgs">[</span><span class="n" id="mwAgw">k</span><span class="p" id="mwAg0">]</span> <span class="o" id="mwAg4">*</span> <span class="n" id="mwAg8">trans_prob</span><span class="p" id="mwAhA">[</span><span class="n" id="mwAhE">k</span><span class="p" id="mwAhI">][</span><span class="n" id="mwAhM">end_st</span><span class="p" id="mwAhQ">]</span> <span class="k" id="mwAhU">for</span> <span class="n" id="mwAhY">k</span> <span class="ow" id="mwAhc">in</span> <span class="n" id="mwAhg">states</span><span class="p" id="mwAhk">)</span> <span class="c1" id="mwAho"># Backward part of the algorithm</span> <span class="n" id="mwAhs">bkw</span> <span class="o" id="mwAhw">=</span> <span class="p" id="mwAh0">[]</span> <span class="k" id="mwAh4">for</span> <span class="n" id="mwAh8">i</span><span class="p" id="mwAiA">,</span> <span class="n" id="mwAiE">observation_i_plus</span> <span class="ow" id="mwAiI">in</span> <span class="nb" id="mwAiM">enumerate</span><span class="p" id="mwAiQ">(</span><span class="nb" id="mwAiU">reversed</span><span class="p" id="mwAiY">(</span><span class="n" id="mwAic">observations</span><span class="p" id="mwAig">[</span><span class="mi" id="mwAik">1</span><span class="p" id="mwAio">:]</span> <span class="o" id="mwAis">+</span> <span class="p" id="mwAiw">(</span><span class="kc" id="mwAi0">None</span><span class="p" id="mwAi4">,))):</span> <span class="n" id="mwAi8">b_curr</span> <span class="o" id="mwAjA">=</span> <span class="p" id="mwAjE">{}</span> <span class="k" id="mwAjI">for</span> <span class="n" id="mwAjM">st</span> <span class="ow" id="mwAjQ">in</span> <span class="n" id="mwAjU">states</span><span class="p" id="mwAjY">:</span> <span class="k" id="mwAjc">if</span> <span class="n" id="mwAjg">i</span> <span class="o" id="mwAjk">==</span> <span class="mi" id="mwAjo">0</span><span class="p" id="mwAjs">:</span> <span class="c1" id="mwAjw"># base case for backward part</span> <span class="n" id="mwAj0">b_curr</span><span class="p" id="mwAj4">[</span><span class="n" id="mwAj8">st</span><span class="p" id="mwAkA">]</span> <span class="o" id="mwAkE">=</span> <span class="n" id="mwAkI">trans_prob</span><span class="p" id="mwAkM">[</span><span class="n" id="mwAkQ">st</span><span class="p" id="mwAkU">][</span><span class="n" id="mwAkY">end_st</span><span class="p" id="mwAkc">]</span> <span class="k" id="mwAkg">else</span><span class="p" id="mwAkk">:</span> <span class="n" id="mwAko">b_curr</span><span class="p" id="mwAks">[</span><span class="n" id="mwAkw">st</span><span class="p" id="mwAk0">]</span> <span class="o" id="mwAk4">=</span> <span class="nb" id="mwAk8">sum</span><span class="p" id="mwAlA">(</span><span class="n" id="mwAlE">trans_prob</span><span class="p" id="mwAlI">[</span><span class="n" id="mwAlM">st</span><span class="p" id="mwAlQ">][</span><span class="n" id="mwAlU">l</span><span class="p" id="mwAlY">]</span> <span class="o" id="mwAlc">*</span> <span class="n" id="mwAlg">emm_prob</span><span class="p" id="mwAlk">[</span><span class="n" id="mwAlo">l</span><span class="p" id="mwAls">][</span><span class="n" id="mwAlw">observation_i_plus</span><span class="p" id="mwAl0">]</span> <span class="o" id="mwAl4">*</span> <span class="n" id="mwAl8">b_prev</span><span class="p" id="mwAmA">[</span><span class="n" id="mwAmE">l</span><span class="p" id="mwAmI">]</span> <span class="k" id="mwAmM">for</span> <span class="n" id="mwAmQ">l</span> <span class="ow" id="mwAmU">in</span> <span class="n" id="mwAmY">states</span><span class="p" id="mwAmc">)</span> <span class="n" id="mwAmg">bkw</span><span class="o" id="mwAmk">.</span><span class="n" id="mwAmo">insert</span><span class="p" id="mwAms">(</span><span class="mi" id="mwAmw">0</span><span class="p" id="mwAm0">,</span><span class="n" id="mwAm4">b_curr</span><span class="p" id="mwAm8">)</span> <span class="n" id="mwAnA">b_prev</span> <span class="o" id="mwAnE">=</span> <span class="n" id="mwAnI">b_curr</span> <span class="n" id="mwAnM">p_bkw</span> <span class="o" id="mwAnQ">=</span> <span class="nb" id="mwAnU">sum</span><span class="p" id="mwAnY">(</span><span class="n" id="mwAnc">start_prob</span><span class="p" id="mwAng">[</span><span class="n" id="mwAnk">l</span><span class="p" id="mwAno">]</span> <span class="o" id="mwAns">*</span> <span class="n" id="mwAnw">emm_prob</span><span class="p" id="mwAn0">[</span><span class="n" id="mwAn4">l</span><span class="p" id="mwAn8">][</span><span class="n" id="mwAoA">observations</span><span class="p" id="mwAoE">[</span><span class="mi" id="mwAoI">0</span><span class="p" id="mwAoM">]]</span> <span class="o" id="mwAoQ">*</span> <span class="n" id="mwAoU">b_curr</span><span class="p" id="mwAoY">[</span><span class="n" id="mwAoc">l</span><span class="p" id="mwAog">]</span> <span class="k" id="mwAok">for</span> <span class="n" id="mwAoo">l</span> <span class="ow" id="mwAos">in</span> <span class="n" id="mwAow">states</span><span class="p" id="mwAo0">)</span> <span class="c1" id="mwAo4"># Merging the two parts</span> <span class="n" id="mwAo8">posterior</span> <span class="o" id="mwApA">=</span> <span class="p" id="mwApE">[]</span> <span class="k" id="mwApI">for</span> <span class="n" id="mwApM">i</span> <span class="ow" id="mwApQ">in</span> <span class="nb" id="mwApU">range</span><span class="p" id="mwApY">(</span><span class="nb" id="mwApc">len</span><span class="p" id="mwApg">(</span><span class="n" id="mwApk">observations</span><span class="p" id="mwApo">)):</span> <span class="n" id="mwAps">posterior</span><span class="o" id="mwApw">.</span><span class="n" id="mwAp0">append</span><span class="p" id="mwAp4">({</span><span class="n" id="mwAp8">st</span><span class="p" id="mwAqA">:</span> <span class="n" id="mwAqE">fwd</span><span class="p" id="mwAqI">[</span><span class="n" id="mwAqM">i</span><span class="p" id="mwAqQ">][</span><span class="n" id="mwAqU">st</span><span class="p" id="mwAqY">]</span> <span class="o" id="mwAqc">*</span> <span class="n" id="mwAqg">bkw</span><span class="p" id="mwAqk">[</span><span class="n" id="mwAqo">i</span><span class="p" id="mwAqs">][</span><span class="n" id="mwAqw">st</span><span class="p" id="mwAq0">]</span> <span class="o" id="mwAq4">/</span> <span class="n" id="mwAq8">p_fwd</span> <span class="k" id="mwArA">for</span> <span class="n" id="mwArE">st</span> <span class="ow" id="mwArI">in</span> <span class="n" id="mwArM">states</span><span class="p" id="mwArQ">})</span> <span class="k" id="mwArU">assert</span> <span class="n" id="mwArY">p_fwd</span> <span class="o" id="mwArc">==</span> <span class="n" id="mwArg">p_bkw</span> <span class="k" id="mwArk">return</span> <span class="n" id="mwAro">fwd</span><span class="p" id="mwArs">,</span> <span class="n" id="mwArw">bkw</span><span class="p" id="mwAr0">,</span> <span class="n" id="mwAr4">posterior</span>

The function fwd_bkw takes the following arguments: x is the sequence of observations, e.g. ['normal', 'cold', 'dizzy']; states is the set of hidden states; a_0 is the start probability; a are the transition probabilities; and e are the emission probabilities.

For simplicity of code, we assume that the observation sequence x is non-empty and that a[i][j] and e[i][j] is defined for all states i,j.

In the running example, the forward-backward algorithm is used as follows:

<span id="mwAs4"></span><span class="k" id="mwAs8">def</span><span class="w" id="mwAtA"> </span><span class="nf" id="mwAtE">example</span><span class="p" id="mwAtI">():</span> <span class="k" id="mwAtM">return</span> <span class="n" id="mwAtQ">fwd_bkw</span><span class="p" id="mwAtU">(</span> <span class="n" id="mwAtY">observations</span><span class="p" id="mwAtc">,</span> <span class="n" id="mwAtg">states</span><span class="p" id="mwAtk">,</span> <span class="n" id="mwAto">start_probability</span><span class="p" id="mwAts">,</span> <span class="n" id="mwAtw">transition_probability</span><span class="p" id="mwAt0">,</span> <span class="n" id="mwAt4">emission_probability</span><span class="p" id="mwAt8">,</span> <span class="n" id="mwAuA">end_state</span><span class="p" id="mwAuE">,</span> <span class="p" id="mwAuI">)</span>

<span id="mwAuU"></span><span class="gp" id="mwAuY">>>> </span><span class="k" id="mwAuc">for</span> <span class="n" id="mwAug">line</span> <span class="ow" id="mwAuk">in</span> <span class="n" id="mwAuo">example</span><span class="p" id="mwAus">():</span> <span class="gp" id="mwAuw">... </span> <span class="nb" id="mwAu0">print</span><span class="p" id="mwAu4">(</span><span class="o" id="mwAu8">*</span><span class="n" id="mwAvA">line</span><span class="p" id="mwAvE">)</span> <span class="gp" id="mwAvI">... </span> <span class="go" id="mwAvM">{'Healthy': 0.3, 'Fever': 0.04000000000000001} {'Healthy': 0.0892, 'Fever': 0.03408} {'Healthy': 0.007518, 'Fever': 0.028120319999999997}</span> <span class="go" id="mwAvQ">{'Healthy': 0.0010418399999999998, 'Fever': 0.00109578} {'Healthy': 0.00249, 'Fever': 0.00394} {'Healthy': 0.01, 'Fever': 0.01}</span> <span class="go" id="mwAvU">{'Healthy': 0.8770110375573259, 'Fever': 0.1229889624426741} {'Healthy': 0.623228030950954, 'Fever': 0.3767719690490461} {'Healthy': 0.2109527048413057, 'Fever': 0.7890472951586943}</span>