Forward–backward algorithm

1 Problem setup and model structure

1.1 Hidden Markov models and latent state sequences

A hidden Markov model (HMM) models a system whose true state is not directly observed. Instead, the model assumes there is an underlying sequence of latent states that evolves over time, while observations are generated by emitting signals that depend on the current latent state. The forward–backward algorithm uses this layered structure to infer posterior distributions over the hidden states given the observed data.

1.2 Observation sequences and emission probabilities

Let an observation sequence be denoted by \(O = (O_1,\dots,O_T)\). Each observation \(O_t\) is assumed to be drawn from a distribution conditioned on the latent state at time \(t\). The emission probabilities (or emission density, in continuous cases) capture how likely each observation is under each latent state, enabling likelihood computations that depend on the entire observation sequence.

1.3 Markov assumptions and conditional independence

The defining simplification is the Markov property for the latent states: the next latent state depends only on the current latent state. Coupled with conditional independence of observations, this implies that once the latent state at time \(t\) is fixed, the observation \(O_t\) does not depend on other time steps. These assumptions factorize the joint probability of latent and observed sequences into a chain of local terms, which is what makes dynamic programming efficient.

1.4 Notation, conventions, and time indexing

Time is typically indexed from \(t=1\) to \(t=T\). The latent state space is usually finite with states indexed by \(i \in \{1,\dots,N\}\). Transition probabilities are written as \(a_{ij}\) for moving from state \(i\) at time \(t\) to state \(j\) at time \(t+1\). Initial state probabilities are written as \(\pi_i\). Emission likelihoods are denoted \(b_i(O_t)\), meaning the probability (or density) of observing \(O_t\) when the latent state is \(i\).

2 Forward pass (filtering)

2.1 Forward variables and their meaning

The forward pass computes quantities that summarize evidence up to each time step. These are usually called forward variables, often denoted \(\alpha_t(i)\), representing the probability of being in state \(i\) at time \(t\) together with producing the observed prefix up to \(t\).

2.1.1 Initialization at time t=1

At \(t=1\), the forward variable is computed from the initial distribution and the emission likelihood for the first observation: \[ \alpha_1(i) = \pi_i\, b_i(O_1). \] This initializes the recursion by combining prior belief about the starting state with how consistent that state is with the first observation.

2.2 Recursive computation over time

For each subsequent time \(t=2,\dots,T\), the algorithm aggregates contributions from all previous states: \[ \alpha_t(j) = b_j(O_t)\sum_{i=1}^N \alpha_{t-1}(i)\, a_{ij}. \] Conceptually, the sum accumulates all ways the latent chain could have arrived at state \(j\), weighted by the probability of the observation prefix.

2.3 Interpreting forward probabilities as partial likelihoods

Under the standard interpretation, \(\alpha_t(i)\) is an unnormalized belief about state \(i\) after seeing observations up to time \(t\). If one divides by the total probability of the observed prefix, the result becomes the filtering posterior \(P(X_t=i\mid O_{1:t})\). In practice, the unnormalized form is convenient for recursion and later combination with backward information.

2.4 Computational considerations and complexity

The recursion for each time step requires summing over \(i\) for each destination state \(j\), leading to \(O(N^2)\) operations per time step and \(O(TN^2)\) time overall. Memory can be reduced by storing only the previous \(\alpha_{t-1}(\cdot)\) and current \(\alpha_t(\cdot)\) values, giving \(O(N)\) space for the forward pass itself.

3 Backward pass (smoothing components)

3.1 Backward variables and their meaning

The backward pass provides complementary evidence from the future. Backward variables, often denoted \(\beta_t(i)\), represent the probability of observing the suffix from time \(t+1\) to \(T\), given that the latent state at time \(t\) is \(i\).

3.1.1 Initialization at final time T

A common convention is: \[ \beta_T(i)=1 \quad \text{for all } i, \] because once the sequence ends, there is no future observation to account for. This sets the boundary condition for backward recursion.

3.2 Recursive computation backward in time

For \(t=T-1,\dots,1\), the recursion is: \[ \beta_t(i)=\sum_{j=1}^N a_{ij}\, b_j(O_{t+1})\, \beta_{t+1}(j). \] Here, one sums over possible next states \(j\), weighting by the transition from \(i\) to \(j\), the likelihood of the next observation under \(j\), and the probability of the remaining future suffix.

3.3 Relation to partial likelihoods and future evidence

