Mixed-Precision Training

Mixed-precision training runs most of a network’s training arithmetic and storage in a 16-bit floating-point format, FP16 (IEEE half precision) or BF16 (bfloat16), and more recently in 8-bit FP8, while keeping the few steps that need more range or precision in 32-bit FP32. In NVIDIA’s guide the term means the combined use of different numerical precisions in one computation. Three techniques make it work: an FP32 master copy of the weights that receives the updates, loss scaling, which shifts small gradients into the range of FP16, and FP32 accumulation of the products inside matrix multiplications and convolutions. Frameworks automate it as automatic mixed precision (AMP); in PyTorch that is torch.autocast together with torch.amp.GradScaler. The neural network training article places the technique in the training workflow; this article explains the mechanism, the formats as they affect training, the framework implementation, FP8 and the evidence.

The technique is distinct from the mixed-precision quantization of a trained model for inference, which gives different layers different bit widths, and from running a trained model in FP16 at inference time.

Purpose

Training in FP32 spends memory, memory bandwidth and arithmetic on precision that most of the computation does not need. NVIDIA’s guide lists three gains from 16-bit formats: less memory, which allows larger models or larger mini-batches; less memory traffic, since half precision halves the bytes moved in layers limited by bandwidth; and faster arithmetic, since its GPUs offer up to 8 times more half-precision than single-precision throughput for math-limited layers. On a V100, the guide states, Tensor Cores speed up matrix multiplications and convolutions by up to 8 times in FP16 over FP32. He et al. (2019) quote the V100 at 14 TFLOPS in FP32 and over 100 TFLOPS in FP16.

Simply casting everything to FP16 does not train every network. FP16 has 5 exponent bits and 10 fraction bits: its normal numbers run from $2^{-14}$ to 65,504 and its subnormals down to $2^{-24}$, a dynamic range of 40 powers of two against 264 for FP32 (NVIDIA guide, §2.1). Small gradients round to zero, and small weight updates vanish against large weights. Mixed precision keeps the speed and memory gains of 16-bit arithmetic while restoring FP32 behaviour where these failures occur.

When to Use It

Mixed precision helps when training is limited by arithmetic on matrix multiplications and convolutions, by memory bandwidth, or by memory capacity, on hardware with fast low-precision units. Its benefit depends on the share of time spent in those operations. The guide’s example, an instance of Amdahl’s law: if half the training time is in Tensor Core operations and they become 5 times faster, the whole step becomes $1/(0.5 + 0.5/5) \approx 1.67$ times faster. Very small linear and convolution layers gain little, since there is not enough arithmetic to use Tensor Cores fully, and the guide advises first finding whether a workflow is limited by data input or by computation.

Which form to use depends on the hardware. FP16 needs loss scaling for many networks. BF16 has FP32’s range, so, in our reading, loss scaling is not needed. Google’s TPUs multiply in BF16 and accumulate in FP32 by default, and PyTorch’s own CPU training example uses BF16 autocasting with no gradient scaler. FP8 needs hardware that supports it and a scale for each tensor (see Variants).

Precision is not free. A network whose weights were trained in BF16 may not train in FP16: PyTorch’s GradScaler warns that models pretrained in bf16 are often incompatible with AMP in FP16 because of the different dynamic ranges.

How It Works

The quantization article’s section on number formats tabulates the bits, range and machine epsilon of FP32, FP16, BF16 and the two FP8 encodings. Here the question is how those limits act on the quantities of training. NVIDIA’s blog post on the technique sorts the tensors of training into four kinds: activations, activation gradients, weights and weight gradients. In the post’s experience, activations, weights and weight gradients fit in the range of FP16. Activation gradients of some networks fall below it.

FP32 master weights

Each gradient descent step adds an update $-\eta g$ to every weight. Weight gradients are often much smaller than the weights, especially after multiplication by the learning rate (Micikevicius, 2017). A floating-point addition aligns the exponents of its operands, and when the update is less than half the spacing between the weight and its neighbouring representable number, half a unit in the last place (ulp), the sum rounds back to the weight. For an FP16 weight between 1 and 2, the spacing is $2^{-10}$, so every update smaller than $2^{-11} \approx 4.9 \times 10^{-4}$ is lost, however many steps repeat it. In BF16, with 7 fraction bits, the threshold there is $2^{-8} \approx 3.9 \times 10^{-3}$.

The remedy keeps a master copy of the weights in FP32 and applies the updates to it. Each iteration rounds the master weights to a 16-bit copy for the forward and backward passes, converts the weight gradients to FP32, and updates the master copy (Micikevicius, 2017). Small updates then accumulate in FP32 until they move the 16-bit copy. The blog post notes that not every network needs this; recurrent networks often did.

Loss scaling

Activation gradients are the vulnerable tensor. NVIDIA’s guide shows a histogram of the activation gradients of the SSD detector (VGG-D backbone) collected across all layers during FP32 training. Its caption gives 66.8% of the values as zero; the text says that converting to FP16 turns 31% of the values into zeros, leaving only 5.3% nonzero, which made training diverge. As the guide states them, these three shares do not add up. Most of FP16’s range above these values was unused, so the gradients can be shifted up. By the chain rule, multiplying the loss by a factor $S$ before backpropagation multiplies every gradient in the backward pass by $S$, at the cost of one multiplication:

