Low-Precision GPU Programming

A value of 500 goes through a scaled FP8 round trip and comes back as 10. No exception fires, no infinity gets written, and no NaN flag lights up. An earlier batch's scale made the new spike too large to encode, so saturating conversion clamped it to the format ceiling. The consumer GPU kernel can accumulate the changed input in FP32 without recovering the discarded information.

Smaller types can save storage and expose faster matrix instructions, but each format transition has its own hardware rate. Halving element width halves the element payload, not necessarily total traffic: scales, padding, conversions, and the output still occupy bytes. TF32 is an immediate counterexample to a storage-saving claim because its operands remain stored in FP32.

Shrinking bits isn't as simple as flipping a dtype argument in PyTorch. Doing low-precision math safely demands four explicit contracts:

1. Encoding: which exact bit patterns represent normal numbers, subnormals, zeros, infinities, and NaNs?
2. Scaling architecture: which values share a scale factor, when is that scale measured, and where does it live in memory?
3. Arithmetic and accumulation: which precision enters the multiplier ALUs, which type accumulates dot products, and where does dequantization happen?
4. Kernel application binary interface (ABI): which layouts, byte alignments, architecture targets, and swizzle patterns does the hardware instruction require?

Our projection $Y = XW$ has two rows of 32 activations and one 32-element weight column. One row contains values as small as $0.001$; the other contains a spike of $500$. These are constructed fixtures, not measured activation distributions. Predict whether protecting the small row also improves the spike's reconstruction. The two-row experiment will let you check both answers.

Format taxonomy: range, spacing, and bitfields

A floating-point format name reveals how many exponent and fraction bits it provides. Exponent bits determine dynamic range: the distance between the smallest subnormal and the largest finite number. Fraction bits (the mantissa) determine precision: the spacing between representable numbers within any given binary octave.

