gemm_4bit (via matmul_4bit / MatMul4Bit) and dequantize_4bit + F.linear differ by up to 1 ulp for M <= 256 on sm_120 (RTX 5070 Laptop), fp16 and bf16 - intended, and can the path be pinned? #2121

Summary

For the same packed NF4 weight, bnb.matmul_4bit(x, packed.t(), quant_state=state) and F.linear(x, dequantize_4bit(packed, state)) do not return bit-identical results on an RTX 5070 Laptop GPU (sm_120). The difference only appears for the row counts where gemm_4bit selects its custom kernel (M = 1 up to 256 for a 4096 x 14336 weight, up to 64 for 4096 x 4096 and 14336 x 4096 in my sweep), and it is exactly zero where gemm_4bit falls back to dequantize + GEMM.

This is probably a rounding-order difference and not a defect: the largest difference I see is exactly one ulp of the output dtype at the largest output magnitude, and both results are within about one ulp (at most 1.5 in all runs I did) of an fp32 product of the same dequantized weights. I am reporting it because the choice between the two paths depends on M, so the same weights give slightly different numbers when a caller changes batch size or sequence length, and because exact comparisons between a Linear4bit model and one that dequantizes explicitly cannot succeed. Two questions at the end.

Not a crash and not an accuracy loss; if this is expected behaviour, a sentence in the documentation would already help.

Environment

• GPU: NVIDIA GeForce RTX 5070 Laptop GPU (sm_120), 8123 MiB
• Driver: 616.92
• torch: 2.14.0+cu130 (CUDA 13.0)
• bitsandbytes: 0.50.2 (PyPI wheel, libbitsandbytes_cuda130.dll)
• Python 3.12.10, Windows-11-10.0.26200-SP0

I could not test other architectures (no sm_80/86/89/90 card here), so I do not know whether this is specific to sm_120.

Reproduction

No model download and no dependency beyond torch and bitsandbytes. One CUDA GPU.

"""NF4: bnb.matmul_4bit (gemm_4bit / MatMul4Bit) vs dequantize_4bit + F.linear."""

import math

import bitsandbytes as bnb
import torch
import torch.nn.functional as F
from bitsandbytes.functional import dequantize_4bit, quantize_4bit

print("torch", torch.__version__, "cuda", torch.version.cuda, torch.cuda.get_device_name(0),
      "sm_%d%d" % torch.cuda.get_device_capability(0))
print("bitsandbytes", bnb.__version__, bnb.__file__)


def max_abs(a, b):
    return (a.float() - b.float()).abs().max().item()


K, N = 4096, 14336  # Llama-3-8B MLP up/gate projection
torch.manual_seed(0)
for dtype, mantissa_bits in ((torch.float16, 10), (torch.bfloat16, 7)):
    weight = torch.randn(N, K, device="cuda", dtype=dtype)
    packed, state = quantize_4bit(weight, blocksize=64, quant_type="nf4", compress_statistics=True)
    dense = dequantize_4bit(packed, state).to(dtype)
    print(f"\n{dtype}  K={K} N={N}   (max abs differences; ulp = spacing of {dtype} at max|output|)")
    print("    M  |fused-dequant|  |autograd-fused|  |fused-fp32|  |dequant-fp32|    ulp")
    for rows in (1, 4, 12, 32, 33, 64, 256, 512, 2048):
        x = torch.randn(rows, K, device="cuda", dtype=dtype)
        with torch.no_grad():
            fused = bnb.matmul_4bit(x, packed.t(), quant_state=state)  # torch.ops.bitsandbytes.gemm_4bit
            dequant = F.linear(x, dense)  # dequantize, then cuBLAS GEMM
            exact = F.linear(x.float(), dense.float())  # fp32 reference (TF32 is off by default)
        x_grad = x.clone().requires_grad_(True)
        autograd = bnb.matmul_4bit(x_grad, packed.t(), quant_state=state)  # MatMul4Bit.apply
        ulp = 2.0 ** (math.floor(math.log2(exact.abs().max().item())) - mantissa_bits)
        print(f"{rows:>5}  {max_abs(fused, dequant):>14.4g}  {max_abs(autograd, fused):>15.4g}  "
              f"{max_abs(fused, exact):>11.4g}  {max_abs(dequant, exact):>13.4g}   {ulp:>6g}")

Output (bitsandbytes 0.50.2, installed from PyPI)

torch 2.14.0+cu130 cuda 13.0 NVIDIA GeForce RTX 5070 Laptop GPU sm_120
bitsandbytes 0.50.2 <site-packages>\bitsandbytes\__init__.py