$$\frac{\partial (S\,\mathcal{L})}{\partial \theta} = S\,\frac{\partial \mathcal{L}}{\partial \theta}.$$

The weight gradients are divided by $S$ before the update. For SSD, multiplying by 8 was enough to match FP32 accuracy, and multiplying by 32K would have recovered all but 0.1% of the values lost to zero without overflow (guide, §2.3). The guide calls it simplest to unscale right after the backward pass, before gradient clipping or other gradient computations, so that no hyperparameter such as the clipping threshold or weight decay has to change.

The factor can be chosen two ways (guide, §2.3.2):

• Static scaling fixes $S$. Across the networks NVIDIA trained, the required factors ranged from 8 to 32K, and many networks needed none. When gradient statistics are known, $S$ is chosen so that $S$ times the largest absolute gradient stays below 65,504.
• Dynamic scaling starts with a large $S$ and adjusts it during training. If a gradient contains an infinity or NaN, the step is skipped and $S$ is reduced; if no overflow has occurred for $N$ iterations, $S$ is increased. NVIDIA trained networks with $N = 2000$, a growth factor of 2 and a reduction factor of 0.5, and found that skipping updates infrequently does not require changing the training schedule.

Overflow is detected by checking the weight gradients for an infinity or NaN (guide, §2.3.2).

FP32 accumulation and FP32 operations

A matrix product or convolution sums many products. Volta Tensor Cores compute $D = AB + C$ on $4 \times 4$ FP16 matrices $A$ and $B$, with $C$ and $D$ in FP16 or FP32 (guide, §2.2). NVIDIA found accumulation into single precision critical to good training results; the accumulated values are converted to half precision before being written to memory (Micikevicius, 2017). The product of two FP16 numbers has at most 22 significant bits, so it fits exactly in FP32; the sum is where precision is lost.

Some operations stay in FP32 altogether. The guide names values computed by large reductions, such as the mean and variance of batch normalization and the softmax; batch normalization may still read and write FP16 tensors. NVIDIA’s AMP sorts operations into three lists: an allow list of operations that benefit from Tensor Cores (convolutions, fully connected layers), a deny list of operations that may overflow FP16 or need FP32 precision (large reductions, cross-entropy loss, L1 loss, exponentials), and an infer list of operations safe in either format (element-wise additions, multiplication by a constant).

BF16

BF16 has 1 sign bit, 8 exponent bits and 7 mantissa bits with bias 127 (OFP8 specification, §4.1), the same exponent field as FP32. Its range equals FP32’s (Google Cloud documentation), so gradients that underflow FP16 survive in BF16 without scaling, but its spacing is eight times coarser than FP16’s, as the quantization article notes, and by the half-ulp rule above small updates stall sooner. On TPUs, conversion from FP32 to BF16 rounds to nearest even and overflows to infinity, subnormals are flushed to zero, and a conversion back to FP32 pads the mantissa with zeros, since FP32 represents every BF16 value exactly.

Memory

By our arithmetic, the master copy means that mixed precision does not halve the memory of the weights: 2 bytes for the 16-bit copy plus 4 for the FP32 master weights is 6 bytes per parameter, against 4 in FP32. In our reading, the savings come from activations, which a network stores for the backward pass, and, in recipes that keep gradients in 16 bits, from the gradients; under PyTorch autocast the gradients of FP32 parameters are themselves FP32. In such a 16-bit-gradient recipe, with Adam and its two moment estimates kept in FP32, the model states take 16 bytes per parameter, the accounting the distributed training article uses for sharding.

How to Apply It

Automatic mixed precision in PyTorch

PyTorch’s AMP follows the recipe above, with one difference of form: the model’s parameters themselves stay in FP32 and serve as the master weights. The following describes the PyTorch 2.14.1 source.