The Open Compute Project (OCP) OFP8 specification standardizes these FP8 encodings.^{[1]Reference 1OCP 8-bit Floating Point Specification (OFP8), Revision 1.0https://www.opencompute.org/documents/ocp-8-bit-floating-point-specification-ofp8-revision-1-0-2023-12-01-pdf-1} A common hybrid training recipe uses E4M3 for forward operands and E5M2 for gradients, trading precision for range:^{[2]Reference 2Using FP8 and FP4 with Transformer Enginehttps://docs.nvidia.com/deeplearning/transformer-engine/examples/fp8_primer.html}

• E4M3's extra fraction bit: Spacing in $[1,2)$ is 0.125, versus E5M2's 0.25. Relative spacing varies within a binade; 12.5% isn't a constant relative rounding error. E4M3 reaches 448, has no infinities, and reserves 0x7F and 0xFF for NaNs.
• E5M2's wider range: Its smallest subnormal is $2^{-16}$ and largest finite magnitude is 57,344. This gives more room for a scaled gradient distribution, not immunity to underflow. An unscaled $10^{-7}$ still falls below its representable floor. Format choice and scaling need evidence from the actual tensors.^{[3]Reference 3FP8 Formats for Deep Learning.https://arxiv.org/abs/2209.05433}

Normalization isn't a guarantee that every activation lies in $[-20,20]$. It constrains a statistic, and learned affine parameters or later operations can change range. For example, before affine parameters or epsilon, LayerNorm applied to one nonzero value among 1,024 entries gives that entry a magnitude of $\sqrt{1023}\approx31.98$.

Beyond 8-bit floats sits the sub-8-bit frontier:

• OCP Microscaling (MX): The specified formats use 32-element blocks with an 8-bit E8M0 scale. Elements retain their own narrower encoding: MXFP8 uses E4M3 or E5M2, MXFP6 uses E3M2 (maximum 28) or E2M3 (7.5), and MXFP4 uses E2M1 (6).^{[4]Reference 4OCP Microscaling Formats (MX) Specification, Version 1.0https://www.opencompute.org/documents/ocp-microscaling-formats-mx-v1-0-spec-final-pdf}
• NVIDIA NVFP4: E2M1 elements use an E4M3 scale per 16-element block and an additional FP32 tensor scale. This differs from MXFP4's 32-element/E8M0 contract. Transformer Engine's training recipe can share one weight scale across a $16\times16$ block; its activation and gradient scales remain one-dimensional. Smaller weight storage doesn't imply the attention KV cache uses the same format.^{[5]Reference 5Transformer Engine NVFP4 Recipehttps://docs.nvidia.com/deeplearning/transformer-engine/features/low_precision_training/nvfp4/nvfp4.html}

Grid spacing increases at binade boundaries: E4M3's spacing is 0.125 in $[1,2)$ and 32 from 256 up to its finite ceiling. Actual error also depends on position within the grid and can be zero. Scaling manages range; a power-of-two scale doesn't add fraction bits or improve relative precision for values that stay normal.

Scale ownership and the GEMM dataflow

Let's derive the exact mathematics of how scale factors travel through a matrix multiplication. Consider a standard GEMM $Y = XW$, with activation matrix $X \in \mathbb{R}^{M \times K}$ and weight matrix $W \in \mathbb{R}^{K \times N}$.

Every element is represented by a positive scale factor and a quantized low-precision code:

Here $\tilde{X}_{ik}$ and $\tilde{W}_{kj}$ are low-precision FP8 or FP4 values, while $S_{X, i}$ and $S_{W, j}$ are higher-precision (FP32 or E8M0) scale factors. Substitute these into the dot product definition for output element $Y_{ij}$:

Because row scale $S_{X, i}$ and column scale $S_{W, j}$ don't depend on the contraction index $k$, both scales factor completely out of the inner sum:

These are real-arithmetic identities for the reconstructed operands, denoted by the hat. Input quantization has already changed the operation, and finite-precision execution needn't preserve a bitwise equality when operations are reordered. For a per-row/per-column FP8 kernel, work can divide into two phases:

1. Matrix multiply-accumulate: The inner sum consumes the low-precision operands with the chosen instruction and accumulator mode. Storage in a 32-bit accumulator doesn't prove every intermediate step has full FP32 precision.
2. Epilogue Dequantization: Dequantization scale $(S_{X, i} \cdot S_{W, j})$ is applied in the kernel's epilogue stage right before writing the final result back to global memory or piping it into a fused bias and activation function.

This raises a critical architecture question: who owns quantization and dequantization?

• A possible producer fusion: RMSNorm, an activation, or a preceding epilogue can also generate codes and scales, when its range reduction and output ABI permit that fusion. A separate conversion kernel is another valid choice.
• A possible consumer fusion: Applying scales inside GEMM avoids materializing a wider intermediate. Standalone dequantization remains useful for debugging, fallback operators, or an unsupported consumer. An FP32 intermediate contains four bytes per FP8 element, but total traffic and end-to-end performance depend on all reads, writes, scales, and kernels.

If scaling happens block-wise along the contraction axis $K$ (as in 32-element microscaling with block index $b$):

The block scale product stays inside the outer sum, with a partial final block zero-padded when required. This equation specifies meaning, not one mandatory instruction sequence. OCP leaves internal dot-product precision and operation order implementation-defined; supported block-scaled matrix instructions can consume scale operands directly.^{[4]Reference 4OCP Microscaling Formats (MX) Specification, Version 1.0https://www.opencompute.org/documents/ocp-microscaling-formats-mx-v1-0-spec-final-pdf[6]Reference 6PTX ISA: Fifth-Generation TensorCore Instructionshttps://docs.nvidia.com/cuda/parallel-thread-execution/#tensorcore-5th-generation-instructions}

Numerical hazards: subnormals, saturation, and rounding

Running low-precision code exposes three distinct numerical pitfalls. Understanding how they manifest keeps models from silently losing fidelity.

Subnormals and flush-to-zero (FTZ)

Normal floating-point numbers have an implicit leading 1 ($1.f \times 2^{e - \text{bias}}$). Subnormal numbers occur when the exponent bits are all zero; the leading bit drops to 0 ($0.f \times 2^{1 - \text{bias}}$). This allows values to shrink smoothly toward zero, albeit with progressively fewer significant bits.

In OCP E4M3, the smallest positive normal is $2^{-6} = 0.015625$. Subnormals extend down to $2^{-9} \approx 0.001953125$. Under standard round-to-nearest-even (RNE):

• An input of $0.75 \times 2^{-9}$ rounds up to the smallest subnormal $2^{-9}$.
• An input of $0.5 \times 2^{-9}$ lands on an exact tie and rounds down to signed zero (+0.0 or -0.0).

Flush-to-zero behavior belongs to the specific instruction and mode. As a concrete example, PyTorch documents input/output denormal flushing for MI200 FP16/BF16 matrix operations and alternate FP16 paths used during backward computation. Don't infer a universal one-cycle penalty or model failure from a dtype name. Test the actual kernel's subnormal behavior.^{[7]Reference 7Numerical Accuracyhttps://docs.pytorch.org/docs/stable/notes/numerical_accuracy.html}

Saturation cliffs

OFP8 conversion first rounds to the destination precision. If that rounded magnitude exceeds the finite range, non-saturating conversion produces an E4M3 NaN or E5M2 infinity. A raw value just above the largest finite value can instead round back to it; distinguish the pre-cast magnitude from the rounded overflow condition.

Saturating conversion limits overflow to $\pm448$ or $\pm57,344$. CUDA exposes explicit saturation modes, and its FP8 constructors use __NV_SATFINITE; framework casts need their own check. OFP8 requires RNE and both saturation modes to be supported, rather than one universal default.^{[1]Reference 1OCP 8-bit Floating Point Specification (OFP8), Revision 1.0https://www.opencompute.org/documents/ocp-8-bit-floating-point-specification-ofp8-revision-1-0-2023-12-01-pdf-1[8]Reference 8FP8 Conversion and Data Movementhttps://docs.nvidia.com/cuda/cuda-math-api/cuda_math_api/group__CUDA__MATH__FP8__MISC.html}

Saturating E4M3 conversion maps both $500$ and $50,000$ to $448$, losing their relative magnitude. The values remain finite, so NaN-only checks miss the change. Model quality needs a separate evaluation. Backpropagation depends on the implemented estimator: the exact clipping function's flat region has zero derivative, but training often uses a surrogate or straight-through rule.

Rounding: RNE versus stochastic rounding

OFP8 requires support for RNE and permits additional rounding modes. Consider repeated updates to a stored low-precision weight, with no wider master value or residual to retain the lost increment.

If $W+\Delta W$ stays within W's nearest-value interval, rounding the updated weight returns W. The update isn't necessarily zero as a separately represented number. Neighboring intervals can have different spacing at binade boundaries:

Repeating the same absorbed increment can leave that stored weight unchanged. This is a conditional failure mode, not a claim that every mixed-precision optimizer freezes its weights.

Stochastic rounding breaks this trap by rounding up or down probabilistically based on the distance to the adjacent representable points:

Here l and u are neighboring representable floating-point values, not the integer floor and ceiling. For finite $l<x<u$ with these probabilities, $\mathbb{E}[\operatorname{SR}(x)]=x$; an exactly representable value stays unchanged. Clipping outside the finite range breaks this unbiasedness. Stochastic rounding can help preserve expected increments, but it doesn't guarantee convergence or make wider optimizer state unnecessary. Transformer Engine's NVFP4 recipe applies it to gradients, not every tensor.^{[5]Reference 5Transformer Engine NVFP4 Recipehttps://docs.nvidia.com/deeplearning/transformer-engine/features/low_precision_training/nvfp4/nvfp4.html}

from fractions import Fraction

l, u, x = Fraction(1), Fraction(9, 8), Fraction(101, 100)
p_upper = (x - l) / (u - l)
nearest = l if x - l < u - x else u
assert nearest == 1 and p_upper == Fraction(2, 25)
assert p_upper * u + (1 - p_upper) * l == x
clipped = min(Fraction(500), Fraction(448))
assert clipped != 500
print(f"neighbors={float(l):.3f},{float(u):.3f}; RNE={float(nearest):.3f}")
print(f"upper probability={float(p_upper):.3f}; expected rounded value={float(x):.3f}")
print("clipped 500 -> 448: expectation can't equal the original")

neighbors=1.000,1.125; RNE=1.000
upper probability=0.080; expected rounded value=1.010
clipped 500 -> 448: expectation can't equal the original

Accumulator type isn't the whole arithmetic contract

Prefer an accumulator appropriate for reduction length and input range, and verify the supported instruction. FP16 accumulation is legal for some Tensor Core operations; full FP32 isn't a universal hardware-enforced rule.^{[9]Reference 9Parallel Thread Execution ISAhttps://docs.nvidia.com/cuda/parallel-thread-execution/}

The preceding GEMM lesson shows absorption of a positive 0.25 update at 1024 in FP16 and its preservation in FP32. It also shows why negative updates near that binade boundary need different spacing. K=4096 alone doesn't establish a partial sum's magnitude or when absorption happens.

Even a declared FP32 accumulator can involve reduced internal precision or intermediate reductions. Inspect the operation's guarantees and framework settings, including FP16/BF16 reduced-precision reduction controls. The ABI must distinguish input, product/internal accumulation, accumulator storage, and output types.^{[7]Reference 7Numerical Accuracyhttps://docs.pytorch.org/docs/stable/notes/numerical_accuracy.html}

Hands-on: auditing FP8 conversion boundaries

Test subnormal rounding, an explicit finite-input clamp, and byte-level decoding on CPU. This isn't a test of GPU matrix arithmetic or CUDA's non-saturating conversion mode.

The plain cast and explicit clamp agree for this fixture in the recorded PyTorch build. That observation contradicts a blanket claim that a plain CPU cast must yield NaN at 500. Rerun it on your installed build rather than inferring the overflow mode from the dtype tag.

import torch

torch.set_num_threads(1)

print("PyTorch:", torch.__version__, "device: cpu")
d = 2.0 ** -9
values = torch.tensor([0.5 * d, 0.75 * d, -0.5 * d, 1.0625, 1.1875, 500.0])

# Observe this build's plain CPU cast; don't assume CUDA __NV_NOSAT semantics.
raw = values.to(torch.float8_e4m3fn).float()

# Saturating conversion: clamps finite values to the E4M3 ceiling of 448
saturated = values.clamp(-448.0, 448.0).to(torch.float8_e4m3fn).float()

print("plain cast:", raw.tolist())
print("saturating finite inputs:", saturated.tolist())
print("negative half-step keeps sign:", bool(torch.signbit(raw[2])))

assert raw[:5].tolist() == [0.0, d, -0.0, 1.0, 1.25]
assert saturated[-1] == 448.0
assert torch.isfinite(saturated).all()

# The exact same byte pattern means different values across specs
codes = torch.tensor([0x38, 0x80], dtype=torch.uint8)
print("OCP E4M3 bytes:", codes.view(torch.float8_e4m3fn).float().tolist())
print("FNUZ E4M3 bytes:", codes.view(torch.float8_e4m3fnuz).float().tolist())

PyTorch: 2.14.1+cu130 device: cpu
plain cast: [0.0, 0.001953125, -0.0, 1.0, 1.25, 448.0]
saturating finite inputs: [0.0, 0.001953125, -0.0, 1.0, 1.25, 448.0]
negative half-step keeps sign: True
OCP E4M3 bytes: [1.0, -0.0]
FNUZ E4M3 bytes: [0.5, nan]

Look at the byte decoding at the end: byte 0x38 decodes to 1.0 in OCP E4M3, but becomes 0.5 in AMD's FNUZ E4M3 because their exponent biases differ (7 versus 8). Byte 0x80 represents negative zero in OCP E4M3, but represents NaN in FNUZ because FNUZ drops signed zeros. Feeding an OCP-encoded buffer into a kernel expecting FNUZ silently alters standard numeric values throughout your model.^{[10]Reference 10Tensor Attributeshttps://docs.pytorch.org/docs/stable/tensor_attributes.html[11]Reference 11Data Types and Precision Supporthttps://rocm.docs.amd.com/en/docs-10.0.0/reference/precision-support.html}

Worked failure: stale delayed scale

Suppose delayed scaling calculates an E4M3 scale factor from a previous iteration's maximum absolute value ($\text{amax} = 10.0$):

In the current iteration, an activation outlier spikes to $500.0$. The scaled magnitude is:

Under saturating conversion, $22,400.0$ clamps directly to the format maximum of $448$. During consumer dequantization, reconstruction returns:

The reconstructed value is finite, positive, and shows no NaN errors, but its relative error is $98\%$. The FP32 accumulator faithfully adds the corrupted product into the running sum, masking the numerical collapse from monitoring systems.

Scaling architectures: granularity and delayed scaling

The granularity of your scale factor sets which values compete for the same dynamic range exponent.

One letter, two reciprocal conventions

This chapter uses S for reconstruction: $\hat{x}=S q$ and $q=Q(x/S)$. Many APIs call the reciprocal $\alpha=1/S$ the “scale”: $q=Q(\alpha x)$, followed by $\hat{x}=q/\alpha$. Read the formula before interpreting a field named scale or scale_inv.

For finite amax > 0 and an ideal unrounded E4M3 reconstruction scale, $S=\text{amax}/448$, whereas $\alpha=448/\text{amax}$. A zero-only group needs an explicit positive-scale policy, such as S=1; division by amax isn't valid there.

Current scaling: where the reduction happens

Current scaling derives a scale from the tensor being encoded:

Copying a device amax to the host for a dependent decision adds transfer and synchronization costs. Their duration depends on the system and path, not a universal 10–50 microsecond penalty per layer. A device reduction and dependent cast can stay ordered on the device.

Transformer Engine 2.20.2's current-scaling documentation describes two tensor reads during quantization: one for amax and one for scaling/casting. GEMM subsequently consumes the encoded tensor. Producer fusion can change that cost, but only if the reduction and layout dependencies permit it.^{[12]Reference 12Transformer Engine FP8 Current Scalinghttps://docs.nvidia.com/deeplearning/transformer-engine/features/low_precision_training/fp8_current_scaling/fp8_current_scaling.html}

Delayed scaling in Transformer Engine

Delayed scaling uses prior amax observations to choose the current scale. With a maximum-over-history rule, no margin, and positive history amax:^{[13]Reference 13Transformer Engine FP8 Delayed Scalinghttps://docs.nvidia.com/deeplearning/transformer-engine/features/low_precision_training/fp8_delayed_scaling/fp8_delayed_scaling.html}

The quantization pass also records the current amax for future scale updates. Compared with its current-scaling recipe, Transformer Engine describes one tensor read rather than two. This doesn't eliminate all kernels or reductions: distributed recipes may still reduce amax across ranks, and history/scale updates have their own ordering.

The risk is scale lag. If an activation spike hits after a distribution shift, the historical scale is too small, causing massive saturation clipping. If the historical maximum was an unusual outlier, the scale remains too large, pushing typical activations down toward the subnormal floor.

Online scaling isn't automatically fused into TMA

DeepSeek-V3's report measures amax online for $1\times128$ activation groups and $128\times128$ weight blocks. It reports promoting partial FP8 GEMM results to CUDA-core FP32 accumulation every 128 contraction elements to address its observed H800 internal-precision limitation.^{[14]Reference 14DeepSeek-V3 Technical Report.https://arxiv.org/abs/2412.19437}

The report explicitly describes separate BF16 reads, quantization, FP8 writes, and subsequent MMA reads. Fusing FP8 casting with TMA is its proposed future optimization, not the implemented mechanism. The pinned DeepGEMM README likewise tells callers to handle FP8 casting and transposition separately or fuse them into preceding kernels. Online scales reduce history lag, but don't prove no clipping, no conversion pass, or no synchronization.^{[14]Reference 14DeepSeek-V3 Technical Report.https://arxiv.org/abs/2412.19437[15]Reference 15DeepGEMM Source Snapshot, July 15 2026https://github.com/deepseek-ai/DeepGEMM/tree/559d79fb6994a58b8a15b4b93bf13ccc16edf247}

fp8_limit, old_amax, spike = 448.0, 10.0, 500.0
S, alpha = old_amax / fp8_limit, fp8_limit / old_amax
assert old_amax / S == old_amax * alpha == fp8_limit
reconstructed = min(spike / S, fp8_limit) * S
error_percent = 100 * abs(reconstructed - spike) / spike
partials, scale_products = [2.0, 3.0], [5.0, 8.0]
scaled_sum = sum(p * s for p, s in zip(partials, scale_products))
wrong_sum = scale_products[0] * sum(partials)
assert (reconstructed, error_percent, scaled_sum, wrong_sum) == (10.0, 98.0, 34.0, 25.0)
print(f"S={S:.8f}, alpha={alpha:.1f}; peak code={old_amax / S:.1f}")
print(f"stale reconstruction: {spike:.1f} -> {reconstructed:.1f} ({error_percent:.1f}% error)")
print(f"block-scale sum={scaled_sum:.1f}, incorrectly pulled-out scale={wrong_sum:.1f}")

S=0.02232143, alpha=44.8; peak code=448.0
stale reconstruction: 500.0 -> 10.0 (98.0% error)
block-scale sum=34.0, incorrectly pulled-out scale=25.0

Microscaling: the two-row experiment

An E8M0 scale byte c from 0 through 254 represents $2^{c-127}$; 255 is NaN. It has no zero, sign, or infinity encoding. A power-of-two scale can be applied through exponent adjustment in real arithmetic, but that doesn't establish a chip's exact datapath or instruction cost.^{[4]Reference 4OCP Microscaling Formats (MX) Specification, Version 1.0https://www.opencompute.org/documents/ocp-microscaling-formats-mx-v1-0-spec-final-pdf}

A 32-value block holds 256 bits of element data plus 8 bits of shared scale, giving an effective payload of:

import math

def decode_e8m0(code):
    if type(code) is not int or not 0 <= code <= 255:
        raise ValueError("an unsigned scale byte is required")
    return math.nan if code == 255 else 2.0 ** (code - 127)

scale_codes = [-8 + 127, 1 + 127]
assert [decode_e8m0(c) for c in scale_codes] == [2.0 ** -8, 2.0]
assert decode_e8m0(0) > 0 and math.isnan(decode_e8m0(255))
assert (64 + 4, 64 + 2) == (68, 66)
print(f"E8M0 scale bytes: {[hex(c) for c in scale_codes]}; 0x00 is positive; 0xff is NaN")
print("packed X payload: tensor=68 B, MX=66 B; scale overhead=6.25%/3.125%")

E8M0 scale bytes: ['0x77', '0x80']; 0x00 is positive; 0xff is NaN
packed X payload: tensor=68 B, MX=66 B; scale overhead=6.25%/3.125%

These packed-payload counts exclude W, container padding, duplicated orientations, history, and runtime allocations. A tensor's four scale bytes aren't zero overhead. NVFP4 similarly has block and tensor scales; four-bit elements don't promise a 50% reduction in total storage.

Test the two-row matrix on CPU. We model conservative power-of-two reconstruction scales by rounding amax/448 upward, as described in Transformer Engine's current MXFP8 guide. OCP also describes another scale-selection algorithm; identical format names don't guarantee identical quantization recipes.^{[16]Reference 16MXFP8https://docs.nvidia.com/deeplearning/transformer-engine/features/low_precision_training/mxfp8/mxfp8.html[4]Reference 4OCP Microscaling Formats (MX) Specification, Version 1.0https://www.opencompute.org/documents/ocp-microscaling-formats-mx-v1-0-spec-final-pdf}

The script stores scales as FP32 tensors, reconstructs values to FP32, and uses CPU FP32 matrix multiplication. It doesn't pack E8M0 scale bytes, implement a GPU scale swizzle, or execute MXFP8 hardware. In particular, two rows don't satisfy Transformer Engine's documented multiple-of-32 constraint for the nonfinal dimensions.

Predict what happens to $0.001$ when it shares one global scale with $500.0$, versus when each row gets its own 32-element block scale:

import math
import torch

torch.set_num_threads(1)
torch.set_printoptions(precision=6, sci_mode=False)

pattern_small = [0.001, 0.03, 0.117, 0.5, 1.0, -0.25, 0.0625, -0.015625]
pattern_large = [4.0, 16.0, 64.0, 500.0, -2.0, -8.0, -32.0, -250.0]
weight_pattern = [0.5, -0.25, 0.125, -0.0625, 0.03125, -0.5, 0.25, -0.125]

x = torch.tensor([pattern_small * 4, pattern_large * 4], dtype=torch.float32)
w = torch.tensor(weight_pattern * 4, dtype=torch.float32).reshape(32, 1)
fp8_max = torch.finfo(torch.float8_e4m3fn).max

def roundtrip(values: torch.Tensor, scale: torch.Tensor) -> torch.Tensor:
    if not torch.isfinite(values).all() or not torch.isfinite(scale).all() or (scale <= 0).any():
        raise ValueError("finite values and positive finite scales required")
    scaled = values / scale
    return scaled.clamp(-fp8_max, fp8_max).to(torch.float8_e4m3fn).float() * scale

# 1. Tensor-wide scaling: single scale derived from global maximum
tensor_scale_x = x.abs().max() / fp8_max
tensor_scale_w = w.abs().max() / fp8_max
x_tensor = roundtrip(x, tensor_scale_x)
w_tensor = roundtrip(w, tensor_scale_w)

# 2. Block-level scaling: conservative power-of-two scale per 32-value row
def conservative_power_of_two_scale(block: torch.Tensor) -> float:
    ratio = block.abs().max().item() / fp8_max
    if ratio == 0:
        return 1.0
    exponent = max(-127, math.ceil(math.log2(ratio)))
    return 2.0 ** exponent

mx_scales_x = torch.tensor(
    [conservative_power_of_two_scale(row) for row in x], dtype=torch.float32
).reshape(2, 1)
mx_scale_w = torch.tensor(conservative_power_of_two_scale(w), dtype=torch.float32)

x_mx = roundtrip(x, mx_scales_x)
w_mx = roundtrip(w, mx_scale_w)

reference_y = (x @ w).squeeze()
tensor_y = (x_tensor @ w_tensor).squeeze()
mx_y = (x_mx @ w_mx).squeeze()

print(f"tensor_scale_x={tensor_scale_x.item():.6f}")
print(f"mxfp8_scales_x={mx_scales_x.squeeze().tolist()}")
print(f"x[0,0]: ref={x[0, 0].item():.8f} tensor_fp8={x_tensor[0, 0].item():.8f} mxfp8={x_mx[0, 0].item():.8f}")
print(f"x[1,3]: ref={x[1, 3].item():.8f} tensor_fp8={x_tensor[1, 3].item():.8f} mxfp8={x_mx[1, 3].item():.8f}")
print(f"reference_y=[{reference_y[0]:.6f}, {reference_y[1]:.6f}]")
print(f"tensor_fp8_y=[{tensor_y[0]:.6f}, {tensor_y[1]:.6f}]")
print(f"mxfp8_y=[{mx_y[0]:.6f}, {mx_y[1]:.6f}]")
print(f"tensor_abs_error={(tensor_y - reference_y).abs().tolist()}")
print(f"mxfp8_abs_error={(mx_y - reference_y).abs().tolist()}")

tensor_scale_x=1.116071
mxfp8_scales_x=[0.00390625, 2.0]
x[0,0]: ref=0.00100000 tensor_fp8=0.00000000 mxfp8=0.00097656
x[1,3]: ref=500.00000000 tensor_fp8=500.00000000 mxfp8=512.00000000
reference_y=[0.600812, 7.750000]
tensor_fp8_y=[0.583104, 7.568359]
mxfp8_y=[0.601562, 7.750000]
tensor_abs_error=[0.017708778381347656, 0.181640625]
mxfp8_abs_error=[0.0007500052452087402, 0.0]

Follow the numbers closely:

1. Tensor-wide scale ($S = 1.116071$): The smallest nonzero reconstructed E4M3 step is about 0.00218. Our 0.001 falls below half that step and rounds to zero. The row's dot-product error is about 3% in this fixture, including quantization of both X and W; it isn't entirely attributable to this one small activation.
2. Under block-level scaling ($S_0 = 2^{-8} = 0.00390625$): The first block's scale shifts the dynamic range downward. Its smallest subnormal step becomes $2^{-8} \times 2^{-9} = 2^{-17} \approx 7.63\times10^{-6}$. The small value $0.001$ survives accurately as $0.00097656$ (just $2.3\%$ relative error).
3. Large row scale ($S_1=2$): The spike reconstructs as 512, a 2.4% error, although the row's dot product happens to equal its reference. Local scaling protects the other row; it doesn't make every value or result more accurate. Try changing W to observe how cancellation affects that dot product.

The zero projection error in the large row isn't evidence that its operands are accurate. Keep the encoded activations unchanged and replace W with a column that selects only the spike:

selector = torch.zeros_like(w)
selector[3, 0] = 1.0
selected_reference = (x @ selector).squeeze()
selected_mx = (x_mx @ selector).squeeze()

assert mx_y[1] == reference_y[1] == 7.75
assert selected_reference.tolist() == [0.5, 500.0]
assert selected_mx.tolist() == [0.5, 512.0]
print(f"original projection: reference={reference_y[1]:.2f}, MX={mx_y[1]:.2f}")
print(f"spike-only projection: reference={selected_reference[1]:.1f}, MX={selected_mx[1]:.1f}")
print(f"spike-only absolute error={abs(selected_mx[1] - selected_reference[1]):.1f}")

original projection: reference=7.75, MX=7.75
spike-only projection: reference=500.0, MX=512.0
spike-only absolute error=12.0

Cancellation in the original dot product hid the reconstruction error. A different consumer exposes it without any new quantization. Validate reconstructed operands and representative downstream computations; one matching output doesn't establish a general accuracy guarantee.

For a fixed logical matrix, changing from rowwise 1-by-32 groups to columnwise 32-by-1 groups changes scale membership. Transposing element bytes doesn't generally create the second quantization. Transformer Engine derives both orientations from the original wider data to avoid another quantization loss. Special equal-scale cases can agree, and transposing a representation while carrying its original group meaning isn't the same as regrouping it for a different contraction.^{[16]Reference 16MXFP8https://docs.nvidia.com/deeplearning/transformer-engine/features/low_precision_training/mxfp8/mxfp8.html}

Hardware support and kernel ABI contracts

Hardware claims must always be tied to specific GPU microarchitectures, compute capabilities, and driver toolkits.

This table reflects documentation inspected on October 10, 2026, including ROCm's versioned 10.0.0 architecture table. In PyTorch 2.14, use torch.backends.cuda.matmul.fp32_precision = "tf32" or "ieee" to permit TF32 or request IEEE precision. Don't mix these controls with legacy allow_tf32 settings.^{[7]Reference 7Numerical Accuracyhttps://docs.pytorch.org/docs/stable/notes/numerical_accuracy.html[16]Reference 16MXFP8https://docs.nvidia.com/deeplearning/transformer-engine/features/low_precision_training/mxfp8/mxfp8.html[5]Reference 5Transformer Engine NVFP4 Recipehttps://docs.nvidia.com/deeplearning/transformer-engine/features/low_precision_training/nvfp4/nvfp4.html[11]Reference 11Data Types and Precision Supporthttps://rocm.docs.amd.com/en/docs-10.0.0/reference/precision-support.html}

Before launching a low-precision kernel in production, write out its complete numeric contract:

kernel_contract = {
  operation: GEMM_Y_equals_XW,
  operand_x: { dtype: OCP_E4M3, layout: row_major, scale_granularity: per_token_1x32, scale_dtype: E8M0 },
  operand_w: { dtype: OCP_E4M3, layout: col_major_packed, scale_granularity: per_channel_32x1, scale_dtype: E8M0 },
  accumulator: FP32,
  scale_application: { location: inside_block_scaled_contraction, varies_along_K: true },
  epilogue: { bias: FP32, activation: SiLU, output_dtype: BF16 },
  rounding_mode: round_nearest_even,
  overflow_policy: saturate,
  alignment_and_scale_layout: record_requirements_of_the_selected_kernel,
  internal_accumulation: record_instruction_guarantees_and_framework_reduction_modes
}

Treat framework, compiler, and kernel boundaries as strict serialization interfaces. Passing dtype tags, scales, orientation, padding, and accumulator modes together prevents silent numerical drift.

Troubleshooting numerical failures

When low-precision outputs diverge from an FP32 reference run, use this diagnostic roadmap:

A reliable incident trace must capture input histograms, amax per scale group, selected scale factors, clip counts, zero counts, accumulator modes, and the profiler-confirmed kernel symbol. Checking only the final output relative error hides where the contract failed.

Try three changes to the two-row fixture: make the small row all zeros, move 500 into its first 32-value group, and replace W with a column selecting only the spike. Predict which values lose precision and whether the projection error reveals that loss. Keep the code/scale bytes and reconstruction formula together; an accurate dot product can hide inaccurate operands through cancellation.