With this definition, \(\beta_t(i)\) can be seen as how strongly state \(i\) at time \(t\) is supported by the observations that follow. When later combined with the forward variables, this future support sharpens the posterior distribution over hidden states, producing smoothed estimates that use all observations.

3.4 Numerical stability issues

Probabilities in \(\alpha\) and \(\beta\) can underflow because they often involve products of many terms less than one. Standard solutions include scaling the variables at each time step or performing calculations in log space. Scaling is especially common in implementations because it preserves the recursion form while maintaining representable magnitudes.

4 Combining forward and backward information

4.1 Posterior state probabilities (smoothing)

Smoothing refers to estimating hidden states using both past and future evidence. The fundamental quantity is the posterior probability that the latent state at time \(t\) equals \(i\) given the full observation sequence.

4.1.1 Computing γ_t(i) from α and β

A typical expression is: \[ \gamma_t(i) = P(X_t=i\mid O) = \frac{\alpha_t(i)\beta_t(i)}{\sum_{k=1}^N \alpha_t(k)\beta_t(k)}. \] The denominator normalizes across states at time \(t\). The numerator combines the likelihood of the observed prefix compatible with \(i\) and the likelihood of the observed suffix compatible with \(i\).

4.2 Posterior transition probabilities

Beyond marginal state posteriors, many learning procedures require expectations over pairs of adjacent hidden states. This is captured by posterior transition probabilities.

4.2.1 Computing ξ_t(i,j) for parameter updates

Define: \[ \xi_t(i,j)=P(X_t=i, X_{t+1}=j\mid O). \] A standard computation is: \[ \xi_t(i,j)=\frac{\alpha_t(i)\, a_{ij}\, b_j(O_{t+1})\, \beta_{t+1}(j)}{\sum_{r=1}^N\sum_{s=1}^N \alpha_t(r)\, a_{rs}\, b_s(O_{t+1})\, \beta_{t+1}(s)}. \] These pairwise posteriors identify which state transitions are most supported by the observed data.

4.3 Consistency checks and marginalization identities

The combined quantities satisfy useful identities. For instance, marginalizing the pairwise posterior recovers the state posterior: \[ \gamma_t(i)=\sum_{j=1}^N \xi_t(i,j), \] and also \[ \gamma_{t+1}(j)=\sum_{i=1}^N \xi_t(i,j). \] These relationships help validate implementations and ensure that the normalization across time is handled correctly.

5 Likelihood and normalization

5.1 Computing the sequence likelihood P(O | model)

The forward variables can be used to obtain the likelihood of the entire observation sequence under the model. One common formula is: \[ P(O\mid \text{model}) = \sum_{i=1}^N \alpha_T(i). \] This works because \(\alpha_T(i)\) already encodes the probability of reaching state \(i\) at the end while generating the full observation sequence.

5.2 Normalized variants to prevent underflow

To avoid numerical underflow, one may scale \(\alpha_t(i)\) by a factor at each time step. If \(c_t\) denotes the scale, then one often computes: \[ \hat{\alpha}_t(i)=\frac{\alpha_t(i)}{c_t}. \] The scaling factors can be accumulated to recover the overall likelihood. Similarly, backward variables can be scaled using the same or compatible factors so that \(\gamma_t(i)\) and \(\xi_t(i,j)\) remain correctly normalized.

5.3 Scaling factors and log-space computations

Another approach is to compute in log space. For example, storing \(\log \alpha_t(i)\) helps manage underflow and allows use of the log-sum-exp trick to evaluate sums safely. While log-space introduces additional computation overhead, it offers robustness when sequence lengths are large or emission/transition probabilities are extremely small.

5.4 Verification using forward-only and forward–backward identities

A useful consistency check is that the likelihood computed from the forward pass matches the quantity implied by the combined forward–backward terms. With correct scaling and normalization, expressions like \(\sum_i \alpha_t(i)\beta_t(i)\) should relate consistently to \(P(O\mid \text{model})\) across different \(t\). These checks are practical safeguards when debugging.

6 Parameter estimation with Baum–Welch (EM)

6.1 Overview of the EM framework

The Baum–Welch algorithm is an expectation–maximization (EM) procedure specialized for HMMs. It iteratively improves parameter estimates by alternating between computing expected hidden-state usage under the current parameters (E-step) and updating parameters to maximize the expected complete-data log-likelihood (M-step). The forward–backward algorithm supplies the needed expectations.

6.2 E-step: expected sufficient statistics