• torch.autocast(device_type, dtype=None) wraps the forward pass and the loss computation. Inside it, each operation runs in an operation-specific type. The default is torch.float16 on CUDA and torch.bfloat16 on CPU. The docstring warns against calling .half() or .bfloat16() on the model or its inputs, and recommends exiting autocast before backward(); backward operations run in the type autocast chose for the matching forward operation.
• The operation lists live in ATen’s autocast_mode.h and .cpp. On CUDA, convolutions, matmul, mm, bmm, addmm, linear, einsum, recurrent cells and scaled_dot_product_attention are cast to the lower-precision type. Operations cast to FP32 include exp, log, pow, reciprocal, rsqrt, softplus, layer_norm, group_norm, rms_norm, most losses (NLL, L1, smooth L1, MSE, KL divergence, binary_cross_entropy_with_logits), cdist, logsumexp and the upsampling operators. softmax, log_softmax, sum, prod, cumsum and the vector and matrix norms run in FP32 unless the caller sets an output type. Operations with several inputs, such as addcmul, index_put and scatter_add, run in the widest input type. Inputs are cast before autograd records them, so the lower-precision ops save their inputs for the backward pass in 16 bits; the source comments note that this can significantly reduce a model’s memory footprint.
• binary_cross_entropy is banned under autocast. The error message explains that many models put a sigmoid right before it and recommends binary_cross_entropy_with_logits, which is safe to autocast.
• torch.amp.GradScaler implements dynamic loss scaling with the defaults init_scale $= 2^{16}$, growth_factor $= 2$, backoff_factor $= 0.5$ and growth_interval $= 2000$. scaler.scale(loss) multiplies the loss by the current scale. scaler.step(optimizer) unscales the gradients, checks them for infinities and NaNs, and skips optimizer.step() if any are found; scaler.update() then multiplies the scale by the backoff factor, or by the growth factor after 2,000 consecutive unskipped steps. The docstring expects skipped steps in the first few iterations, while the scale calibrates, and rarely afterwards. scaler.unscale_(optimizer) unscales early, for clipping the true gradients. The scaler refuses to unscale FP16 gradients, which is why the parameters must stay in FP32.

A training loop then reads as follows, in the docstrings’ pattern. The snippet follows the PyTorch 2.14.1 API and was not run for this article, since it needs a CUDA GPU; model, loss_fn, optimizer and loader are the usual training objects.

scaler = torch.amp.GradScaler("cuda")
for images, targets in loader:
    optimizer.zero_grad()
    with torch.autocast(device_type="cuda", dtype=torch.float16):
        loss = loss_fn(model(images), targets)
    scaler.scale(loss).backward()
    scaler.unscale_(optimizer)                      # optional: clip the true gradients
    torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0)
    scaler.step(optimizer)                          # skipped if a gradient is inf or NaN
    scaler.update()

With dtype=torch.bfloat16 the scaler can be left out. The guide adds practical rules for Tensor Cores: mini-batch sizes, linear layer widths and convolution channel counts that are multiples of 8.

Example: reduced precision in an emulated training run

The example emulates the formats in NumPy, which shows what the hardware hides. It trains logistic regression on synthetic $16 \times 16$ images by full-batch gradient descent and rounds every stored tensor to the format under test. Logistic regression is convex, so its exact optimum can be found by Newton’s method and every run can be measured against it.

Each image holds a blurred blob (Gaussian, $\sigma = 2$ pixels) at a random position plus Gaussian pixel noise ($\sigma = 0.5$); the label says whether the blob’s centre lies right of the vertical midline. There are 4,096 training and 2,048 test images. A rule that needs no learning, comparing the brightness of the two image halves, checks that the labels can be predicted from the pixels. The noise makes the classes overlap, so the optimum is finite (its loss, printed below, is well above zero). The loss is the mean cross-entropy; its gradient with respect to the logit of image $i$, the activation gradient, is $dz_i = (p_i - y_i)/N$, and the weight gradient is $A^\top dz$, where $A$ holds the images with a bias column. The step size is $1/L$, with $L$ the largest eigenvalue of $A^\top A/(4N)$, and every run takes 3,000 steps from zero.

One function, quantize, rounds to any binary floating-point format, given its significand bits, smallest normal exponent and largest value, with round-to-nearest-even and subnormals. The script checks it against NumPy’s own float16 conversion; the same function is used for BF16 and for the FP8 formats E5M2 and E4M3 of the OFP8 specification. Additions in FP16 are emulated as an exact float64 addition followed by one rounding.

The script has three parts:

• Part B, training. It trains in six regimes, which separate the master copy from the loss scale: FP32 throughout; FP16 weights updated in FP16; FP16 arithmetic with an FP32 master copy, with and without loss scaling; BF16 weights updated in BF16; and BF16 arithmetic with an FP32 master copy. In the reduced-precision regimes, inputs, weights, logits, activation gradients and weight gradients are stored in the format, while products are accumulated in FP32 and the sigmoid is computed in FP32. The loss scale follows GradScaler’s rules and defaults. The printout gives the excess training loss over the optimum, the test accuracy, the share of nonzero updates that left the stored weight unchanged at the last step, the number of skipped steps and the final scale.
• Part A, which gradients survive. At two points of the FP32 run, after 10 and 3,000 steps, it rounds the activation gradients to each format and reports the share of nonzero values that become zero or infinite, and the relative error of the resulting weight gradient. FP16 is swept over loss scales $S = 2^0, 2^4, \dots, 2^{32}$. BF16 is unscaled. FP8 is shown without scaling and with a per-tensor scale that maps the largest absolute value to the format’s maximum. The weight gradient is computed from inputs in the row’s 16-bit format (BF16 for the FP8 rows), accumulated in FP32 and stored in the 16-bit format.
• Part C, accumulation. It sums the 262,144 pixels of scikit-image’s bundled camera image (values in $[0, 1]$, rounded to FP16 first) and the 4,096 terms of each weight-gradient entry at step 3,000, sequentially and pairwise, in FP16 and in FP32. Pairwise is printed next to sequential because the order of summation is a choice that the error depends on; in our reading, GPU kernels sum in blocks and trees rather than in one long chain.

