Introduction
When working with probabilities, likelihoods, or any computation involving many small numbers, you’ll quickly run into floating-point underflow — where numbers become so small that Python rounds them to zero. The solution is to work in the log domain, where multiplication becomes addition and tiny numbers become manageable negative values.
The Problem: Floating-Point Underflow
Python’s 64-bit floats (IEEE 754 double precision) can represent numbers as small as about 5e-324. Below that, they underflow to zero:
import numpy as np
# Numbers near the underflow boundary
a = 1e-323
b = 1e-324 # underflows to 0.0!
print(a, b)
# => 1e-323 0.0
# Multiplying many small probabilities
probs = [0.1] * 500
product = 1.0
for p in probs:
product *= p
print(product) # => 0.0 (underflowed!)
This is a real problem in machine learning, NLP, and statistics — computing the probability of a sequence of events (like words in a sentence) involves multiplying many small probabilities together.
The Solution: Log Domain Arithmetic
The key insight: log(a * b) = log(a) + log(b)
Instead of multiplying probabilities, add their logarithms:
import numpy as np
# Instead of multiplying
probs = [0.1] * 500
log_probs = [np.log(p) for p in probs]
log_product = sum(log_probs)
print(log_product) # => -1151.29... (no underflow!)
print(np.exp(log_product)) # => 1e-500 (the actual value)
Log Arithmetic Rules
| Operation | Normal Domain | Log Domain |
|---|---|---|
| Multiply | a * b |
log(a) + log(b) |
| Divide | a / b |
log(a) - log(b) |
| Power | a^n |
n * log(a) |
| Add | a + b |
logsumexp([log(a), log(b)]) |
Practical Example: Language Model Probability
import numpy as np
from scipy.special import logsumexp
# Probability of each word in a sentence
word_probs = [0.05, 0.12, 0.08, 0.03, 0.15, 0.07]
# WRONG: underflows for long sequences
sentence_prob = 1.0
for p in word_probs:
sentence_prob *= p
print(f"Direct: {sentence_prob}") # => 1.2096e-08 (ok for 6 words, fails for 500)
# CORRECT: log domain
log_sentence_prob = sum(np.log(p) for p in word_probs)
print(f"Log domain: {log_sentence_prob}") # => -18.23...
print(f"Recovered: {np.exp(log_sentence_prob)}") # => 1.2096e-08
The logsumexp Trick
Adding probabilities in log space is tricky because log(a + b) ≠ log(a) + log(b). Use logsumexp:
from scipy.special import logsumexp
# Problem: log(1.25e-500 + 3.54e-500)
# Direct computation underflows:
a = 1.25e-500
b = 3.54e-500
print(np.log(a + b)) # => -inf (both underflow to 0.0)
# Log domain solution:
log_a = np.log(1.25) - 500 * np.log(10) # log(1.25e-500)
log_b = np.log(3.54) - 500 * np.log(10) # log(3.54e-500)
log_sum = logsumexp([log_a, log_b])
print(log_sum) # => -1149.726...
# Recover the actual sum
actual_sum = np.exp(log_sum + 500 * np.log(10))
print(actual_sum) # => 4.79 (= 1.25 + 3.54)
How logsumexp Works
The log-sum-exp trick avoids overflow/underflow by factoring out the maximum:
log(sum(exp(x_i))) = max(x) + log(sum(exp(x_i - max(x))))
Since x_i - max(x) <= 0, all exponentials are in [0, 1] — no overflow. The maximum is added back exactly.
def logsumexp_manual(x):
c = np.max(x)
return c + np.log(np.sum(np.exp(x - c)))
# Test
x = np.array([1000, 1001, 1002])
print(logsumexp_manual(x)) # => 1002.408...
print(logsumexp(x)) # => 1002.408... (same)
2D Arrays: Row-wise and Column-wise
from scipy.special import logsumexp
a = np.array([[1, 2, 3],
[4, 5, 6]])
# Sum along rows (axis=1)
print(logsumexp(a, axis=1))
# => [3.40760596 6.40760596]
# Row 0: log(exp(1) + exp(2) + exp(3)) = 3.4076...
# Row 1: log(exp(4) + exp(5) + exp(6)) = 6.4076...
# Sum along columns (axis=0)
print(logsumexp(a, axis=0))
# => [4.31326169 5.31326169 6.31326169]
Softmax in Log Space
Softmax is exp(x_i) / sum(exp(x)). Computing it directly overflows for large inputs:
# WRONG: overflows for large logits
def softmax_naive(x):
return np.exp(x) / np.sum(np.exp(x))
logits = np.array([1000, 1001, 1002])
print(softmax_naive(logits)) # => [nan nan nan] (overflow!)
# CORRECT: numerically stable
def log_softmax(x):
return x - logsumexp(x)
def softmax(x):
return np.exp(log_softmax(x))
print(softmax(logits)) # => [0.09003057 0.24472847 0.66524096]
print(softmax(logits).sum()) # => 1.0
Normalizing Log-Probabilities
# Unnormalized log-probabilities (e.g., from a model)
log_scores = np.array([-1.0, -2.0, -0.5, -3.0, -1.5])
# Normalize: subtract log(sum(exp(scores)))
log_probs = log_scores - logsumexp(log_scores)
# Verify they sum to 1 in probability space
print(np.exp(log_probs).sum()) # => 1.0
print(log_probs)
# => [-1.56... -2.56... -1.06... -3.56... -2.06...]
Comparing Log-Probabilities
When comparing which of two sequences is more likely, compare log-probabilities directly:
# Which sentence is more likely?
log_prob_sentence1 = -45.3
log_prob_sentence2 = -38.7
# Compare directly in log space (no need to exponentiate)
if log_prob_sentence1 > log_prob_sentence2:
print("Sentence 1 is more likely")
else:
print("Sentence 2 is more likely") # => this one
Common Pitfalls
# Pitfall 1: log(0) = -inf
np.log(0) # => -inf (not an error, but propagates)
# Handle with np.log1p for values near 0
np.log1p(1e-300) # log(1 + 1e-300) — more stable than log(1e-300)
# Pitfall 2: exp of large positive numbers overflows
np.exp(1000) # => inf
# Pitfall 3: Forgetting to use logsumexp for log-domain addition
log_a = -500.0
log_b = -501.0
# WRONG
log_sum = np.log(np.exp(log_a) + np.exp(log_b)) # => -inf (underflow)
# CORRECT
log_sum = logsumexp([log_a, log_b]) # => -499.69...
Comments