The Learning Library
Contents

Module 02 — Implementing the Inference Engine From Scratch

Part I · Foundations & the Inference Engine · Status: Draft v0.1 Scope: From-scratch forward, scaling, log-space, and Viterbi · Prerequisites: 00, 01

Overview

Raw forward probabilities shrink with every time step until your floats give up. You fix this with scaling constants or log-space arithmetic before touching any library. This module implements both fixes from scratch in numpy.

You reuse THE SHARED TOY PROBLEM from Module 01: Rainy and Sunny states, Walk/Shop/Clean observations, and the same start, transition, and emission matrices. Every checkable number here must agree with Module 01 exactly.

By the end you verify your scaled forward pass against brute-force enumeration on a tiny example. Two independent implementations printing the same number is your first real inference-engine test.

How It Works

The forward pass multiplies one probability per time step, so the running product decays geometrically. Scaling divides each day’s scores by their sum and banks that sum as a constant. The true likelihood is the product of the banked constants.

How to read this: raw scores flow down the left branch toward zero, while the right branch rescales each day and accumulates the log-likelihood safely.

flowchart TD
    RAW[Raw forward scores] --> SHRINK[Scores shrink each step]
    SHRINK --> UF[Long runs hit zero]
    SCALE[Scaled forward scores] --> BANK[Bank daily sums]
    BANK --> LOGLIKE[Sum logs for likelihood]
    LOGV[Log-space Viterbi] --> STABLE[Stable decoding]

You implement the scaled forward first, then the log-space alternative, then log-space Viterbi. Each routine stays under thirty lines and runs with numpy only. No hmmlearn appears until Part II.

[!NOTE] Scaling and log-space solve the same underflow problem two ways. Scaling keeps probabilities normalized per step; log-space turns products into sums.

Scaling the Forward Pass

Scaling normalizes each day’s forward scores to sum to one. You store each day’s normalizer, called the scaling constant, before dividing. The log-likelihood is then the sum of the logs of those constants.

The relationship in words: the log of the total likelihood equals the sum over days of the log of each day’s scaling constant. In code the loop banks one constant per day.

import numpy as np

start_prob = np.array([0.6, 0.4])
trans_mat = np.array([[0.7, 0.3], [0.4, 0.6]])
emit_mat = np.array([[0.1, 0.4, 0.5], [0.6, 0.3, 0.1]])
obs_seq = [0, 2]  # Walk then Clean, shared toy sequence

num_steps = len(obs_seq)
fwd_scaled = np.zeros((num_steps, 2))
scale_const = np.zeros(num_steps)
# Day 1: raw scores, then bank their sum and normalize.
raw_day1 = start_prob * emit_mat[:, obs_seq[0]]
scale_const[0] = raw_day1.sum()
fwd_scaled[0] = raw_day1 / scale_const[0]
# Day 2: same rescale step keeps values near 1.0.
raw_day2 = (fwd_scaled[0] @ trans_mat) * emit_mat[:, obs_seq[1]]
scale_const[1] = raw_day2.sum()
fwd_scaled[1] = raw_day2 / scale_const[1]
log_likelihood = float(np.log(scale_const).sum())
print("scaled total:", np.exp(log_likelihood))
print("log-likelihood:", log_likelihood)

The worked numbers confirm Module 01. Day 1 raw scores are 0.06 and 0.24, so the first constant is 0.30. Day 2 raw-from-scaled scores are 0.23 and 0.054, so the second constant is 0.284. Their product is 0.30 times 0.284 = 0.0852, matching the hand total.

That product means scaling changes nothing mathematically; it only keeps intermediate values near 1.0. Your likelihood stays exact while your floats stay comfortable.

Log-Space Viterbi

Log-space replaces every multiplication with an addition of logs. You convert the start, transition, and emission tables once, then add instead of multiply. Underflow becomes impossible because even tiny probabilities are just moderate negative numbers.

The relationship in words: the log of a product equals the sum of the logs, so the Viterbi max over products becomes a max over sums. In code the max moves into log-space directly.

import numpy as np

start_prob = np.array([0.6, 0.4])
trans_mat = np.array([[0.7, 0.3], [0.4, 0.6]])
emit_mat = np.array([[0.1, 0.4, 0.5], [0.6, 0.3, 0.1]])
obs_seq = [0, 2]  # Walk then Clean

# Convert once: products below become additions here.
log_start = np.log(start_prob)
log_trans = np.log(trans_mat)
log_emit = np.log(emit_mat)
vit_log = np.zeros((2, 2))
vit_log[:, 0] = log_start + log_emit[:, obs_seq[0]]
for state_now in range(2):
    # Best predecessor in log-space, then add emission.
    prev_best = max(
        vit_log[pred, 0] + log_trans[pred, state_now]
        for pred in range(2)
    )
    vit_log[state_now, 1] = prev_best + log_emit[state_now, 1]