The script runs with NumPy 2.3.5 and scikit-image 0.25.2, in about 5 seconds.

import os
for v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"):
    os.environ[v] = "1"                                   # fixed summation order in BLAS
import numpy as np
from skimage import data

def quantize(x, p, emin, xmax, saturate=False):
    """Round to a binary float format: p significand bits, smallest normal exponent emin,
    largest finite value xmax; round to nearest, ties to even, subnormals kept."""
    x = np.asarray(x, dtype=np.float64)
    _, e = np.frexp(x)                                    # |x| = m 2^e with m in [0.5, 1)
    step = np.ldexp(1.0, np.maximum(e - 1, emin) - (p - 1))
    r = np.round(x / step) * step                         # np.round: half to even
    return np.where(np.abs(r) > xmax, np.copysign(xmax if saturate else np.inf, x), r)

FP16, BF16 = (11, -14, 65504.0), (8, -126, (2 - 2.0 ** -7) * 2.0 ** 127)
E5M2, E4M3 = (3, -14, 57344.0), (4, -6, 448.0)
f16 = lambda x: quantize(x, *FP16)
bf16 = lambda x: quantize(x, *BF16)

# Checks: the rounding matches NumPy's float16 cast; updates below half an ulp vanish.
rng = np.random.default_rng(1)
t = rng.choice([-1, 1], 200000) * 2.0 ** rng.uniform(-30, 17, 200000)
t = np.concatenate([t, 1 + (2 * np.arange(1, 500) + 1) * 2.0 ** -11])   # exact ties
with np.errstate(over="ignore"):                         # values past 65504 overflow to inf
    print("quantize vs np.float16 mismatches:", np.sum(f16(t) != t.astype(np.float16)))
print("fp16: 1 + 2^-12 ->", f16(1 + 2.0 ** -12), " 1 + 2^-11 ->", f16(1 + 2.0 ** -11),
      " 1 + 3*2^-12 ->", f16(1 + 3 * 2.0 ** -12), "  fp32: 1 + 2^-12 ->", np.float32(1) + np.float32(2.0 ** -12))

# Synthetic images: a blurred blob; label 1 if its centre is right of the midline.
def make_images(n, rng, size=16, sigma=2.0, noise=0.5):
    cx, cy = rng.uniform(2.0, size - 3.0, n), rng.uniform(2.0, size - 3.0, n)
    yy, xx = np.mgrid[0:size, 0:size]
    img = np.exp(-((xx - cx[:, None, None]) ** 2 + (yy - cy[:, None, None]) ** 2) / (2 * sigma ** 2))
    img += noise * rng.standard_normal(img.shape)
    return np.hstack([img.reshape(n, -1), np.ones((n, 1))]), (cx > (size - 1) / 2).astype(float)

rng = np.random.default_rng(0)
A, y = make_images(4096, rng)                             # last column: bias input
Ate, yte = make_images(2048, rng)
N = len(y)
half_rule = lambda X: X[:, :-1].reshape(-1, 16, 16)[:, :, 8:].sum((1, 2)) > X[:, :-1].reshape(-1, 16, 16)[:, :, :8].sum((1, 2))
print(f"class-1 share {y.mean():.3f} / {yte.mean():.3f}; right-half-brighter rule, test accuracy "
      f"{np.mean(half_rule(Ate) == yte):.4f}")

sigmoid = lambda z: 0.5 * (1 + np.tanh(z / 2))
loss = lambda w: np.mean(np.logaddexp(0, A @ w) - y * (A @ w))

w = np.zeros(A.shape[1])                                  # Newton's method: the exact optimum
for _ in range(40):
    p = sigmoid(A @ w); g = A.T @ (p - y) / N
    d = np.linalg.solve((A * (p * (1 - p))[:, None]).T @ A / N, g); s = 1.0
    while loss(w - s * d) > loss(w) - 1e-4 * s * (g @ d) and s > 1e-10: s /= 2
    w -= s * d
L_star = loss(w)
eta = 4 / np.linalg.eigvalsh(A.T @ A / N).max()           # step 1/L
print(f"optimum loss {L_star:.6f}, step size {eta:.4f}")

