7  Computation and numerical practice

The algorithms of Chapter 5 are exact in exact arithmetic. In floating point they are not, and a naive implementation of the forward recursion fails on sequences of moderate length. This chapter covers what has to be done to make the mathematics survive contact with a computer, together with the complexity facts that determine what is feasible.

Zucchini et al. (2016) give these matters unusual and welcome attention; most of what follows makes their advice explicit and adds the error analysis.

By the end of this chapter you will be able to

• explain why the forward recursion underflows, and fix it by scaling or by working in logarithms;
• state the time and memory cost of each algorithm, and exploit sparsity in \(\boldsymbol{\Gamma}\) when it is there;
• vectorise the recursions so the inner loop is linear algebra rather than scalar arithmetic;
• test an implementation using four invariants that hold regardless of the data.

7.1 Underflow and the two remedies

7.1.1 The problem

The unnormalised forward variable \(\tilde\alpha_t(i) = p(S_t = i, y_{1:t})\) is a joint density over \(t\) observations. Each step multiplies by a transition probability and an emission density, so its magnitude decays geometrically. If the typical per-step factor is \(10^{-2}\), then \(\tilde\alpha_t \sim 10^{-2t}\).

Double precision holds normalised numbers down to about \(10^{-308}\). Underflow therefore begins around \(t \approx 154\) and is complete shortly after: the vector becomes exactly zero, the likelihood is reported as zero, and the log-likelihood as \(-\infty\). Sequences of length a few hundred are entirely ordinary, so this is not an edge case.

Note that the problem is intrinsic to the quantity, not to the algorithm. The true value of \(\tilde\alpha_{500}(i)\) genuinely is around \(10^{-1000}\), and no rearrangement of the arithmetic will make it representable. Both remedies work by never forming it.

7.1.2 Scaling

Normalise the vector at each step, as in Proposition 5.1: propagate \(\hat{\boldsymbol{\alpha}}_t\) summing to one and accumulate \(\ell_T = \sum_t \log c_t\).

Cost per step is one vector–matrix product (\(K^2\) multiply–adds), \(K\) multiplications for the emission reweighting, \(K\) additions for the sum, \(K\) divisions, and one logarithm. The overhead beyond the unavoidable \(K^2\) work is \(\mathcal{O}(K)\), so scaling is close to free.

Accuracy is good. Each \(c_t\) is a predictive density of order one, so its logarithm is computed with full relative precision, and the total error in \(\ell_T\) accumulates as \(\mathcal{O}(T\varepsilon_{\text{mach}})\) in absolute terms. For \(T = 10^6\) that is around \(10^{-10}\) — negligible against any statistical uncertainty.

Scaling is the default. It is what Zucchini et al. (2016) implement, and it produces the one-step predictive densities needed for residual diagnostics as a by-product.

7.1.3 The log domain

Store logarithms throughout and combine with log-sum-exp: \[ a_{t+1}(j) = \log f_j(y_{t+1}) + \operatorname{logsumexp}_i \bigl\{ a_t(i) + \log\gamma_{ij} \bigr\}, \qquad \operatorname{logsumexp}(z) = m + \log\!\sum_i e^{z_i - m},\ \ m = \max_i z_i. \] Subtracting the maximum guarantees the largest exponential is \(e^0 = 1\); terms that underflow to zero were negligible anyway.

Cost per step is \(K^2\) additions, \(K\) maxima over \(K\) elements, \(K^2\) exponentials and \(K\) logarithms. Transcendental functions are perhaps an order of magnitude more expensive than multiplication, so this is several times slower than scaling.

What it buys is robustness in cases where scaling struggles: emission densities that vary over an enormous dynamic range within one time step, and structural zeros in \(\boldsymbol{\Gamma}\), which become \(-\infty\) and propagate correctly without special-casing. It is also the natural representation when the same code must run Viterbi, since that already works additively.

TipWhich to use

Scaling for routine work: faster, accurate, and it yields the predictive densities you want anyway. The log domain when the model has structural zeros, when emissions can be extremely small under some states, or when robustness matters more than speed. Implementing both and checking they agree to a few digits on a test problem is a good use of an afternoon.

7.1.4 Backward and Viterbi

The backward variables need the same treatment, and the cleanest option is the synchronised scaling of Chapter 5: divide \(\hat{\boldsymbol{\beta}}_t\) by the forward scale factor \(c_{t+1}\), so that \(\gamma_t(i) = \hat\alpha_t(i)\hat\beta_t(i)\) exactly with no further normalisation.