Given current parameters, the E-step computes \(\gamma_t(i)\) and \(\xi_t(i,j)\). These quantities act as expected counts:

• \(\gamma_t(i)\) measures the expected occupancy of state \(i\) at time \(t\).
• \(\xi_t(i,j)\) measures the expected number of transitions from \(i\) to \(j\) around time \(t\).

For discrete emissions, one can further accumulate expected emission counts by summing \(\gamma_t(i)\) over times where each observation value occurs.

6.3 M-step: re-estimation of transition and emission parameters

In the M-step, parameters are updated using these expectations, typically with normalization constraints. For transitions: \[ \hat{a}_{ij} = \frac{\sum_{t=1}^{T-1} \xi_t(i,j)}{\sum_{t=1}^{T-1} \gamma_t(i)}. \] For initial probabilities: \[ \hat{\pi}_i = \gamma_1(i). \] For discrete emissions with parameters \(b_i(k)=P(O_t=k\mid X_t=i)\), one commonly updates: \[ \hat{b}_i(k)=\frac{\sum_{t: O_t=k} \gamma_t(i)}{\sum_{t=1}^{T} \gamma_t(i)}. \] Continuous-emission variants replace these with distribution-specific weighted updates.

6.4 Convergence criteria and stopping rules

EM is typically stopped when improvement in the log-likelihood falls below a threshold, when parameter changes become negligible, or after a fixed number of iterations. Because the likelihood is guaranteed to not decrease under exact EM updates, monitoring the log-likelihood provides a direct diagnostic for progress.

6.5 Common pitfalls (local maxima, initialization)

Baum–Welch can converge to local optima. Parameter initialization therefore has practical importance; multiple restarts are often used to improve chances of finding a better solution. Additionally, when states receive near-zero posterior probability early in training, updates can become unstable or overly confident, sometimes motivating smoothing priors or minimum-probability floors.

7 Extensions and related formulations

7.1 Discrete HMMs vs. continuous emissions

For discrete emissions, \(b_i(O_t)\) is a probability mass function, and learning updates can be expressed in terms of counts. With continuous observations, \(b_i(\cdot)\) is typically a density (e.g., Gaussian or mixture models), and the M-step becomes a weighted parameter estimation problem using \(\gamma_t(i)\) as weights. The forward–backward structure remains the same, while emission updates change.

7.2 Extensions to other chain-structured graphical models

The method generalizes to any model whose latent structure forms a directed or undirected chain with local dependencies compatible with dynamic programming. In such cases, “states” may represent more complex latent variables, but the computational principle remains: factorize the joint distribution and reuse intermediate results.

7.3 Connection to belief propagation on factor graphs

An HMM can be represented as a factor graph. Forward–backward corresponds to belief propagation (sum-product) on that graph, where messages propagate along the chain in both directions. This perspective clarifies why the algorithm is efficient: it performs exact marginalization on a graph with treewidth 1 (a chain).

7.4 Relation to Kalman filtering/smoothing (conceptual)

Although Kalman filtering is designed for linear dynamical systems with Gaussian noise, it shares the conceptual role of combining past and future information to infer latent variables. Where Kalman uses continuous Gaussian beliefs and linear algebra, forward–backward uses discrete (or non-Gaussian) probability computations. The analogy helps develop intuition, even though the mathematical machinery differs.

8 Practical implementation details

8.1 Algorithm pseudocode and data structures

A typical implementation has two phases: compute forward variables, then compute backward variables, then combine them to produce \(\gamma_t(i)\) and (if needed) \(\xi_t(i,j)\). Arrays of size \(N\) store \(\alpha_t\) and \(\beta_t\) at each time; if memory allows, one may store full time histories, otherwise use rolling buffers for \(\alpha\). Transition matrices of size \(N\times N\) and emission parameters are held in structured arrays.

8.2 Complexity analysis (time and memory)

Time complexity is dominated by the \(O(TN^2)\) nested summations in the recursions, both for forward and backward passes. For memory, a rolling implementation can achieve \(O(N)\) storage for forward and backward vectors, plus optional \(O(TN)\) storage for \(\gamma\) and \(\xi\) when required for parameter updates. When emissions are discrete, additional bookkeeping may be needed for expected emission counts.

8.3 Handling missing observations

If some \(O_t\) are missing, one can adjust the emission likelihood terms. A common strategy sets the emission likelihood to a constant that does not favor any state at that time, effectively integrating out the missing observation. Formally, this corresponds to using \(b_i(O_t)=P(\text{missing}\mid X_t=i)\) with a likelihood that is identical across \(i\), or by summing emission probabilities over plausible observed values if a missingness model is specified.