# ---- Part B: gradient descent with weights and arithmetic in different formats --------------
def train(fmt, master, scaling, T=3000, keep=()):
    q = (lambda x: np.asarray(x, np.float32).astype(np.float64)) if fmt is None else fmt
    Aq = q(A).astype(np.float32)                          # inputs stored in the compute format
    w = np.zeros(A.shape[1])                              # float64 holder of the stored weights
    S, tracker, skipped, saved = (2.0 ** 16 if scaling else 1.0), 0, 0, {}
    for step in range(1, T + 1):
        wc = q(w).astype(np.float32)                      # compute copy of the weights
        z = q(Aq @ wc)                                    # products accumulated in float32
        p = sigmoid(z.astype(np.float32))
        dz = q(S * (p - y) / N)                           # activation gradient, stored in the format
        gs = q(Aq.T @ dz.astype(np.float32))              # weight gradient, stored in the format
        if not np.all(np.isfinite(gs)):                   # overflow: skip the step, halve S
            skipped += 1; S *= 0.5; tracker = 0; continue
        g = (gs.astype(np.float32) / np.float32(S)).astype(np.float64)
        if scaling:
            tracker += 1
            if tracker == 2000: S *= 2.0; tracker = 0
        if master or fmt is None:                         # float32 weights (master copy)
            u = np.float32(eta) * g.astype(np.float32)
            new = (w.astype(np.float32) - u).astype(np.float64)
        else:                                             # weights kept in the low format
            u = q(eta * g); new = q(w - u)
        lost = np.mean(new[u != 0] == w[u != 0]) if np.any(u != 0) else 0.0
        w = new
        if step in keep: saved[step] = w.copy()
    acc = np.mean(((Ate @ w) > 0) == yte)
    return w, loss(w) - L_star, acc, lost, skipped, S, saved

print("\nregime                          excess loss  test acc  lost updates  skipped  final S")
w32, *res, saved = train(None, True, False, keep=(10, 3000))
print(f"{'FP32':30s}  {res[0]:11.3e}  {res[1]:8.4f}  {res[2]:12.3f}  {res[3]:7d}  {res[4]:7g}")
for name, fmt, master, scaling in [("FP16 weights, loss scaling", f16, False, True),
                                   ("FP16 + FP32 master, scaling", f16, True, True),
                                   ("FP16 + FP32 master, no scaling", f16, True, False),
                                   ("BF16 weights", bf16, False, False),
                                   ("BF16 + FP32 master", bf16, True, False)]:
    _, *res, _ = train(fmt, master, scaling)
    print(f"{name:30s}  {res[0]:11.3e}  {res[1]:8.4f}  {res[2]:12.3f}  {res[3]:7d}  {res[4]:7g}")

# ---- Part A: which activation gradients survive the cast, against the loss scale -----------
for step in (10, 3000):
    dz = (sigmoid(A @ saved[step]) - y) / N; nz = dz != 0
    g_true = A.T @ dz
    print(f"\nafter {step} FP32 steps: |dz| from {np.abs(dz[nz]).min():.1e} to {np.abs(dz).max():.1e}")
    print("  dz format, scale     zero    inf   gradient rel. error")
    def row(label, dzq, S, q16):                          # q16: format of the inputs and the output
        with np.errstate(invalid="ignore"):              # inf * 0 when dz overflows
            gq = q16(q16(A).astype(np.float32).T @ dzq.astype(np.float32)) / S   # fp32 accumulation
        err = np.linalg.norm(gq - g_true) / np.linalg.norm(g_true)
        print(f"  {label:18s} {np.mean(dzq[nz] == 0):6.4f} {np.mean(np.isinf(dzq[nz])):6.4f}   {err:.2e}")
    for k in range(0, 33, 4):
        row(f"fp16, S = 2^{k}", f16(2.0 ** k * dz), 2.0 ** k, f16)
    row("bf16, S = 1", bf16(dz), 1.0, bf16)
    for name, fmt in (("E5M2", E5M2), ("E4M3", E4M3)):
        s = fmt[2] / np.abs(dz).max()                     # per-tensor scale: amax -> format max
        row(f"{name}, unscaled", quantize(dz, *fmt, saturate=True), 1.0, bf16)
        row(f"{name}, per-tensor", quantize(s * dz, *fmt, saturate=True), s, bf16)

# ---- Part C: accumulating sums in fp16 and in fp32 ------------------------------------------
def adder(fmt):                                           # one correctly rounded addition
    if fmt == np.float16:                                 # exact in float64, then one rounding
        def add(a, b):
            with np.errstate(over="ignore"):              # sums past 65504 become inf
                return (a.astype(np.float64) + b).astype(np.float16)
        return add
    return lambda a, b: a + b                             # native float32 addition

def seq_sum(x, fmt):
    x, add = x.astype(fmt), adder(fmt)
    acc = np.zeros(x.shape[1:], fmt)
    for v in x:
        acc = add(acc, v)
    return acc.astype(np.float64)

def pair_sum(x, fmt):
    x, add = x.astype(fmt), adder(fmt)
    while len(x) > 1:
        if len(x) % 2: x = np.concatenate([x, np.zeros((1,) + x.shape[1:], fmt)])
        x = add(x[0::2], x[1::2])
    return x[0].astype(np.float64)

cam = f16(data.camera().astype(np.float64) / 255).ravel()   # 262,144 pixels in [0, 1]
exact = cam.sum()
print(f"\ncamera pixel sum: exact {exact:.4f}")
for name, fmt in (("fp16", np.float16), ("fp32", np.float32)):
    for how, fn in (("sequential", seq_sum), ("pairwise", pair_sum)):
        s = fn(cam, fmt)
        print(f"  {name} {how:10s} {s:12.4f}  rel. error {abs(s - exact) / exact:.2e}")