Viterbi needs no special handling beyond working in logarithms, since \(\max\) commutes with \(\log\) and no log-sum-exp arises. It is the numerically best-behaved of the three recursions.

7.2 Complexity

7.2.1 Single sequence

Cost of exact inference

Operation	Time	Memory
Forward, likelihood only	\(\mathcal{O}(K^2T)\)	\(\mathcal{O}(K)\)
Forward–backward, smoothing	\(\mathcal{O}(K^2T)\)	\(\mathcal{O}(KT)\)
Viterbi	\(\mathcal{O}(K^2T)\)	\(\mathcal{O}(KT)\)
FFBS, one trajectory draw	\(\mathcal{O}(K^2T)\)	\(\mathcal{O}(KT)\)

The \(K^2\) factor is the vector–matrix product and the \(T\) factor is the sweep; both are unavoidable for exact inference in a dense model. The memory asymmetry matters: computing the likelihood needs only the current vector, but anything requiring a backward pass must retain the forward quantities.

For very long sequences, checkpointing stores the forward vector every \(\sqrt{T}\) steps and recomputes intermediate segments during the backward pass, giving \(\mathcal{O}(K\sqrt{T})\) memory for roughly twice the time. Rarely needed, but the right answer when it is.

7.2.2 Estimation

One EM iteration is one forward–backward pass plus M-step updates that cost \(\mathcal{O}(K^2T)\) for the transition counts and \(\mathcal{O}(KT)\) for typical emission families. So an iteration is \(\mathcal{O}(K^2T)\), and a full fit is that times the number of iterations — which, as Chapter 6 noted, can be large when states overlap.

With \(N\) independent sequences of average length \(T\), cost per iteration is \(\mathcal{O}(NK^2T)\), and the sequences are independent so the E-step parallelises exactly. This is the easiest speedup available and should be the first one taken.

The quadratic growth in \(K\) is what keeps state numbers small. Going from \(K=5\) to \(K=20\) multiplies the cost sixteenfold, and the parameter count — \(K(K-1)\) transition parameters alone — grows just as fast, so the statistical problem worsens in step with the computational one.

7.2.3 Exploiting sparsity

If \(\boldsymbol{\Gamma}\) has many structural zeros, as in left-to-right models for speech or profile HMMs for sequence alignment, the recursion touches only the non-zero entries and the cost falls to \(\mathcal{O}(|\mathcal{E}|T)\) where \(|\mathcal{E}|\) is the number of permitted transitions. For a banded matrix allowing transitions only to nearby states this is \(\mathcal{O}(KT)\) — linear in \(K\). Sparse representations are worth the implementation effort whenever the topology is constrained.

7.3 Vectorisation and hardware

The forward recursion is inherently sequential in \(t\): step \(t+1\) needs the output of step \(t\). Within a step, however, everything is a dense linear algebra operation, and that is where the performance is.

• Vectorise over states. Write the step as \(\hat{\boldsymbol{\alpha}}_t^\top = \hat{\boldsymbol{\alpha}}_{t-1}^\top \boldsymbol{\Gamma} \mathbf{Q}(y_t)\) and let a BLAS routine do it, rather than writing loops over \(i\) and \(j\). In R or Python this alone can be a hundredfold improvement over the naive double loop.
• Batch over sequences. With \(N\) sequences, stack their forward vectors into an \(N \times K\) matrix and advance all of them with one matrix–matrix product. This converts a memory-bound operation into a compute-bound one and is where GPUs help.
• Precompute emission matrices. The values \(f_i(y_t)\) for all \((i,t)\) can be computed once as a \(T \times K\) matrix before the recursion begins. For expensive emission densities this is a large saving, since EM recomputes them at every iteration otherwise.

7.4 Approximate inference

Exact inference is cheap enough that approximation is needed only outside the finite-state, moderate-\(K\) setting.

Beam search for decoding. Retain only the \(B\) highest-scoring partial paths at each step, giving \(\mathcal{O}(BKT)\). The result is not guaranteed optimal, and the error is unbounded in principle, but with a generous beam the loss is usually negligible. Standard in large-vocabulary speech recognition, where \(K\) runs to millions.

Pruning. Drop states whose forward probability falls below a threshold relative to the current maximum. Cheaper than beam search and easier to implement, with the same caveat.

