Logarithms and underflow

Why you need this

Scoring a 300-column level under a hidden Markov model means multiplying about 600 probabilities, each well under 1. The true product is smaller than any number a double can hold, so the computer reports 0 for every level and cannot rank them. Logarithms turn that product into a sum that lands near -700, an ordinary number. The forward algorithm, Viterbi, and Baum-Welch all run in log space or with rescaling for exactly this reason.

The idea

What a logarithm is

$\log x$ is the power you raise $e \approx 2.718$ to in order to get $x$: $e^{\log x} = x$ and $\log(e^y) = y$. This cluster uses the natural log (base $e$) unless it says otherwise. The log is defined only for $x > 0$.

The rules you need:

Rule	Statement	With numbers
product	$\log(ab) = \log a + \log b$	$\log(0.5 \times 0.1) = -0.693 + (-2.303) = -2.996$
quotient	$\log(a / b) = \log a - \log b$	$\log(0.5 / 0.1) = -0.693 + 2.303 = 1.609$
power	$\log(a^n) = n \log a$	$\log(0.1^{300}) = 300 \times (-2.303) = -690.8$
anchors	$\log 1 = 0$, $\log e = 1$

Two facts about signs. The log of a number in $(0, 1)$ is negative, and the closer to 0 the number, the more negative: $\log 0.5 \approx -0.693$, $\log 0.1 \approx -2.303$, $\log 0.01 \approx -4.605$. And the log is increasing: if $a < b$ then $\log a < \log b$. So comparing two logs gives the same answer as comparing the two originals, which is the whole reason log space is safe to work in.

Underflow