dz = f16(2.0 ** 16 * (sigmoid(A @ saved[3000]) - y) / N)   # scaled fp16 gradients, step 3000
terms = f16(f16(A) * dz[:, None])                           # fp16 products, rounded to fp16
exact = terms.sum(axis=0)                                   # exact (float64) sum of the same terms
print("weight gradient A^T dz (4,096 terms per entry):")
for name, fmt in (("fp16", np.float16), ("fp32", np.float32)):
    for how, fn in (("sequential", seq_sum), ("pairwise", pair_sum)):
        s = fn(terms, fmt)
        print(f"  {name} {how:10s} rel. error {np.linalg.norm(s - exact) / np.linalg.norm(exact):.2e}")

quantize vs np.float16 mismatches: 0
fp16: 1 + 2^-12 -> 1.0  1 + 2^-11 -> 1.0  1 + 3*2^-12 -> 1.0009765625   fp32: 1 + 2^-12 -> 1.0002441
class-1 share 0.491 / 0.502; right-half-brighter rule, test accuracy 0.9375
optimum loss 0.079612, step size 0.8550

regime                          excess loss  test acc  lost updates  skipped  final S
FP32                              1.214e-04    0.9365         0.000        0        1
FP16 weights, loss scaling        2.598e-03    0.9395         1.000        0   131072
FP16 + FP32 master, scaling       1.213e-04    0.9365         0.000        0   131072
FP16 + FP32 master, no scaling    1.239e-04    0.9365         0.000        0        1
BF16 weights                      1.778e-02    0.9468         1.000        0        1
BF16 + FP32 master                1.230e-04    0.9365         0.000        0        1

after 10 FP32 steps: |dz| from 4.1e-07 to 2.2e-04
  dz format, scale     zero    inf   gradient rel. error
  fp16, S = 2^0      0.0000 0.0000   2.64e-04
  fp16, S = 2^4      0.0000 0.0000   2.36e-04
  fp16, S = 2^8      0.0000 0.0000   2.34e-04
  fp16, S = 2^12     0.0000 0.0000   2.34e-04
  fp16, S = 2^16     0.0000 0.0000   2.34e-04
  fp16, S = 2^20     0.0000 0.0000   2.34e-04
  fp16, S = 2^24     0.0000 0.0000   inf
  fp16, S = 2^28     0.0000 0.0000   inf
  fp16, S = 2^32     0.0000 0.5750   nan
  bf16, S = 1        0.0000 0.0000   2.01e-03
  E5M2, unscaled     0.2178 0.0000   5.95e-02
  E5M2, per-tensor   0.0000 0.0000   1.96e-02
  E4M3, unscaled     1.0000 0.0000   1.00e+00
  E4M3, per-tensor   0.0000 0.0000   1.22e-02

after 3000 FP32 steps: |dz| from 1.2e-15 to 2.4e-04
  dz format, scale     zero    inf   gradient rel. error
  fp16, S = 2^0      0.4690 0.0000   3.29e-02
  fp16, S = 2^4      0.2688 0.0000   1.44e-02
  fp16, S = 2^8      0.1187 0.0000   1.43e-02
  fp16, S = 2^12     0.0347 0.0000   1.44e-02
  fp16, S = 2^16     0.0066 0.0000   1.44e-02
  fp16, S = 2^20     0.0015 0.0000   1.44e-02
  fp16, S = 2^24     0.0002 0.0000   1.44e-02
  fp16, S = 2^28     0.0000 0.0000   1.44e-02
  fp16, S = 2^32     0.0000 0.1365   nan
  bf16, S = 1        0.0000 0.0000   1.09e-01
  E5M2, unscaled     0.8279 0.0000   4.50e+00
  E5M2, per-tensor   0.0017 0.0000   3.02e+00
  E4M3, unscaled     1.0000 0.0000   1.00e+00
  E4M3, per-tensor   0.1885 0.0000   1.14e+00

camera pixel sum: exact 132674.8857
  fp16 sequential    2048.0000  rel. error 9.85e-01
  fp16 pairwise            inf  rel. error inf
  fp32 sequential  132773.5000  rel. error 7.43e-04
  fp32 pairwise    132674.8750  rel. error 8.03e-08
weight gradient A^T dz (4,096 terms per entry):
  fp16 sequential rel. error 2.08e-01
  fp16 pairwise   rel. error 3.14e-02
  fp32 sequential rel. error 2.68e-05
  fp32 pairwise   rel. error 4.06e-06

The checks at the top confirm the rounding and the half-ulp rule. An increment of $2^{-12}$ to 1 vanishes in FP16, $2^{-11}$ is an exact tie and rounds to the even neighbour, which is 1, and only the larger increment moves the weight to the next FP16 number, $1 + 2^{-10}$. FP32 keeps the small increment.