8.4 Batch processing for multiple sequences

When training on many sequences, the algorithm is run for each sequence and the resulting expected statistics are aggregated. Because sequences may have different lengths, implementations typically loop over sequences, compute per-sequence \(\gamma\) and \(\xi\), and then sum expectations for parameter updates. Vectorized implementations can improve throughput, especially when many sequences share the same state count.

8.5 Testing, debugging, and unit checks

Robust verification includes:

• Checking normalization: \(\sum_i \gamma_t(i)=1\) for each \(t\).
• Checking marginalization: \(\gamma_t(i)=\sum_j \xi_t(i,j)\).
• Recomputing likelihood in consistent ways using scaling factors.
• Ensuring boundary conditions: \(\beta_T(i)=1\) under the unscaled convention.

These tests catch indexing errors, normalization mistakes, and dimension mismatches.

9 Worked examples and intuition

9.1 Small HMM example with explicit calculations

Consider a two-state HMM (\(N=2\)) with \(T=3\) and a short observation sequence. Computing \(\alpha_1\) uses \(\pi_i\) and the emission likelihood for \(O_1\). Then \(\alpha_2\) and \(\alpha_3\) follow the recursive sum over predecessor states. The backward pass initializes \(\beta_3(i)=1\), then produces \(\beta_2\) and \(\beta_1\) via sums over successor states. Finally, \(\gamma_t(i)\) is obtained by multiplying \(\alpha_t(i)\) and \(\beta_t(i)\) and normalizing. Even with small numbers, the example illustrates how uncertainty propagates through time in both directions.

9.2 Visualizing forward/backward evidence flow

Intuitively, the forward pass pushes information from left to right: at time \(t\), \(\alpha_t(i)\) represents compatibility with having arrived at state \(i\) after matching the seen observations. The backward pass pulls information from right to left: \(\beta_t(i)\) represents compatibility with the remaining future observations if state \(i\) holds now. The smoothing posterior \(\gamma_t(i)\) reflects the intersection of these two evidence streams.

9.3 Interpreting posteriors and uncertainty

Posterior state probabilities express graded belief rather than a single discrete guess. When emissions are informative, \(\gamma_t(i)\) becomes peaked, indicating strong evidence for particular states. When emissions overlap significantly across states, the posterior remains diffuse, reflecting ambiguity that cannot be resolved from the data at hand. The forward–backward algorithm therefore naturally quantifies uncertainty.

9.4 Comparing naive enumeration vs. dynamic programming

Naively, one could enumerate all \(N^T\) possible hidden-state paths and compute their probabilities, then normalize to obtain posteriors. Forward–backward avoids this explosion by exploiting the chain structure: it reuses intermediate results \(\alpha_t(i)\) and \(\beta_t(i)\) so that each local transition contribution is counted many times implicitly. The result is a dramatic reduction from exponential to polynomial time.

10 Common variants and optimizations

10.1 Viterbi vs. forward–backward (MAP vs. marginals)

Forward–backward computes marginal posteriors of states and transitions (smoothing). The Viterbi algorithm instead finds the single most probable latent path, a maximum a posteriori (MAP) sequence. While both rely on dynamic programming, Viterbi uses maximization instead of summation, producing a different output: a best path rather than full probability distributions.

10.2 Beam/approximate methods in large state spaces

When the number of states is very large, exact forward–backward can be computationally expensive due to \(O(N^2)\) transitions. Approximate methods such as beam search restrict attention to a subset of states at each time step, trading exactness for speed. These approaches can be useful when emissions are strong enough that most states have negligible posterior probability.

10.3 Parallelization opportunities

The recursions over destination states involve sums over predecessor states, which can be parallelized across \(j\) (forward) or across \(i\) (backward). On modern hardware, vectorized linear algebra or GPU kernels can accelerate matrix-vector operations, especially when transition matrices are dense and state counts are moderate.

10.4 Memory optimization strategies (rolling arrays)

To lower memory use, implementations often store only the current and previous forward vectors and reuse space for backward computations. When training with Baum–Welch, one may need \(\gamma_t(i)\) and \(\xi_t(i,j)\) across time; this can be addressed by accumulating expected sufficient statistics on the fly without storing all intermediate arrays. Rolling arrays reduce peak memory while preserving correctness.