Variational approximation. Restrict \(q\) in the bound 6.1 to a tractable family and maximise. For a plain HMM this is unnecessary — the exact posterior over states is already tractable — but for the factorial and switching models of Chapter 12 it is the main option (Ghahramani and Hinton 2000; Wainwright and Jordan 2008).

Sequential Monte Carlo. For continuous state spaces, or discrete spaces too large for \(\mathcal{O}(K^2)\) steps.

Online EM. For streaming data, update sufficient statistics recursively with a decreasing gain \(\eta_t\) rather than sweeping the whole series each iteration. Convergence rests on the ergodicity assumed in Chapter 8. Cappé et al. (2005, ch. 11) develops the theory.

7.5 Implementation checks

Bugs in HMM code are easy to write and hard to see, because wrong answers still look like probabilities. Four checks catch most of them.

1. Likelihood consistency across \(t\). By Theorem 5.4, \(\sum_i \tilde\alpha_t(i)\beta_t(i)\) equals \(L_T\) for every \(t\). Compute it at several \(t\) and confirm agreement. This catches scaling errors immediately.
2. Marginal consistency of pairwise probabilities. Proposition 5.2 requires \(\sum_j \xi_t(i,j) = \gamma_t(i)\) and \(\sum_i \xi_t(i,j) = \gamma_{t+1}(j)\).
3. Monotonicity of EM. By Theorem 6.1 the log-likelihood cannot decrease. If it does, the M-step is wrong, or a constraint such as stationarity is being imposed incorrectly — see the warning in Chapter 6.
4. Recovery on simulated data. Simulate from known parameters with \(T\) large and confirm the estimates converge to them, up to a permutation of the state labels. This is the only check that tests the model as a whole, and it is worth automating.

A fifth, for samplers: run with no data and confirm the output reproduces the prior.

7.6 Summary

• Forward variables underflow within a few hundred time steps. Scaling (fast, accurate, gives predictive densities) and log-domain arithmetic (slower, more robust, handles structural zeros) both fix it exactly.
• Exact inference costs \(\mathcal{O}(K^2T)\) time; memory is \(\mathcal{O}(K)\) for the likelihood alone and \(\mathcal{O}(KT)\) once a backward pass is needed.
• Sparsity in \(\boldsymbol{\Gamma}\) reduces the cost to \(\mathcal{O}(|\mathcal{E}|T)\) and is worth exploiting when the topology is constrained.
• The recursion is sequential in \(t\) but vectorises perfectly over states and batches over sequences.
• Four cheap invariants — likelihood consistency across \(t\), marginal consistency of \(\xi_t\), EM monotonicity, and recovery on simulated data — catch most implementation errors.

7.7 Further reading

Zucchini et al. (2016, chs. 3–5 and appendices) give practical guidance and working R code with scaling handled exactly as described here. Cappé et al. (2005, chs. 5, 11) cover numerical considerations and online algorithms rigorously. Murphy (2012, ch. 17) discusses the numerical issues alongside the algorithms. For approximate methods, Wainwright and Jordan (2008) on variational inference and Chopin and Papaspiliopoulos (2020) on sequential Monte Carlo.

Cappé, Olivier, Éric Moulines, and Tobias Rydén. 2005. Inference in Hidden Markov Models. Springer Series in Statistics. Springer. https://doi.org/10.1007/0-387-28982-8.

Chopin, Nicolas, and Omiros Papaspiliopoulos. 2020. An Introduction to Sequential Monte Carlo. Springer. https://doi.org/10.1007/978-3-030-47845-2.

Ghahramani, Zoubin, and Geoffrey E. Hinton. 2000. “Variational Learning for Switching State-Space Models.” Neural Computation 12 (4): 831–64. https://doi.org/10.1162/089976600300015619.

Murphy, Kevin P. 2012. Machine Learning: A Probabilistic Perspective. MIT Press.

Wainwright, Martin J., and Michael I. Jordan. 2008. “Graphical Models, Exponential Families, and Variational Inference.” Foundations and Trends in Machine Learning 1 (1–2): 1–305. https://doi.org/10.1561/2200000001.

Zucchini, Walter, Iain L. MacDonald, and Roland Langrock. 2016. Hidden Markov Models for Time Series: An Introduction Using R. 2nd ed. Chapman; Hall/CRC. https://doi.org/10.1201/b20790.