Training. With an FP32 master copy, all three reduced-precision regimes end within about 2% of the FP32 run’s excess loss, with the same test accuracy. Without one, both 16-bit regimes freeze: at the last step every nonzero update was lost against its weight, and the excess loss is about 21 times FP32’s for FP16 and about 150 times for BF16, whose spacing is eight times coarser. Loss scaling made little difference here. The run without it ends slightly worse than the run with it, while the scaler, starting at $2^{16}$, never overflowed and doubled once after 2,000 steps. In our reading the reason is that the activation gradients that underflow are those of the best-classified images, which contribute least to the gradient, and that deep networks, whose activation gradients can shrink from layer to layer (as DALL·E’s did, under Variants), are where loss scaling matters. The frozen 16-bit regimes reach slightly higher test accuracy than FP32, though: on this data a lower training loss is not a better classifier, and, in our reading, stalling acts like early stopping. The point of the master copy is fidelity to the FP32 optimization, not accuracy for its own sake.

Underflow against the scale. After 10 steps no activation gradient underflows even at $S = 1$, since the smallest, $4.1 \times 10^{-7}$, lies above FP16’s smallest subnormal, $2^{-24} \approx 6 \times 10^{-8}$. After 3,000 steps, near the optimum, 47% of the nonzero activation gradients round to zero in FP16. Each factor of 16 in the scale recovers part of them, and by $2^{28}$ none is lost. Overflow comes at $2^{32}$, where the largest values pass 65,504; this is the trade-off that dynamic scaling negotiates. The error column shows an earlier limit: after 10 steps the stored weight gradient overflows at $2^{24}$ already, although no activation gradient does, because the weight gradient, a sum of many scaled terms, is larger. This is the overflow that GradScaler detects in the weight gradients. Once underflow is removed, the FP16 error after 3,000 steps levels off at about 1.4%, much larger than after 10 steps. In our reading, near the optimum the weight gradient is a small difference of large terms, so the relative rounding error of the inputs is amplified. BF16 underflows nothing but has a larger error, set by its 7-bit mantissa.

FP8. Without scaling, E5M2 loses 22% and then 83% of the activation gradients, and E4M3, whose smallest subnormal is $2^{-9}$, loses all of them. A per-tensor scale that maps the largest value to the format’s maximum recovers nearly all for E5M2 and, after 10 steps, for E4M3. Late in training, E4M3 still loses 19%: in our reading, because the activation gradients then span far more binades than its 18. Near the optimum, the weight gradient computed from FP8 activation gradients has errors larger than the gradient itself. The example puts one tensor of a nearly converged model through FP8, an extreme case; the recipes under Variants use E5M2 for gradients, scales per tensor or per block, and higher precision for sensitive layers.

Accumulation. Summing the camera image’s pixels one at a time in FP16 stops at exactly 2,048. From there the spacing of FP16 numbers is 2, and adding any pixel value of at most 1 rounds back to the running sum. The FP32 running sum drifts by 0.07%, and pairwise FP32 summation is accurate to $10^{-7}$. We expected pairwise FP16 summation to come far closer to the true sum than sequential FP16; that check failed, because the true sum, about 132,675, exceeds FP16’s maximum of 65,504, so the last levels of the tree overflow to infinity. An FP16 accumulator fails on range as well as on precision. In the weight gradient, whose sums stay within FP16’s range, pairwise FP16 does come closer: FP16 accumulation errs by 21% sequentially and 3% pairwise, against $3 \times 10^{-5}$ and $4 \times 10^{-6}$ in FP32. These are the reasons for accumulating into FP32 inside matrix multiplications and for running large reductions in FP32.

Variants

• Pure 16-bit training. NVIDIA’s guide notes that many networks match FP32 when all tensors are stored in FP16, while some need the FP32 master copy.
• BF16 mixed precision. BF16 compute with FP32 master weights and no loss scaling, the default of PyTorch’s CPU autocast and of TPU matrix units. Its range removes the scaler; in our reading, its 7-bit mantissa makes FP32 master weights and accumulation more necessary, not less.
• Per-layer scaling. DALL·E stored most parameters, Adam moments and activations in 16 bits (Ramesh et al., 2021, §2.4). One loss scale suffices when the smallest and largest activation gradients fit in FP16’s exponent range together. In DALL·E’s 12-billion-parameter transformer they did not: the activation-gradient norms fell from earlier to later residual blocks, and the later blocks underflowed. Ramesh et al. (2021, §2.4) kept a separate gradient scale for each residual block, as the DALL·E article describes.
• FP8 with per-tensor scaling. The OFP8 specification defines two encodings. E4M3 has a maximum of 448, no infinities and 18 binades of range; E5M2 has a maximum of 57,344, infinities and 32 binades. It requires conversions that round to nearest even, with saturating and non-saturating modes, and places the selection of scaling factors outside its scope. NVIDIA’s Transformer Engine primer uses E4M3 for activations and weights in the forward pass, which need more precision, and E5M2 for gradients, which need more range. Its authors note that FP8’s range suffices for any single tensor but not for all at once, so a single loss scale is infeasible and each tensor gets its own scale. Just-in-time scaling derives the scale from the current tensor’s maximum absolute value (amax), which the primer calls infeasible in practice because it needs several passes over the data. Delayed scaling uses the amaxes of previous iterations, kept as a history. Transformer Engine also offers a current-scaling recipe. The FP8 formats paper of Micikevicius et al. (2022), not read for this article, is summarized in the quantization article.
• Block-scaled formats. MXFP8 gives every block of 32 consecutive values its own power-of-two scale (E8M0), which lets all tensors use E4M3. NVFP4 stores E2M1 values with an E4M3 scale per block of 16 and a per-tensor FP32 scale. Transformer Engine’s NVFP4 training recipe adds stochastic rounding of gradients, $16 \times 16$ blocks for weights, random Hadamard transforms against outliers, and recommends higher precision for the last few layers (primer, as of November 2026).