A double-precision number (JavaScript's only number type) has a smallest normal value of about $2.2 \times 10^{-308}$. Below that the format enters denormals, which trade precision for range and bottom out near $4.9 \times 10^{-324}$. Anything smaller becomes exactly 0. That is underflow.

Now multiply 300 probabilities each around 0.1. The product is $0.1^{300} = 10^{-300}$, already in the region where a few more factors push it to zero. A realistic level has a transition and an emission factor per column, so 600 factors, and the product is well below $10^{-324}$: the computer says 0. Every candidate level scores 0. Ties everywhere, and no way to say which level was more likely.

In log space the same product is $300 \log 0.1 \approx -690.8$. Nothing special happens near that number; doubles handle magnitudes up to about $10^{308}$ with room to spare.

Comparing two paths in log space

Suppose two candidate paths through a chain have per-step probabilities

• path P: 0.5, 0.3, 0.2, 0.4
• path Q: 0.6, 0.1, 0.5, 0.5

Their log-probabilities are sums:

$$\log P = -0.693 - 1.204 - 1.609 - 0.916 = -4.423$$ $$\log Q = -0.511 - 2.303 - 0.693 - 0.693 = -4.200$$

Q is larger (less negative), so Q is more probable. Check with the raw products: P is 0.012 and Q is 0.015. The difference of logs is the log of the ratio, $-4.200 - (-4.423) = 0.223$, and $e^{0.223} \approx 1.25$: Q is 1.25 times as likely as P. With four steps you could have multiplied; with four hundred you could not.

Costs instead of probabilities

Negate the log and you get a non-negative cost: $-\log 0.5 = 0.693$, $-\log 0.1 = 2.303$. Rare steps cost more. A product of probabilities becomes a sum of costs, and "most probable path" becomes "cheapest path". That is why the Viterbi algorithm looks exactly like a shortest-path search on a grid.

Adding in log space: log-sum-exp

Multiplying is easy in log space (add the logs). Adding is not: there is no simple rule for $\log(a + b)$. You need it whenever you sum probabilities, and the forward algorithm does that at every step.

The naive route, exponentiate back, add, take the log, underflows: $e^{-1000}$ is 0 on a machine. The fix is to pull out the largest term first. With $\ell_k = \log a_k$ and $m = \max_k \ell_k$,

$$\log \sum_k e^{\ell_k} = m + \log \sum_k e^{\ell_k - m}$$

Every exponent $\ell_k - m$ is at most 0, and the largest is exactly 0, so the sum inside contains a 1 and can never underflow to 0 or overflow.

Three numbers: $\ell = (-1000, -1001, -1003)$. Then $m = -1000$, the shifted values are $0, -1, -3$, their exponentials are $1, 0.368, 0.050$, the sum is $1.418$, and $\log 1.418 = 0.349$. Answer: $-1000 + 0.349 = -999.651$. The direct computation returns $\log 0 = -\infty$.

Sanity check on numbers small enough to do both ways: $a = 0.5$, $b = 0.25$, so the answer should be $\log 0.75 = -0.288$. Logs are $-0.693$ and $-1.386$; $m = -0.693$; shifted $0$ and $-0.693$; exponentials $1$ and $0.5$; sum $1.5$; $\log 1.5 = 0.405$; result $-0.693 + 0.405 = -0.288$. Correct.

Worked example

Six columns, each with probability $p_t$ under some chain. Watch the running product shrink while the running log-sum just drifts down.

$t$	$p_t$	$\log p_t$	running product	running sum of logs
1	0.5	-0.693	0.5	-0.693
2	0.2	-1.609	0.1	-2.303
3	0.1	-2.303	0.01	-4.605
4	0.3	-1.204	0.003	-5.809
5	0.1	-2.303	0.0003	-8.112
6	0.2	-1.609	0.00006	-9.721

Check the last row: $\log 0.00006 = -9.721$. Six columns cost about 1.6 in log units each. At that rate 300 columns land near $-486$, a product of about $10^{-211}$, still representable but uncomfortably close; 500 columns land near $-810$, a product of $10^{-352}$, which a double rounds to 0. The log column has no such cliff.

function logSumExp(logs) {
  const m = Math.max(...logs);
  if (m === -Infinity) return -Infinity;   // every input was log(0)
  let s = 0;
  for (const l of logs) s += Math.exp(l - m);
  return m + Math.log(s);
}
console.log(logSumExp([-1000, -1001, -1003]));            // -999.6509877832318
console.log(Math.log(Math.exp(-1000) + Math.exp(-1001))); // -Infinity: the naive way underflows
console.log(Math.pow(0.1, 300), Math.pow(0.1, 330));      // 1.0000000000000166e-300 0
console.log(Number.MIN_VALUE, 2 ** -1022);                // 5e-324 2.2250738585072014e-308

The third line shows both halves of the story: 300 factors of 0.1 survive with a little rounding noise in the last digits, and 330 do not survive at all.

In a map generator

• Layout. Scoring a column sequence with the forward algorithm sums probabilities at every column, so it runs in log space with log-sum-exp or keeps a per-column scale factor. Recovering hidden sections with Viterbi needs only max and add, so it is pure cost arithmetic. Fitting tables with Baum-Welch accumulates expected counts in log space.
• Comparing levels. "How typical is this level of the training set" is a log-likelihood. Divide by the column count to compare levels of different lengths.
• Tuning. $-\log p$ is a "surprise" score. Printing the surprise of each column of a generated level shows you exactly where the chain did something rare.
• Regime over time. Regime sequences are short and rarely scored, so plain probabilities usually suffice there.

Common mistakes

• Taking $\log 0$. It is $-\infty$, and $-\infty - (-\infty)$ is NaN. Symptom: one impossible transition turns a whole table to NaN. Guard zeros, or floor probabilities at a tiny epsilon before taking logs.
• Multiplying first, logging after. Once the product is 0, its log is $-\infty$ for every candidate. Symptom: Viterbi always returns the first path it tried; every level ties.
• Mixing bases. One table in $\log_2$, another in natural log. Symptom: two models disagree by a constant factor of about 1.44 and nobody knows why.
• Summing with exp-add-log inside the forward algorithm. Symptom: the forward variable becomes 0 somewhere around column 150 and stays there.
• Comparing paths of different lengths by raw total. Every extra step adds a negative number, so longer always loses. Symptom: the generator prefers the shortest possible level.
• Feeding log-probabilities to a sampler as weights. They are negative. Symptom: the sampler returns index 0 forever or throws.

Cost

Taking the log of each entry of $A$ ($N \times N$, with $N$ states) and $B$ ($N \times M$, with $M$ observation symbols) once costs $O(N^2 + NM)$ time and the same memory, done at load. Log-sum-exp over $N$ terms is $O(N)$ with $N$ calls to exp. The forward algorithm performs one such sum per state per column, so $O(T N^2)$ exponentials for $T$ columns, a constant factor over plain multiplication. For $N$ under 30 you will not notice. The alternative, rescaling each column to sum to 1 and keeping the log of the scale factors, avoids the exponentials entirely and is the usual choice for long sequences.

Going further

• Forward algorithm: log-sum-exp or rescaling in action.
• Viterbi algorithm: costs, max, and add, no log-sum-exp needed.
• Baum-Welch: expected counts accumulated in log space.
• Entropy: the expected value of $-\log p$.
• The IEEE 754 double-precision format, for the exact limits quoted above.