print("log-scores day 2:", vit_log[:, 1])
print("best prob:", np.exp(vit_log[:, 1].max()))

The worked numbers agree with Module 01. The day-2 log-scores are log(0.048) = -3.037 for Rainy and log(0.0144) = -4.241 for Sunny. Exponentiating the max returns 0.048, the same Viterbi winner. That match means your log-space port is faithful.

[!WARNING] Never mix log and linear values in one expression. You convert all three tables up front, or a stray raw probability silently corrupts the path.

Testing Your Implementation

You demonstrate the underflow danger first, then prove your engine correct. The snippet below shows what a 200-step joint probability looks like in float64. A per-step probability near 0.05 compounds fast.

import numpy as np

per_step_prob = 0.05
for horizon in [50, 200, 250]:
    # Repeated multiplication is what the raw forward does.
    joint_prob = per_step_prob ** horizon
    print(horizon, joint_prob)
print("float64 floor near 5e-324; min normal 2.2e-308")

The expected output prints about 8.9e-66 for 50 steps and about 6.2e-261 for 200 steps. At 250 steps Python prints exactly 0.0, a true underflow since the real value near 10^-326 sits below the smallest subnormal. That zero means any raw likelihood past that horizon is garbage, while the scaled and log-space versions keep working.

Your correctness proof enumerates all 16 state paths of a tiny T=4, K=2 example and compares the sum against your scaled forward. Both must print the same number.

import numpy as np
import itertools

start_prob = np.array([0.6, 0.4])
trans_mat = np.array([[0.7, 0.3], [0.4, 0.6]])
emit_mat = np.array([[0.1, 0.4, 0.5], [0.6, 0.3, 0.1]])
obs_seq = [0, 1, 2, 0]  # Walk Shop Clean Walk

# Brute force: score each of the 16 paths explicitly.
brute_total = 0.0
for path in itertools.product([0, 1], repeat=4):
    path_prob = start_prob[path[0]] * emit_mat[path[0], obs_seq[0]]
    for step in range(1, 4):
        prev, now = path[step - 1], path[step]
        path_prob *= trans_mat[prev, now] * emit_mat[now, obs_seq[step]]
    brute_total += path_prob

# Scaled forward on the same observations.
fwd = start_prob * emit_mat[:, obs_seq[0]]
scales = [fwd.sum()]
fwd = fwd / scales[0]
for step in range(1, 4):
    fwd = (fwd @ trans_mat) * emit_mat[:, obs_seq[step]]
    scales.append(fwd.sum())
    fwd = fwd / scales[-1]
scaled_total = float(np.prod(scales))
print("brute force:", brute_total)
print("scaled fwd :", scaled_total)
assert abs(brute_total - scaled_total) < 1e-12
print("engine matches brute force")

The expected output shows both lines printing the same number (about 0.0090888) before the confirmation. If they differ beyond 1e-12, your transition indexing or emission lookup is wrong. You do not proceed to Baum-Welch until this passes.

Hands-On Project

Extend the brute-force check to the shared 2-day Walk-then-Clean sequence with all four paths enumerated by hand. Write each path probability explicitly, sum them, and confirm 0.0852. Your reproduction card records the shared matrices, numpy only, and seed not applicable.

Your acceptance check has two parts. Your hand-enumerated four-path sum must equal your scaled-forward total to twelve decimals. Your log-space Viterbi best probability must equal the max single-path probability from your enumeration. If the Viterbi number exceeds the total, you have a bug, because no single path can beat the sum over all paths.

[!TIP] Print each of the four path probabilities separately. The largest one must equal 0.048 from Module 01.

Key Takeaways

  • Raw forward products decay geometrically and underflow to zero on long sequences.
  • At 0.05^250 float64 prints 0.0; scaling and log-space avoid this entirely.
  • Scaling banks one constant per day; the likelihood is the product of the constants.
  • Log-space turns Viterbi products into additions with identical winners.
  • Your scaled forward must match brute-force enumeration exactly on tiny examples.
  • No single Viterbi path can ever outscore the total over all paths.
  • You earn the right to use hmmlearn only after this from-scratch engine passes.

References

  • Mechanics: L. R. Rabiner, “A tutorial on hidden Markov models,” Proc. IEEE 77(2), 1989, scaling discussion, doi:10.1109/5.18626 — https://courses.physics.illinois.edu/ece417/fa2017/rabiner89.pdf, retrieved 2026-09-17.
  • Background: hmmlearn developer documentation on numerical stability conventions — https://hmmlearn.readthedocs.io/, retrieved 2026-09-17.
  • Exclusion: pre-0.x pomegranate HMM tutorials were excluded because their API predates the verified PyTorch rewrite and their code no longer runs.