Best Practices

• Use the framework’s AMP, as NVIDIA’s guide recommends, rather than hand-casting. In PyTorch it applies the operation lists, keeps the FP32 parameters as master weights, and caches casts of the weights.
• Keep the loss and its reductions in FP32. Compute losses from logits (binary_cross_entropy_with_logits, not a sigmoid followed by binary_cross_entropy), and leave softmax, normalization statistics and large sums to FP32, as the operation lists do.
• Unscale before touching gradients. Clip, inspect or penalize gradients only after unscale_, so that thresholds keep their FP32 meaning (NVIDIA guide; GradScaler docstring).
• Prefer dynamic loss scaling for FP16, and expect a few skipped steps at the start; GradScaler’s docstring expects skipping afterwards to be rare, once every few hundred or thousand iterations.
• Consider BF16 where the hardware supports it. In our reading it removes the loss scale and its failure modes; keep FP32 master weights, which its coarse spacing makes more necessary.
• Respect Tensor Core shapes. Multiples of 8 for batch sizes, widths and channel counts on FP16 hardware (NVIDIA guide, §4.1); multiples of 16 for FP8 linear layers in Transformer Engine.

Pitfalls

• Lost updates without a master copy. With 16-bit weights, updates below half an ulp vanish. The example’s 16-bit-only runs froze completely.
• Unscaled clipping and weight decay. Clipping scaled gradients clips at the wrong threshold.
• Overflow in reductions. A sum, a norm or an exponential in FP16 can overflow even when its inputs are small; softmax’s exp overflows above about 11 in float16.
• A scale that never settles. If the scale keeps falling, gradients overflow for a reason scaling cannot fix. GradScaler’s docstring warns that a scale below 1 with NaNs in the loss or gradients means something is wrong, and names bf16-pretrained models in FP16 as an example.
• Casting the model by hand under autocast. Calling .half() on the model removes the FP32 master weights; GradScaler then refuses to unscale the FP16 gradients.
• Measuring the wrong thing. Speed-ups depend on the share of time spent in Tensor Core operations; a model limited by data loading or small layers may see little gain.

Evidence and Impact

Results come from the hardware vendor and from users of the technique:

• Accuracy (Micikevicius, 2017). NVIDIA’s blog post reports ILSVRC12 validation top-1 accuracy in FP32 and mixed precision: AlexNet 56.77% and 56.93%, VGG-D 65.40% and 65.43%, GoogLeNet 68.33% and 68.43%, Inception v1 70.03% and 70.02%, ResNet-50 73.61% and 73.75%. None of these needed loss scaling. On PASCAL VOC 2007 test, Faster R-CNN reached 69.1% mAP in FP32, 68.6% in mixed precision without loss scaling and 69.7% with it. Multibox SSD reached 76.9% in FP32, did not train in mixed precision without loss scaling, and reached 77.1% with a scale of 8. The bigLSTM language model needed a scale of 128 and otherwise diverged.
• Accuracy and throughput (NVIDIA guide, Table 3). On DGX-1V and DGX-2V systems:
  • ResNet-50 v1.5 in MXNet reached 76.67% ImageNet top-1 in FP32 and 76.49% in mixed precision, at 2,957 against 10,263 images per second (3.47 times).
  • SSD with a ResNet-50 backbone in TensorFlow reached 0.268 and 0.269 COCO 2017 mAP at 1.32 times the throughput.
  • U-Net Industrial gained 1.10 times.
  • Tacotron 2 and WaveGlow gained 1.18 and 1.94 times.
• Speed in a training recipe (He et al., 2019). On a V100, switching from FP32 to FP16 sped training up 2 to 3 times (§3.2). Table 3 reports ResNet-50 at 4.4 minutes per epoch against 13.3. That comparison also changes the batch size, from 256 in FP32 to 1,024 in FP16. In their ablation (Table 4), adding FP16 changed ResNet-50 top-1 from 76.16% to 76.15% at batch 256 and from 76.03% to 76.21% at batch 1,024.
• Scale. In DALL·E, getting a model past one billion parameters to train in 16-bit precision without diverging was, in its authors’ words, the most challenging part of the project; per-resblock gradient scaling fixed the underflow (Ramesh et al., 2021, §2.4).

The ICLR paper by Micikevicius et al. (2018) was not read for this article; the accuracy results above are those NVIDIA’s blog post reports, which refers readers to the paper for full experimental details. Mixed precision combines with sharding and activation checkpointing, which distributed training describes, and was used with them in the largest runs that article reports.