Mixed-precision arithmetic

GPUs and TPUs run several times faster on 16-bit floating point than on 32-bit, and modern hardware has dedicated tensor cores for it. Mixed-precision uses 16-bit storage and compute where possible and falls back to 32-bit only where stability demands it[1]. Knowing which 16-bit format to use and where the failure modes hide is the difference between a 2x speedup and silent training divergence.

Three formats

fp16 (IEEE-754 half precision) has a tiny dynamic range: values above 65,504 or below about $6 \times 10^{-8}$ saturate. Gradients in deep networks routinely escape this range, so they need rescaling.

bf16 (bfloat16) keeps fp32's 8-bit exponent and gives up precision in the mantissa. It has the same range as fp32, so it does not overflow during training, at the cost of about three significant decimal digits of precision. This range advantage made bf16 the default for transformer training.

Where each format fails

Low precision (~$10^{-3}$ for fp16, $10^{-2}$ for bf16) shows up in three characteristic ways:

• Vanishing updates. When the weight magnitude $|w|$ is much larger than the update magnitude $|\eta g|$, the sum $w + \eta g$ rounds back to $w$ in low precision, and the update is lost. The standard fix is keeping a master copy of the weights in fp32, applying gradient updates to it, and casting to low precision only for the forward and backward passes.
• Lossy reductions. Summing many small numbers in fp16 underestimates the true sum, because each addition rounds. Reductions like batch normalization statistics and softmax denominators should accumulate in fp32 even when inputs are fp16.
• Overflow in fp16. Activations and gradients can exceed 65,504 (often the attention logits before softmax). Computing in bf16 sidesteps this; for fp16, loss scaling does.

Loss scaling

Loss scaling lifts small gradients out of the fp16 underflow zone before the backward pass. Multiply the loss by a large factor $S$ so that gradients are scaled by $S$ end-to-end (chain rule). After backprop, divide gradients by $S$ before the optimizer step:

The intermediate fp16 values stay in the representable range; the final unscaled gradient is recovered exactly modulo rounding. Dynamic loss scaling monitors for gradient overflow (any inf or NaN) and halves $S$ when it occurs, doubles $S$ after a clean run of steps: an automatic search for the largest safe scale[2].

Check it on arrays

import numpy as np

print("fp16  eps:", np.finfo(np.float16).eps,  " max:", np.finfo(np.float16).max)
print("fp32  eps:", np.finfo(np.float32).eps,  " max:", np.finfo(np.float32).max)

# fp16 overflow: gradients of moderate logits can blow up
x = np.float16(1e4)
print("x * 10 in fp16:", x * np.float16(10))                # inf

# Vanishing update: tiny gradient added to a moderate weight is lost in fp16
w  = np.float16(1.0)
g  = np.float16(1e-4)
print("fp16: 1 + 1e-4 ==", float(w + g), "  lost ->", w + g == w)
print("fp32: 1 + 1e-4 ==", np.float32(1.0) + np.float32(1e-4))

# Loss scaling fix: scale gradient up, do reduction in fp16, then unscale
g_small = np.full(1000, 1e-5, dtype=np.float32)            # true gradient
naive = g_small.astype(np.float16).mean()                  # underflows to 0
print("naive fp16 mean:", float(naive))                    # 0.0
S = 1024.0
scaled = (S * g_small).astype(np.float16).mean() / S       # loss-scaled
print("loss-scaled fp16 mean:", float(scaled), "  truth:", g_small.mean())

# Reduction precision: sum many fp16 values directly vs accumulating in fp32
v = np.full(10_000, 1e-3, dtype=np.float16)
direct = v.sum(dtype=np.float16)                           # lossy
upcast = v.astype(np.float32).sum()                        # accurate
print("fp16 reduction:", float(direct), "  fp32 accumulator:", float(upcast))

fp16 overflows above 65,504 and silently drops tiny updates, loss scaling rescues the underflowing gradients, and accumulating sums in fp32 fixes lossy reductions.

Exercises

Exercise 1. fp16 has $10$ mantissa bits, so its relative precision near $1$ is about $2^{-10}$. Argue that an update $w + \eta g$ with $w = 1$ is lost in fp16 once $|\eta g|$ falls below roughly $2^{-11}$, and state the standard fix.

Exercise 2. Loss scaling multiplies the loss by $S$ before backprop and divides gradients by $S$ after. Show that the recovered gradient equals the true gradient (modulo rounding), and explain what range constraint picks the largest safe $S$.

Exercise 3 (code). Demonstrate a lossy fp16 reduction: accumulate $5000$ ones one at a time in an fp16 register and show the running sum stalls at $2048$, while an fp32 accumulator reaches $5000$.

import numpy as np

acc16 = np.float16(0.0)
for _ in range(5000):
    acc16 = np.float16(acc16 + np.float16(1.0))

acc32 = np.float32(0.0)
for _ in range(5000):
    acc32 = np.float32(acc32 + np.float32(1.0))

# Once the fp16 running sum reaches 2048, the gap to the next float is 2,
# so adding 1 rounds back and the accumulator freezes.
assert float(acc16) == 2048.0
assert float(acc32) == 5000.0
print("fp16 sequential sum:", float(acc16), "  fp32 sequential sum:", float(acc32))

This closes the A4 numerical-methods track; A5 turns from floating-point math to the data structures and algorithms the runtime is built on.