Rabiner 1989 - Implementation Issues - Scaling, Multiple Sequences and Initialisation

Summary

Forward variables underflow exponentially with $T$, so scale ${\alpha }_t$ by $c_t=1/{\sum }_i{\alpha }_t(i)$; the scale factors cancel exactly in the reestimation formulas and $\log P(O\mid \lambda )=-{\sum }_t\log c_t$. Multiple observation sequences are handled by summing expected counts, each weighted by $1/P_k$. Good initial emission estimates matter, small-sample estimates need interpolation or floors, and model size is chosen by experience.

Overview

Section V is explicitly practical (“seat-of-the-pants experience gained from working with HMMs over the last several years”) but contains exact derivations for scaling.

Main Content

Scaling (eqs 91-103)

${\alpha }_t(i)$ is a sum of terms $\prod a\prod b$ each below 1, so for $t\gtrsim 100$ the dynamic range exceeds any machine precision. Define, with ${\tilde{\alpha }}_t(i)={\sum }_j{\hat{\alpha }}_{t-1}(j)a_{ji}{b}_i(O_t)$ and scale $c_t=1/{\sum }_i{\tilde{\alpha }}_t(i)$ (eq. 91): ${\hat{\alpha }}_t(i)=c_t{\tilde{\alpha }}_t(i)$. Then ${\hat{\alpha }}_t(i)=\left({\prod }_{s\le t}{c}_s\right){\alpha }_t(i)=C_t{\alpha }_t(i)$ and ${\hat{\alpha }}_t(i)={\alpha }_t(i)/{\sum }_j{\alpha }_t(j)$. Scaled backward variables use the same factors: ${\hat{\beta }}_t(i)=c_t{\beta }_t(i)$ (${\hat{\beta }}_t=D_t{\beta }_t$ with $D_t={\prod }_{s\ge t}{c}_s$). In the reestimation of $a_{ij}$ the factor $C_t{D}_{t+1}={\prod }_{s=1}^T{c}_s=C_T$ is independent of $t$ and cancels between numerator and denominator, giving the exact unscaled formula in terms of scaled variables (eq. 98). Since ${\prod }_t{c}_t{\sum }_i{\alpha }_T(i)=1$:

$$P(O\mid \lambda )=\frac{1}{{\prod }_{t=1}^T{c}_t},\log P(O\mid \lambda )=-\sum _{t=1}^T\log c_t.$$

Scaling can be applied only when needed (set $c_t=1$ at other times). Viterbi needs no scaling if done in logs (${\varphi }_t$ recursion, eqs 104-105).

Multiple observation sequences (eqs 106-111)

For left-right models a single sequence cannot train all parameters. With $K$ independent sequences ${\mathbf{O}}^{(k)}$, maximise $P(\mathbf{O}\mid \lambda )={\prod }_{k}P({\mathbf{O}}^{(k)}\mid \lambda )={\prod }_k{P}_k$. Add the per-sequence expected counts:

$${\bar{a}}_{ij}=\frac{{\sum }_k\frac{1}{P_k}{\sum }_{t=1}^{T_k-1}{\alpha }_t^k(i)a_{ij}{b}_j(O_{t+1}^{(k)}){\beta }_{t+1}^k(j)}{{\sum }_k\frac{1}{P_k}{\sum }_{t=1}^{T_k-1}{\alpha }_t^k(i){\beta }_t^k(i)},{\bar{b}}_j(\ell )=\frac{{\sum }_k\frac{1}{P_k}{\sum }_{t:O_t=v_{\ell }}{\alpha }_t^k(j){\beta }_t^k(j)}{{\sum }_k\frac{1}{P_k}{\sum }_t{\alpha }_t^k(j){\beta }_t^k(j)}.$$

${\pi }_i$ is not reestimated for left-right (${\pi }_1=1$). With scaling, the same scale factors cancel in $1/P_k$ (eq. 111). This is the batched analogue of E24’s “padded multi-track HMM” (elk tracks as independent sequences with shared parameters).

Initial estimates (Sec. V.C)

No simple answer for reaching the global maximum. Random (subject to constraints) or uniform initialisation of $\pi$ and $A$ suffices in almost all cases. For $B$, good initial estimates are helpful in the discrete-symbol case and essential in the continuous (mixture) case, obtained by manual segmentation, maximum-likelihood segmentation with averaging, or segmental $k$-means with clustering.

Insufficient training data (Sec. V.D; eqs 112-113; Fig. 10)

(1) Enlarge the training set (often impractical); (2) reduce model size (often impossible); (3) interpolate with a smaller model ${\lambda }^{\prime }$ for which data suffice: $\overline{\stackrel{ˉ}{\lambda }}=ϵ\lambda +(1-ϵ){\lambda }^{\prime }$ (eq. 112), with $ϵ$ estimated by deleted interpolation (Jelinek and Mercer 1980) by treating the interpolation as an expanded HMM with a neutral state $\tilde{s}$ (Fig. 10) and running forward-backward on held-out partitions; (4) floors: require $b_j(k)\ge \delta$ (eq. 113a) or $U_{jk}(r,r)\ge \delta$ (113b), with a post-processing step that corrects violations and rescales.

Choice of model (Sec. V.E). No theoretically correct rule: model type (ergodic or left-right), number of states, observation type (discrete or continuous, number of mixtures) are chosen by the physics of the signal and by experiment. A modern Bayesian answer is cross-validation or integrated LOO over the hidden path (E24).

Why scaling matters (numpy)

import numpy as np
rng = np.random.default_rng(0)
N, T = 3, 2000
A = np.full((N,N), 1/N); b = rng.random((T,N))*0.2             # tiny emission likelihoods
alpha = np.ones(N)/N*b[0]; logp = 0.0
for t in range(1, T):
    alpha = (alpha @ A)*b[t]
    c = alpha.sum(); alpha /= c; logp += np.log(c)             # scale c_t^{-1}; logp = -sum log c_t (eq. 103)
print(logp)   # finite; unscaled alpha would be ~ 0.1^2000 -> underflow to 0.0 in float64

Connections

Log-space forward algorithms (logsumexp) in PyMC/JAX are the modern replacement for explicit scaling: Rabiner 1989 - Problem 1 - Forward and Backward Algorithms. Gradient-based Bayesian fitting avoids Baum-Welch initialisation issues but inherits multimodality as label switching.

See Also

Rabiner 1989 - HMM Tutorial - Overview, Rabiner 1989 - Problem 3 - Baum-Welch Reestimation