The Mathematics of Hidden Markov Models

Mathematics of Hidden Markov Models: filtering, smoothing, and inference via Expectation-Maximization.

Model Definition

A Hidden Markov Model (HMM) consists of a sequence of hidden states $z_{1:T}$ with $z_t \in \{1,\dots,K\}$ evolving over time, each generating an observation $x_t$.

The factorization encodes two key conditional independence assumptions: (1) the evolution of latent states is 1-Markov, and (2) each observation depends only on the contemporaneous latent state.

Inference

Since only $x_{1:T}$ is observed, there are two major inference tasks: (1) filtering $p(z_t=k \mid x_{1:t})$ and (2) smoothing $p(z_{t}=k \mid x_{1:T})$. With smoothing, observations beyond state $z_t$ are used to compute the posterior. From these results, other tasks like future latent states or observations become straightforward. We will first look at filtering.

Filtering: Inferring Current Latent State

Filtering is the task of computing the posterior $p(z_t=k \mid x_{1:t})$ over latent states given observations up to time $t$. Writing

it suffices to track the numerator $p(x_{1:t}, z_{t}=k)$. Using the conditional independence implied by (1), $x_{t} \perp (z_{t-1}, x_{1:t-1}) \mid z_{t}$ and $z_{t} \perp x_{1:t-1} \mid z_{t-1}$, we have

Definition (Forward variable). Define $\alpha_t(k) \equiv p(x_{1:t}, z_t=k)$ for $t \geq 1$.

Proposition (Forward recursion). The forward variables satisfy, for $t \geq 2$:

for $t = 1$:

The marginal likelihood (probability of the data) $p(x_{1:t})$ can be recovered by summing over the latent states $\alpha$:

Which gives us

Smoothing: Inferring an Intermediate Latent State

In the context of HMMs, smoothing refers to finding $p(z_{t} \mid x_{1:T})$, making use of the full observation sequence to infer an intermediate latent state. Intuitively, even observations beyond $z_t$ can still carry information about the state of $z_t$ through the Markov chain. Starting with the object of interest

Focusing on the numerator:

which is using the conditional independence due to the HMM graph, $x_{1:t} \perp x_{t+1:T} \mid z_{t}$:

Definition (Backward variable). Define $\beta_t(i) \equiv p(x_{t+1:T} \mid z_t = i)$.

Proposition (Backward recursion). For $1 \leq t \leq T-1$:

with terminal condition $\beta_T(i) = 1$.

Proof.

and now because $x_{t+2:T} \perp (x_{t+1}, z_{t}) \mid z_{t+1}$ and $x_{t+1} \perp z_t \mid z_{t+1}$:

and for $t=T$, for convenience we define $\beta_{T}(j) = 1$ for $j \in \{1, \ldots, K\}$. $\square$

Definition (Smoother). Define $\gamma_{t}(i) \equiv p(z_{t}=i \mid x_{1:T})$.

Proposition (Smoothing distribution). Putting together (5), for any $t$,

Parameter Fitting

The model parameters consist of the initial state distribution, the transition distribution between latent states, and the emission distribution associated with each latent state.

That is, the model parameters are:

Following the method of maximum likelihood estimation, we'd like to maximize the probability of the data given the parameters.

The sum over latent sequences appears inside the logarithm, which couples the parameters across all possible $\mathbf{z}_n$. This prevents the objective from decomposing into separable terms, and hence no closed-form maximizer exists.

Expectation-Maximization - E-step

Under EM we want to maximize the expected complete-data log-likelihood under the current posterior:

Moving expectations inside:

Definition (Pairwise Smoother). Define $\xi_{t}(j,k) \equiv p(z_{t-1}=j,\, z_t=k \mid x_{1:T})$.

Writing $\xi_{t}(j,k)$ in terms of objects we already have:

To summarize:

We have now computed the expectation of the log-likelihood under the posterior of the latent variables.

Expectation Maximization - M-step

Now, we need to update the parameters of our model $\theta = (\pi, A, B)$. Since $Q$ is the sum of three terms of distinct parameters, we can maximize each term independently. Starting with maximizing the first term $Q_{\pi}$:

Initial distribution.

Setting up the Lagrangian:

Now using the constraint:

and because $\gamma_{1}(k) \equiv p(z_{1}=k \mid x_{1:T})$ the sum across all latent states $\sum_{k} p(z_{1}=k \mid x_{1:T}) = 1$ and hence

Now we take a brief detour to prove a classic result which can be reused.

Lemma. The solution to $\max_p \sum_k w_k \log p_k$ subject to $\sum_k p_k = 1$, $p_k \geq 0$, is

Proof. From the Lagrangian $\mathcal{L}(p, \lambda) = \sum_{k} w_{k} \log p_{k} + \lambda (\sum_{k} p_{k} - 1)$, take the derivative w.r.t. $p_{k}$: $\frac{\partial \mathcal{L}}{\partial p_{k}} = \frac{w_{k}}{p_{k}} + \lambda = 0$. Solving: $\frac{w_{k}}{p_{k}} = -\lambda \implies p_{k} = \frac{w_{k}}{-\lambda} \implies -\lambda = \sum_{k} w_{k}$. The constraint $\sum_k p_k = 1$ then forces $p_{k} = \frac{w_{k}}{\sum_{j} w_{j}}$. $\square$

Transition matrix. The result applies to $Q_{A}$:

Considering a particular $j$ we'd like to maximize

Making use of the Lemma, we see this is maximized by

and because

Therefore, substituting this result into the denominator:

Remark. The updated transition parameters become the normalized estimated/"soft" (from the E-step) counts of transitions from $j \to k$.

The final term to deal with is the emission parameters. The optimization depends on whether we choose emission to be categorical, Gaussian, or some other distribution. For a categorical distribution, we can follow a similar approach using the general result as we did with $Q_{A}$ and $Q_{\pi}$.

Let us also show the derivation for Gaussian emissions.

Emission parameters (Gaussian).

Let's start by finding $\mu_{k}$:

Now, solving for the variance term:

The updates (7) and (8) are the standard Gaussian MLE estimates, except each observation is weighted by the posterior probability that it was generated by state $k$.

EM Algorithm Summary (Vectorized)

Given current parameters $\theta^{\text{old}} = (\pi, A, B)$, define the emission likelihood vector

Let $\odot$ denote elementwise multiplication.

E-step. Compute the forward vectors $\bm{\alpha}_t$ and backward vectors $\bm{\beta}_t$ using

and

Then compute the smoothing probabilities

For $t=2,\dots,T$, compute the pairwise smoothing probabilities in matrix form

where $(\bm{\Xi}_t)_{jk} = \xi_t(j,k)$.

M-step. Update the parameters:

Then define the expected transition count matrix $\mathbf{N} = \sum_{t=2}^{T} \bm{\Xi}_t$, in which each row $j$ encodes the number of expected transitions out of $j$. Then the transition is the row update

Then for the Gaussian parameters, defining $\mathbf{x} \in \mathbb{R}^{T}$ as the vector of observations and $\Gamma \in \mathbb{R}^{T \times K}$ as $\bm{\gamma}_t$ stacked row-wise:

Define a residual matrix $\mathbf{R} \in \mathbb{R}^{T \times K}$ by

Then the variance update is

where $\mathbf{R}^{\odot 2}$ denotes elementwise squaring and division is elementwise.