torch.float16  K=4096 N=14336   (max abs differences; ulp = spacing of torch.float16 at max|output|)
    M  |fused-dequant|  |autograd-fused|  |fused-fp32|  |dequant-fp32|    ulp
    1           0.125                0       0.1389        0.06543     0.25
    4          0.0625                0      0.09744        0.09744     0.25
   12          0.0625                0       0.1237         0.1237     0.25
   32          0.0625                0       0.1216         0.1216     0.25
   33           0.125                0       0.1199         0.1199     0.25
   64            0.25                0       0.1237         0.2322     0.25
  256           0.125                0       0.1257         0.1257     0.25
  512               0                0       0.1247         0.1247     0.25
 2048               0                0        0.126          0.126     0.25

torch.bfloat16  K=4096 N=14336   (max abs differences; ulp = spacing of torch.bfloat16 at max|output|)
    M  |fused-dequant|  |autograd-fused|  |fused-fp32|  |dequant-fp32|    ulp
    1               1                0        1.132         0.8694        2
    4               0                0       0.5573         0.5573        2
   12               0                0       0.9604         0.9604        2
   32               0                0       0.9818         0.9818        2
   33               0                0        1.001          1.001        2
   64               2                0       0.9927          2.361        2
  256               0                0        0.994          0.994        2
  512               0                0       0.9972         0.9972        2
 2048               0                0       0.9994         0.9994        2

Same script against the current main wheel (continuous-release_main, bitsandbytes-1.33.7.preview-py3-none-win_amd64.whl, whose metadata says 0.50.3.dev0, built from a616886 on 2026-10-27; main is three commits past 0.50.2 and none of them touch the CUDA 4-bit kernels): every row is identical to the table above (checked by diffing the two outputs).

I also ran a larger sweep (3 shapes: 4096x4096, 4096x14336, 14336x4096; fp16 and bf16; M = 1, 4, 8, 12, 32, 33, 64, 256, 512, 1536, 1537, 2048; weight std 1.0 and 0.02). Summary: in every (shape, dtype, weight scale) group the largest fused - dequant difference is exactly 1 ulp of the output dtype at the largest output of the row (0.03 to 4 absolute for unit-variance weights, up to 0.06 for std 0.02); it is non-zero only on rows where _gemm_4bit_use_custom_cuda is True and exactly zero on every other row; autograd - fused is zero in all 144 rows; both paths are within 1.5 ulp of the fp32 reference in all of them; and 0.50.3.dev0 (main) reproduces all 144 records exactly.

Observations

1. matmul_4bit with a gradient (MatMul4Bit.apply) and without (gemm_4bit directly) are bit-identical in every row (column autograd-fused is 0). So the difference is between the custom kernel and the dequantize + GEMM fallback inside gemm_4bit, not an artefact of the autograd Function.
2. fused != dequant only on rows where bitsandbytes.backends.cuda.ops._gemm_4bit_use_custom_cuda(...) returns True. On every row where it returns False the two are equal.
3. Each path is deterministic: repeating the same call gives identical bits.
4. Neither path is consistently closer to the fp32 reference: the fused kernel is closer for some shapes and the dequantize path for others.

Expected vs actual

Expected (what I assumed): matmul_4bit returns the same bits as dequantizing and multiplying, or the documentation says that the last bit may differ.
Actual: they differ by one ulp of the output dtype at M <= 256 (shape dependent), identical above.

Why it matters in practice (context)

I maintain a LoRA fine-tuning tool (Soup) that streams NF4 base layers through the GPU, and its test suite compares a streamed model against a resident Linear4bit model with exact equality. On sm_120 that comparison failed with logit differences of 4.9e-04 (fp16) and 3.9e-03 (bf16) because the two sides were running different bitsandbytes kernels, which is how I found this. Related but separate: #2034 (MatMul4Bit keeps B and quant_state on ctx outside save_for_backward, which breaks gradient checkpointing with recycled weight buffers; still open). I also saw #2098 (open PR, sm_120 single-row speed-up of the same kernel); if it lands, the M = 1 rows will change.

Questions

1. Is a one-ulp difference between the custom gemm_4bit kernel and the dequantize + F.linear fallback expected? If so, could the docstring of matmul_4bit / gemm_4bit say so?
2. Is there, or could there be, a supported way to pin the path (always the custom kernel, or always the fallback) so that results do not change with M? Today the only way I found is to bypass matmul_4bit and call dequantize_4bit + F.linear myself, or to call the private _gemm_4bit_use_custom_cuda to mirror the choice.

I am happy to test patches on this card.