diff --git a/tests/pytorch/references/mxfp4_qat_reference.py b/tests/pytorch/references/mxfp4_qat_reference.py
new file mode 100644
index 0000000000..d6e41e2ad5
--- /dev/null
+++ b/tests/pytorch/references/mxfp4_qat_reference.py
@@ -0,0 +1,49 @@
+# Copyright (c) 2022-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+#
+# See LICENSE for license information.
+
+"""Pure-PyTorch reference for MXFP4 weight fake-quantization.
+
+Bit-identical specification of transformer_engine.pytorch.mxfp4_qat
+(tex.mxfp4_fake_quantize): per 1x32 block, scale = 2^clamp(ceil(log2(amax/6)),
+-126, 125) derived from the amax bit pattern, RTNE onto the E2M1 grid,
+saturate at 6, multiply back. Non-finite values NaN-poison their whole block;
+the sign of -0 is preserved; the scale floor/cap keep every grid value exactly
+representable in bf16/fp32.
+"""
+import torch
+
+_E2M1_MAX = 6.0
+_MXFP4_BLOCK = 32
+
+
+def round_to_e2m1_grid(y: torch.Tensor) -> torch.Tensor:
+    """RTNE onto the E2M1 magnitude grid {0, .5, 1, 1.5, 2, 3, 4, 6}; input in [0, 6]."""
+    fine = torch.round(y * 2.0) * 0.5
+    mid = torch.round(y)
+    coarse = torch.round(y * 0.5) * 2.0
+    return torch.where(y <= 2.0, fine, torch.where(y <= 4.0, mid, coarse))
+
+
+def mxfp4_fake_quantize_reference(weight: torch.Tensor) -> torch.Tensor:
+    """Numerical oracle for the mxfp4_fake_quantize kernel (composite torch ops)."""
+    rows, cols = weight.shape
+    w32 = weight.contiguous().to(torch.float32).view(rows, cols // _MXFP4_BLOCK, _MXFP4_BLOCK)
+    amax = w32.abs().amax(dim=-1, keepdim=True)
+    nonfinite = ~torch.isfinite(amax)
+
+    bits = amax.view(torch.int32)
+    exp_field = bits >> 23
+    mantissa = bits & 0x7FFFFF
+    exp = exp_field - 129 + (mantissa > 0x400000).to(torch.int32)
+    exp = torch.where(exp_field > 0, exp, torch.full_like(exp, -126))
+    exp = exp.clamp(min=-126, max=125)
+    scale = torch.ldexp(torch.ones_like(amax), exp)
+    scale = torch.where(amax > 0, scale, torch.ones_like(scale))
+
+    y = (w32 / scale).abs().clamp(max=_E2M1_MAX)
+    q = round_to_e2m1_grid(torch.where(nonfinite, torch.zeros_like(y), y))
+    q = torch.copysign(q, w32)
+    out = q * scale
+    out = torch.where(nonfinite.expand_as(out), torch.full_like(out, float("nan")), out)
+    return out.view(rows, cols).to(weight.dtype)
diff --git a/tests/pytorch/test_mxfp4_qat.py b/tests/pytorch/test_mxfp4_qat.py
new file mode 100644
index 0000000000..e304b49a52
--- /dev/null
+++ b/tests/pytorch/test_mxfp4_qat.py
@@ -0,0 +1,905 @@
+# Copyright (c) 2022-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+#
+# See LICENSE for license information.
+
+"""Per-step and end-to-end tests for the MXFP4 weight-QAT recipes."""
+import ctypes
+import os
+
+for _n in ("libnvrtc.so.13", "libnvrtc.so.12", "libnvrtc.so"):
+    try:
+        ctypes.CDLL(_n, mode=ctypes.RTLD_GLOBAL)
+        break
+    except OSError:
+        continue
+
+import torch
+
+import transformer_engine.pytorch as te
+import transformer_engine_torch as tex
+from transformer_engine.common.recipe import (
+    MXFP4QATFloat8BlockScaling,
+    MXFP4QATMXFP8BlockScaling,
+)
+from transformer_engine.pytorch import fp8_autocast
+from transformer_engine.pytorch.mxfp4_qat import mxfp4_fake_quantize
+from references.mxfp4_qat_reference import mxfp4_fake_quantize_reference
+from transformer_engine.pytorch.tensor.float8_blockwise_tensor import Float8BlockQuantizer
+from transformer_engine.pytorch.tensor.mxfp8_tensor import MXFP8Quantizer
+
+torch.manual_seed(1234)
+DEV = "cuda"
+RUN_PERF = os.environ.get("NVTE_MXFP4_QAT_PERF", "0") == "1"
+E2M1_GRID = torch.tensor([0.0, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0], device=DEV)
+SHAPES = [(128, 256), (192, 448), (512, 1024)]
+
+FP8_TOL_1OP = 0.045
+FP8_TOL_2OP = 0.065
+BF16_TOL = 0.005
+
+
+def make_weight(m, n, dtype=torch.bfloat16, zero_blocks=0, outliers=0):
+    w = (torch.randn(m, n, device=DEV, dtype=torch.float32) * 0.02).to(dtype)
+    if outliers:
+        idx = torch.randint(0, m * n, (outliers,), device=DEV)
+        w.view(-1)[idx] *= 500.0
+    if zero_blocks:
+        for i in range(zero_blocks):
+            r = torch.randint(0, m, (1,)).item()
+            k = torch.randint(0, n // 32, (1,)).item()
+            w[r, k * 32 : (k + 1) * 32] = 0.0
+    return w
+
+
+def rel_err(got, ref):
+    got, ref = got.to(torch.float64), ref.to(torch.float64)
+    return ((got - ref).norm() / ref.norm().clamp(min=1e-30)).item()
+
+
+def fake_quant_ref_fp64(w):
+    """Independent fp64 reference: nearest-with-ties-to-even onto the E2M1 grid."""
+    m, n = w.shape
+    w64 = w.to(torch.float64).view(m, n // 32, 32)
+    amax = w64.abs().amax(dim=-1, keepdim=True)
+    bits = amax.to(torch.float32).view(torch.int32)
+    exp_field = bits >> 23
+    mant = bits & 0x7FFFFF
+    exp = exp_field - 129 + (mant > 0x400000).to(torch.int32)
+    exp = torch.where(exp_field > 0, exp, torch.full_like(exp, -126)).clamp(-126, 125)
+    exp = torch.where(amax > 0, exp, torch.zeros_like(exp))
+    scale = torch.ldexp(torch.ones_like(amax), exp)
+    y = (w64 / scale).clamp(min=-6.0, max=6.0)
+    grid = E2M1_GRID.to(torch.float64)
+    cand = torch.cat([-grid.flip(0)[:-1], grid])
+    d = (y.unsqueeze(-1) - cand).abs()
+    near = d.argmin(dim=-1)
+    even = torch.tensor([abs(v) in (0.0, 1.0, 2.0, 4.0) for v in cand.tolist()], device=w.device)
+    dmin = d.gather(-1, near.unsqueeze(-1)).squeeze(-1)
+    tie_even = d.eq(dmin.unsqueeze(-1)) & even
+    has_even_tie = tie_even.any(dim=-1)
+    near_even = torch.where(has_even_tie, tie_even.to(torch.uint8).argmax(dim=-1), near)
+    q = cand[near_even]
+    return (q * scale).view(m, n).to(torch.float64)
+
+
+def test_fake_quant_grid_scale_rtne():
+    for m, n in SHAPES:
+        w = make_weight(m, n, zero_blocks=4, outliers=4)
+        w_hat = mxfp4_fake_quantize(w)
+        assert w_hat.dtype == w.dtype
+
+        wh32 = w_hat.to(torch.float32).view(m, n // 32, 32)
+        amax = w.to(torch.float32).abs().view(m, n // 32, 32).amax(dim=-1, keepdim=True)
+        bits = amax.view(torch.int32)
+        exp_field = bits >> 23
+        mant = bits & 0x7FFFFF
+        exp = exp_field - 129 + (mant > 0x400000).to(torch.int32)
+        exp = torch.where(exp_field > 0, exp, torch.full_like(exp, -126)).clamp(-126, 125)
+        scale = torch.ldexp(torch.ones_like(amax), exp)
+        scale = torch.where(amax > 0, scale, torch.ones_like(scale))
+
+        vals = (wh32 / scale).abs()
+        assert torch.isin(vals, E2M1_GRID).all(), "values off the E2M1 grid"
+        nz = (amax.double() > 6.0 * 2.0**-126) & (amax.double() <= 6.0 * 2.0**125)
+        assert ((scale.double() * 6.0)[nz] >= amax.double()[nz]).all()
+        assert (((scale.double() / 2.0) * 6.0)[nz] < amax.double()[nz]).all()
+        ref = fake_quant_ref_fp64(w)
+        assert torch.equal(w_hat.to(torch.float64), ref), f"RTNE mismatch {m}x{n}"
+        assert torch.equal(mxfp4_fake_quantize(w_hat), w_hat), "not idempotent"
+
+
+def test_D_lossless_in_bf16():
+    for m, n in SHAPES:
+        w = make_weight(m, n, outliers=4)
+        w_hat_bf16 = mxfp4_fake_quantize(w)
+        w_hat_fp32 = mxfp4_fake_quantize(w.to(torch.float32))
+        assert torch.equal(
+            w_hat_bf16.to(torch.float32), w_hat_fp32
+        ), "MXFP4-grid weight is not exact in bf16"
+
+
+def test_E_mxfp8_rowwise_lossless():
+    for m, n in SHAPES:
+        w = make_weight(m, n, zero_blocks=4, outliers=4)
+        w[0, :32] = torch.tensor(2.0**-127, dtype=torch.bfloat16)
+        w[0, 0] = 2.0**-128
+        w_hat = mxfp4_fake_quantize(w)
+        q = MXFP8Quantizer(fp8_dtype=tex.DType.kFloat8E4M3, columnwise=False).quantize(w_hat)
+        dq = q.dequantize(dtype=torch.float32)
+        wh32 = w_hat.to(torch.float32)
+
+        data = q._rowwise_data.view(torch.uint8)[:m, :n]
+        codes = q._rowwise_scale_inv.view(torch.uint8)[:m, : n // 32].to(torch.int32)
+        raw_vals = data.view(torch.float8_e4m3fn).to(torch.float32).view(m, n // 32, 32)
+        raw_scale = torch.ldexp(torch.ones_like(codes, dtype=torch.float32), codes - 127)
+        raw_dq = (raw_vals * raw_scale.unsqueeze(-1)).view(m, n)
+        assert torch.equal(
+            raw_dq, wh32
+        ), f"RAW MXFP8 encoding of the MXFP4-grid weight is not lossless {m}x{n}"
+
+        assert torch.equal(
+            dq, wh32
+        ), f"MXFP8 rowwise encoding of the MXFP4-grid weight is not lossless {m}x{n}"
+
+
+def test_dequantize_e8m0_extreme_codes():
+    """Real MXFP8 software dequantize with planted UE8M0 codes 0/255."""
+    w = make_weight(32, 64)
+    q = MXFP8Quantizer(fp8_dtype=tex.DType.kFloat8E4M3, columnwise=False).quantize(
+        mxfp4_fake_quantize(w)
+    )
+    data = q._rowwise_data.view(torch.uint8)
+    codes = q._rowwise_scale_inv.view(torch.uint8)
+    data[0, :32] = 56
+    codes[0, 0] = 0
+    data[1, :32] = 56
+    codes[1, 0] = 255
+    dq = q.dequantize(dtype=torch.float32)
+    assert (dq[0, :32] == 2.0**-127).all(), "UE8M0 code 0 must dequantize to 2^-127"
+    assert torch.isnan(dq[1, :32]).all(), "UE8M0 code 255 must dequantize to NaN"
+
+
+def test_G_blockwise_lossless_and_transpose():
+    for m, n in SHAPES:
+        w = make_weight(m, n, zero_blocks=4)
+        w_hat = mxfp4_fake_quantize(w)
+        q = Float8BlockQuantizer(
+            fp8_dtype=tex.DType.kFloat8E4M3,
+            rowwise=True,
+            columnwise=True,
+            force_pow_2_scales=True,
+            block_scaling_dim=2,
+        ).quantize(w_hat)
+        dq = q.dequantize(dtype=torch.float32)
+        assert torch.equal(
+            dq, w_hat.to(torch.float32)
+        ), f"128x128 blockwise encoding of the MXFP4-grid weight is not lossless {m}x{n}"
+        q.update_usage(rowwise_usage=False, columnwise_usage=True)
+        dqc = q.dequantize(dtype=torch.float32)
+        assert torch.equal(
+            dqc, w_hat.to(torch.float32)
+        ), f"columnwise 128x128 dequant differs from the MXFP4-grid weight {m}x{n}"
+
+
+def test_lossless_dequant_to_bf16():
+    """Both weight encodings dequantize to bf16 bit-exactly (bwd override target dtype)."""
+    for m, n in SHAPES:
+        w = make_weight(m, n, zero_blocks=4, outliers=4)
+        w[0, :32] = torch.tensor(2.0**-127, dtype=torch.bfloat16)
+        w_hat = mxfp4_fake_quantize(w)
+
+        q8 = MXFP8Quantizer(fp8_dtype=tex.DType.kFloat8E4M3, columnwise=False).quantize(w_hat)
+        dq8 = q8.dequantize(dtype=torch.bfloat16)
+        assert torch.equal(dq8, w_hat), f"mxfp8 rowwise dequant to bf16 not lossless {m}x{n}"
+
+        w2_hat = mxfp4_fake_quantize(make_weight(m, n, zero_blocks=4))
+        qb = Float8BlockQuantizer(
+            fp8_dtype=tex.DType.kFloat8E4M3,
+            rowwise=True,
+            columnwise=True,
+            force_pow_2_scales=True,
+            block_scaling_dim=2,
+        ).quantize(w2_hat)
+        assert torch.equal(
+            qb.dequantize(dtype=torch.bfloat16), w2_hat
+        ), f"blockwise 128x128 dequant to bf16 not lossless {m}x{n}"
+
+
+def test_C_mxfp8_colwise_bound():
+    for m, n in SHAPES:
+        w_hat = mxfp4_fake_quantize(make_weight(m, n, zero_blocks=4, outliers=4))
+        q = MXFP8Quantizer(fp8_dtype=tex.DType.kFloat8E4M3).quantize(w_hat)
+        q.update_usage(rowwise_usage=False, columnwise_usage=True)
+        dq = q.dequantize(dtype=torch.float32)
+        ref = w_hat.to(torch.float32)
+        blk = ref.view(m // 32, 32, n)
+        bound = blk.abs().amax(dim=1, keepdim=True) * 2.0**-3
+        err = (dq.view(m // 32, 32, n) - blk).abs()
+        assert (err <= bound + 1e-30).all(), f"colwise 32x1 requant error above bound {m}x{n}"
+
+
+_INTREE_KERNELS = {}
+
+
+def _intree_kernel(fast_math=False):
+    """JIT-compile the in-tree kernel; fast_math=True builds the same source with --use_fast_math."""
+    if fast_math not in _INTREE_KERNELS:
+        import os
+        from torch.utils.cpp_extension import load_inline
+
+        src_dir = os.environ.get(
+            "NVTE_MXFP4_QAT_TEST_SRC",
+            os.path.join(
+                os.path.dirname(os.path.abspath(__file__)),
+                "..",
+                "..",
+                "transformer_engine",
+                "common",
+            ),
+        )
+        cuda_src = r"""
+#include <cuda.h>
+#include <torch/extension.h>
+#include <ATen/cuda/CUDAContext.h>
+#include "cast/mxfp4/fake_quantize_mxfp4.cuh"
+
+torch::Tensor mxfp4_fake_quantize_intree(torch::Tensor w) {
+  TORCH_CHECK(w.is_cuda() && w.is_contiguous() && w.numel() % 32 == 0);
+  auto out = torch::empty_like(w);
+  auto stream = at::cuda::getCurrentCUDAStream();
+  if (w.scalar_type() == at::kBFloat16) {
+    transformer_engine::dispatch::mxfp4::fake_quantize_mxfp4_launch<nv_bfloat16>(
+        reinterpret_cast<const nv_bfloat16 *>(w.data_ptr()),
+        reinterpret_cast<nv_bfloat16 *>(out.data_ptr()), w.numel(), stream);
+  } else if (w.scalar_type() == at::kFloat) {
+    transformer_engine::dispatch::mxfp4::fake_quantize_mxfp4_launch<float>(
+        w.data_ptr<float>(), out.data_ptr<float>(), w.numel(), stream);
+  } else {
+    TORCH_CHECK(false, "bf16/fp32 only");
+  }
+  return out;
+}
+"""
+        _INTREE_KERNELS[fast_math] = load_inline(
+            name="nvte_mxfp4_qat_intree_test" + ("_fastmath" if fast_math else ""),
+            cpp_sources="torch::Tensor mxfp4_fake_quantize_intree(torch::Tensor w);",
+            cuda_sources=cuda_src,
+            functions=["mxfp4_fake_quantize_intree"],
+            extra_include_paths=[src_dir],
+            extra_cuda_cflags=[
+                "-O3",
+                "-U__CUDA_NO_HALF_OPERATORS__",
+                "-U__CUDA_NO_HALF_CONVERSIONS__",
+                "-U__CUDA_NO_HALF2_OPERATORS__",
+                "-U__CUDA_NO_BFLOAT16_OPERATORS__",
+                "-U__CUDA_NO_BFLOAT16_CONVERSIONS__",
+            ]
+            + (["--use_fast_math"] if fast_math else []),
+            verbose=False,
+        )
+    return _INTREE_KERNELS[fast_math]
+
+
+def _both_paths(w):
+    got_kernel = _intree_kernel().mxfp4_fake_quantize_intree(w.contiguous())
+    got_torch = mxfp4_fake_quantize_reference(w)
+    return got_kernel, got_torch
+
+
+def _assert_bits_equal(a, b, msg):
+    assert a.dtype == b.dtype, f"{msg}: dtype mismatch"
+    ib = torch.int16 if a.dtype == torch.bfloat16 else torch.int32
+    assert torch.equal(a.view(ib), b.view(ib)), f"{msg}: raw bits differ"
+
+
+def _assert_same_with_nan(a, b, msg):
+    an, bn = torch.isnan(a), torch.isnan(b)
+    assert torch.equal(an, bn), f"{msg}: NaN masks differ"
+    _assert_bits_equal(
+        torch.where(an, torch.zeros_like(a), a),
+        torch.where(bn, torch.zeros_like(b), b),
+        msg,
+    )
+
+
+PARITY_SHAPES = SHAPES + [(3, 416), (1, 32)]
+
+
+def test_kernel_matches_torch_reference():
+    for m, n in PARITY_SHAPES:
+        for dtype in (torch.bfloat16, torch.float32):
+            w = make_weight(m, n, dtype=dtype, zero_blocks=4, outliers=4)
+            gk, gt = _both_paths(w)
+            _assert_bits_equal(gk, gt, f"kernel != torch reference {m}x{n} {dtype}")
+
+
+def test_edge_value_domains():
+    m, n = 64, 128
+    w = make_weight(m, n)
+
+    z = torch.zeros(m, n, device=DEV, dtype=torch.bfloat16)
+    gk, gt = _both_paths(z)
+    assert torch.equal(gk, z) and torch.equal(gt, z)
+
+    wz = w.clone()
+    wz[0, :32] = -0.0
+    gk, gt = _both_paths(wz)
+    assert (gk[0, :32] == 0).all() and torch.signbit(gk[0, :32]).all(), "-0.0 sign lost"
+    _assert_bits_equal(gk, gt, "-0 block")
+
+    ws = w.clone()
+    ws[1, :32] = torch.tensor(2.0**-133, dtype=torch.bfloat16)
+    gk, gt = _both_paths(ws)
+    assert torch.isfinite(gk[1, :32]).all() and torch.equal(gk, gt)
+    ref = fake_quant_ref_fp64(ws)
+    assert torch.equal(gt.to(torch.float64), ref), "subnormal block deviates from fp64 ref"
+
+    wp = w.clone()
+    wp[5, :32] = torch.tensor(2.0**-127, dtype=torch.bfloat16)
+    gk, gt = _both_paths(wp)
+    assert (gk[5, :32].to(torch.float32) == 2.0**-127).all(), "2^-127 grid point lost"
+    _assert_bits_equal(gk, gt, "2^-127 grid point")
+    ref = fake_quant_ref_fp64(wp)
+    assert torch.equal(gt.to(torch.float64), ref), "2^-127 block deviates from fp64 ref"
+
+    wh = w.clone()
+    wh[2, :32] = torch.tensor(3.0e38, dtype=torch.bfloat16)
+    wh[2, 0] = -3.0e38
+    gk, gt = _both_paths(wh)
+    assert torch.isfinite(gk[2, :32]).all() and torch.equal(gk, gt)
+    ref = fake_quant_ref_fp64(wh)
+    assert torch.equal(gt.to(torch.float64), ref), "huge block deviates from fp64 ref"
+
+    for bad in (float("inf"), float("-inf"), float("nan")):
+        wb = w.clone()
+        wb[3, 40] = bad
+        gk, gt = _both_paths(wb)
+        _assert_same_with_nan(gk, gt, f"bad={bad}")
+        assert torch.isnan(gk[3, 32:64]).all(), f"{bad}: block not NaN-poisoned"
+        assert torch.isfinite(gk[3, :32]).all(), f"{bad}: neighbor block affected"
+        clean_rows = torch.ones(m, dtype=torch.bool)
+        clean_rows[3] = False
+        assert torch.isfinite(gk[clean_rows]).all()
+
+    try:
+        mxfp4_fake_quantize(torch.zeros(32, 32, device=DEV, dtype=torch.float16))
+        raise SystemExit("fp16 should be rejected")
+    except ValueError:
+        pass
+
+
+def test_blockwise_recipe_is_128x128():
+    r = MXFP4QATFloat8BlockScaling()
+    assert r.x_block_scaling_dim == 1 and r.w_block_scaling_dim == 2
+    assert r.grad_block_scaling_dim == 1
+    q = Float8BlockQuantizer(
+        fp8_dtype=tex.DType.kFloat8E4M3,
+        rowwise=True,
+        columnwise=True,
+        force_pow_2_scales=True,
+        block_scaling_dim=2,
+    )
+    assert q.block_len == 128
+    w_hat = mxfp4_fake_quantize(make_weight(256, 512))
+    qt = q.quantize(w_hat)
+    assert tuple(qt._rowwise_scale_inv.shape)[0] == 256 // 128
+
+
+def test_blockwise_over_ratio_is_bounded_not_lossless():
+    w = make_weight(128, 128)
+    w[0, :32] = torch.tensor(2.0**-14, dtype=torch.bfloat16)
+    w[64, 0] = 100.0
+    w_hat = mxfp4_fake_quantize(w)
+    q = Float8BlockQuantizer(
+        fp8_dtype=tex.DType.kFloat8E4M3,
+        rowwise=True,
+        columnwise=True,
+        force_pow_2_scales=True,
+        block_scaling_dim=2,
+    ).quantize(w_hat)
+    dq = q.dequantize(dtype=torch.float32)
+    assert torch.isfinite(dq).all()
+    err = (dq - w_hat.to(torch.float32)).abs().max().item()
+    tile_amax = w_hat.to(torch.float32).abs().max().item()
+    assert err <= tile_amax * 2.0**-9, "over-ratio loss larger than FP8 headroom"
+
+
+def test_kernel_perf():
+    if not RUN_PERF:
+        print("  SKIP (set NVTE_MXFP4_QAT_PERF=1 to run benchmarks)")
+        return
+    rows, cols = 8192, 8192
+    w = make_weight(rows, cols)
+    k = _intree_kernel()
+    for _ in range(3):
+        k.mxfp4_fake_quantize_intree(w)
+        mxfp4_fake_quantize_reference(w)
+    torch.cuda.synchronize()
+
+    def timeit(fn, iters=20):
+        start, end = torch.cuda.Event(True), torch.cuda.Event(True)
+        start.record()
+        for _ in range(iters):
+            fn()
+        end.record()
+        torch.cuda.synchronize()
+        return start.elapsed_time(end) / iters * 1e3
+
+    t_kernel = timeit(lambda: k.mxfp4_fake_quantize_intree(w))
+    t_torch = timeit(lambda: mxfp4_fake_quantize_reference(w))
+    moved = rows * cols * 2 * 2
+    print(
+        f"  kernel: {t_kernel:.1f} us ({moved / t_kernel / 1e3:.0f} GB/s), "
+        f"torch: {t_torch:.1f} us ({moved / t_torch / 1e3:.0f} GB/s)"
+    )
+    assert t_kernel < t_torch, "kernel slower than the torch reference"
+
+
+def test_kernel_fast_math_immune():
+    """The kernel stays bit-identical when compiled with --use_fast_math."""
+    k = _intree_kernel(fast_math=True)
+    for m, n in PARITY_SHAPES:
+        for dtype in (torch.bfloat16, torch.float32):
+            w = make_weight(m, n, dtype=dtype, zero_blocks=4, outliers=4)
+            got = k.mxfp4_fake_quantize_intree(w.contiguous())
+            ref = mxfp4_fake_quantize_reference(w)
+            _assert_bits_equal(got, ref, f"fast-math build diverged {m}x{n} {dtype}")
+
+    w = make_weight(64, 128)
+    w[0, :32] = 0.0
+    w[1, 5] = float("inf")
+    w[2, 9] = float("nan")
+    w[3, :32] = torch.tensor(2.0**-133, dtype=torch.bfloat16)
+    w[4, :32] = torch.tensor(3.0e38, dtype=torch.bfloat16)
+    w[5, :32] = torch.tensor(2.0**-127, dtype=torch.bfloat16)
+    got = k.mxfp4_fake_quantize_intree(w.contiguous())
+    ref = mxfp4_fake_quantize_reference(w)
+    _assert_same_with_nan(got, ref, "fast-math edge")
+    assert (
+        got[5, :32].to(torch.float32) == 2.0**-127
+    ).all(), "fast-math build flushed the 2^-127 grid point"
+
+
+def test_exp2f_e8m0_all_codes():
+    """All 256 UE8M0 codes through exp2f/exp2f_rcp extracted verbatim from ptx.cuh."""
+    import re as _re
+    from torch.utils.cpp_extension import load_inline
+
+    src_dir = os.environ.get(
+        "NVTE_MXFP4_QAT_TEST_SRC",
+        os.path.join(
+            os.path.dirname(os.path.abspath(__file__)), "..", "..", "transformer_engine", "common"
+        ),
+    )
+    ptx_src = open(os.path.join(src_dir, "util", "ptx.cuh")).read()
+    m_exp = _re.search(
+        r"__device__ __forceinline__ float exp2f\(e8m0_t biased_exp\) \{.*?\n\}", ptx_src, _re.S
+    )
+    m_rcp = _re.search(
+        r"template <>\n__device__ __forceinline__ float exp2f_rcp<float>\(e8m0_t biased_exp\)"
+        r" \{.*?\n\}",
+        ptx_src,
+        _re.S,
+    )
+    assert m_exp and m_rcp, "could not extract exp2f/exp2f_rcp from ptx.cuh"
+    cuda_src = (
+        "#include <cuda.h>\n#include <torch/extension.h>\n"
+        "#include <ATen/cuda/CUDAContext.h>\n"
+        "namespace teptx {\n"
+        "typedef unsigned char e8m0_t;\n"
+        "constexpr unsigned FP32_MANTISSA_BITS = 23;\n"
+        "template <typename T> __device__ __forceinline__ T exp2f_rcp(e8m0_t);\n"
+        + m_rcp.group(0)
+        + "\n"
+        + m_exp.group(0)
+        + "\n}\n"
+        "__global__ void all_codes_kernel(float *e, float *r) {\n"
+        "  const int b = threadIdx.x;\n"
+        "  e[b] = teptx::exp2f(static_cast<teptx::e8m0_t>(b));\n"
+        "  r[b] = teptx::exp2f_rcp<float>(static_cast<teptx::e8m0_t>(b));\n"
+        "}\n"
+        "void run_all_codes(torch::Tensor e, torch::Tensor r) {\n"
+        "  all_codes_kernel<<<1, 256, 0, at::cuda::getCurrentCUDAStream()>>>(\n"
+        "      e.data_ptr<float>(), r.data_ptr<float>());\n"
+        "}\n"
+    )
+    mod = load_inline(
+        name="nvte_mxfp4_qat_exp2f_test",
+        cpp_sources="void run_all_codes(torch::Tensor e, torch::Tensor r);",
+        cuda_sources=cuda_src,
+        functions=["run_all_codes"],
+        verbose=False,
+    )
+    e = torch.empty(256, dtype=torch.float32, device=DEV)
+    r = torch.empty_like(e)
+    mod.run_all_codes(e, r)
+    torch.cuda.synchronize()
+    codes = torch.arange(256, dtype=torch.int32)
+    exp_ref = torch.ldexp(torch.ones(256), codes - 127)
+    rcp_ref = torch.ldexp(torch.ones(256), 127 - codes)
+    assert torch.isnan(e[255]) and torch.isnan(r[255]), "code 255 must be NaN"
+    assert e.view(torch.int32)[255].item() == 0x7FFFFFFF, "code 255 NaN payload"
+    assert r.view(torch.int32)[255].item() == 0x7FFFFFFF, "code 255 rcp NaN payload"
+    assert torch.equal(e[:255].cpu(), exp_ref[:255]), "exp2f wrong (code 0 == 2^-127?)"
+    assert torch.equal(r[:255].cpu(), rcp_ref[:255]), "exp2f_rcp wrong"
+
+
+def test_ste_identity_gradient():
+    """mxfp4_fake_quantize carries an identity (STE) gradient and matches the reference forward."""
+    for dtype in (torch.bfloat16, torch.float32):
+        w = make_weight(64, 128, dtype=dtype).requires_grad_(True)
+        up = torch.randn_like(w)
+        out = mxfp4_fake_quantize(w)
+        out.backward(up)
+        assert torch.equal(w.grad, up), f"STE gradient is not identity ({dtype})"
+        _assert_bits_equal(out.detach(), mxfp4_fake_quantize_reference(w.detach()), "STE fwd")
+
+
+def test_misaligned_and_noncontiguous():
+    """Misaligned storage-offset views and non-contiguous inputs match an aligned copy bitwise."""
+    for dtype, off in ((torch.bfloat16, 4), (torch.float32, 2)):
+        m, n = 64, 128
+        flat = (torch.randn(m * n + off, device=DEV, dtype=torch.float32) * 0.02).to(dtype)
+        v = flat[off:].view(m, n)
+        assert v.is_contiguous() and v.data_ptr() % 16 != 0
+        got = mxfp4_fake_quantize(v)
+        ref = mxfp4_fake_quantize(v.clone())
+        _assert_bits_equal(got, ref, f"misaligned view {dtype}")
+
+    nc = (torch.randn(64, 256, device=DEV, dtype=torch.float32) * 0.02)[:, ::2]
+    assert not nc.is_contiguous()
+    got = mxfp4_fake_quantize_reference(nc)
+    assert torch.equal(got, mxfp4_fake_quantize_reference(nc.contiguous()))
+
+
+def test_bf16_exhaustive_bit_patterns():
+    """All 65536 bf16 bit patterns: kernel/torch parity, fp64 oracle on finite blocks."""
+    bits = torch.arange(65536, dtype=torch.int32, device=DEV).to(torch.int16)
+    w = bits.view(torch.bfloat16).view(2048, 32)
+    gk, gt = _both_paths(w)
+    _assert_same_with_nan(gk, gt, "bf16 exhaustive")
+    finite = torch.isfinite(w.to(torch.float32)).all(dim=-1)
+    ref = fake_quant_ref_fp64(w[finite])
+    assert torch.equal(gt[finite].to(torch.float64), ref), "bf16 exhaustive vs fp64 ref"
+
+
+def test_fp32_bit_fuzz():
+    """Random raw fp32 bit patterns (subnormals/NaN/Inf included)."""
+    g = torch.Generator().manual_seed(99)
+    for _ in range(4):
+        bits = torch.randint(-(2**31), 2**31 - 1, (4096, 32), generator=g, dtype=torch.int64)
+        w = bits.to(torch.int32).cuda().view(torch.float32)
+        gk, gt = _both_paths(w)
+        _assert_same_with_nan(gk, gt, "fp32 fuzz")
+        finite = torch.isfinite(w).all(dim=-1)
+        if finite.any():
+            ref = fake_quant_ref_fp64(w[finite])
+            assert torch.equal(gt[finite].to(torch.float64), ref), "fp32 fuzz vs fp64 ref"
+
+
+def test_scale_threshold_and_rtne_midpoints():
+    """amax = 1.5*2^k ulp triples and all E2M1 RTNE midpoints vs explicit expected values."""
+    rows = []
+    for k in (-100, -10, 0, 42, 100):
+        a = torch.tensor(1.5 * 2.0**k, dtype=torch.float32)
+        ab = a.view(torch.int32)
+        for delta in (-1, 0, 1):
+            row = torch.zeros(32, dtype=torch.float32)
+            row[0] = (ab + delta).view(torch.float32)
+            rows.append(row)
+    w = torch.stack(rows).cuda()
+    gk, gt = _both_paths(w)
+    _assert_bits_equal(gk, gt, "threshold parity")
+    assert torch.equal(gt.to(torch.float64), fake_quant_ref_fp64(w)), "threshold vs fp64 ref"
+
+    mids = torch.tensor([0.25, 0.75, 1.25, 1.75, 2.5, 3.5, 5.0])
+    expected = torch.tensor([0.0, 1.0, 1.0, 2.0, 2.0, 4.0, 4.0])
+    for t in (-60, 0, 60):
+        sc = 2.0**t
+        block = torch.zeros(32)
+        block[:7] = mids * sc
+        block[7] = 6.0 * sc
+        w = block.view(1, 32).to(torch.float32).cuda()
+        gk, gt = _both_paths(w)
+        _assert_bits_equal(gk, gt, f"midpoints t={t}")
+        exp_row = torch.zeros(32)
+        exp_row[:7] = expected * sc
+        exp_row[7] = 6.0 * sc
+        assert torch.equal(
+            gt[0].cpu().to(torch.float32), exp_row.to(torch.float32)
+        ), f"midpoint values t={t}"
+
+
+def test_deployment_top_cap_domain():
+    """amax <= 6*2^125 is deployment-bitwise; above, satfinite at the cap while TileKernels moves to scale 2^126."""
+    s125 = 6.0 * 2.0**125
+    w = torch.zeros(2, 32, dtype=torch.float32)
+    w[0, 0] = s125
+    just_above = (torch.tensor(s125, dtype=torch.float32).view(torch.int32) + 1).view(torch.float32)
+    w[1, 0] = just_above
+    wc = w.cuda()
+    gk, gt = _both_paths(wc)
+    _assert_bits_equal(gk, gt, "top-cap parity")
+    assert gk[0, 0].item() == s125, "boundary block must be exact"
+    assert gk[1, 0].item() == s125, "above-boundary block must satfinite at 6*2^125"
+    assert torch.equal(gt.to(torch.float64), fake_quant_ref_fp64(wc)), "top-cap vs fp64 ref"
+
+
+def test_blockwise_feasibility_enumeration():
+    """Per-element E4M3-lattice prediction == TE 128x128 quantizer for exponent spreads d=0..16."""
+    payloads = torch.tensor([0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0])
+    boundary = None
+    for d in range(17):
+        w = torch.zeros(128, 128, dtype=torch.float32)
+        w[0, 0] = 6.0
+        vals = payloads * 2.0 ** (-d)
+        w[1, :7] = vals
+        w[1, 7:14] = -vals
+        wc = w.cuda()
+        assert torch.equal(mxfp4_fake_quantize(wc), wc), f"tile not on MXFP4 grid (d={d})"
+        q = Float8BlockQuantizer(
+            fp8_dtype=tex.DType.kFloat8E4M3,
+            rowwise=True,
+            columnwise=False,
+            force_pow_2_scales=True,
+            block_scaling_dim=2,
+        ).quantize(wc)
+        S = q._rowwise_scale_inv.reshape(-1)[0].item()
+        dq = q.dequantize(dtype=torch.float32)
+        pred = ((wc / S).to(torch.float8_e4m3fn).to(torch.float32) * S) == wc
+        got = dq == wc
+        assert torch.equal(pred, got), f"E4M3-lattice prediction mismatch at spread d={d}"
+        if boundary is None and not bool(pred.all()):
+            boundary = d
+    assert (
+        boundary is not None and boundary >= 12
+    ), f"full-payload exactness should hold at least through spread 2^11, got {boundary}"
+    print(
+        f"  full-grid exact through spread d={boundary - 1}; first loss at d={boundary} "
+        "(E4M3 subnormal lattice bound)"
+    )
+
+
+def test_fp16_pipeline_rejected():
+    """fp16 weights and fp16 activation dtypes must fail loudly under MXFP4 QAT."""
+    mod = te.Linear(256, 128, bias=False, params_dtype=torch.float16, device=DEV)
+    x = torch.randn(64, 256, device=DEV, dtype=torch.float16)
+    for recipe in (MXFP4QATMXFP8BlockScaling(), MXFP4QATFloat8BlockScaling()):
+        try:
+            with fp8_autocast(enabled=True, fp8_recipe=recipe):
+                mod(x)
+            raise SystemExit(f"fp16 pipeline must be rejected ({recipe.__class__.__name__})")
+        except (NotImplementedError, ValueError):
+            pass
+
+
+def test_recipe_switch_invalidates_cache():
+    """base<->QAT recipe switches must requantize instead of reusing cached weight workspaces."""
+    import warnings as _warnings
+    from transformer_engine.common.recipe import MXFP8BlockScaling
+
+    torch.manual_seed(3)
+    x = torch.randn(64, 256, device=DEV, dtype=torch.bfloat16)
+    mod = te.Linear(256, 128, bias=False, params_dtype=torch.bfloat16, device=DEV)
+    base, qat = MXFP8BlockScaling(), MXFP4QATMXFP8BlockScaling()
+
+    with _warnings.catch_warnings():
+        _warnings.simplefilter("ignore")
+        with fp8_autocast(enabled=True, fp8_recipe=base):
+            mod(x, is_first_microbatch=True)
+        with fp8_autocast(enabled=True, fp8_recipe=qat):
+            y_switch = mod(x, is_first_microbatch=False)
+            y_fresh = mod(x, is_first_microbatch=True)
+        assert torch.equal(
+            y_switch, y_fresh
+        ), "base->QAT switch reused an UNPROJECTED cached weight"
+        with fp8_autocast(enabled=True, fp8_recipe=base):
+            y2_switch = mod(x, is_first_microbatch=False)
+            y2_fresh = mod(x, is_first_microbatch=True)
+        assert torch.equal(y2_switch, y2_fresh), "QAT->base switch reused a PROJECTED cached weight"
+
+
+def test_recipe_field_controls_qat():
+    """mxfp4_qat_weights field (env NVTE_MXFP4_QAT) turns base recipes into QAT recipes."""
+    import warnings as _warnings
+    from transformer_engine.common.recipe import Float8BlockScaling, MXFP8BlockScaling
+
+    base = MXFP8BlockScaling(mxfp4_qat_weights=True)
+    assert base.mxfp4_qat() and MXFP4QATMXFP8BlockScaling().mxfp4_qat()
+    assert not MXFP8BlockScaling().mxfp4_qat()
+
+    torch.manual_seed(9)
+    x = torch.randn(64, 256, device=DEV, dtype=torch.bfloat16)
+    mod = te.Linear(256, 128, bias=False, params_dtype=torch.bfloat16, device=DEV)
+    with _warnings.catch_warnings():
+        _warnings.simplefilter("ignore")
+        with fp8_autocast(enabled=True, fp8_recipe=base):
+            y_field = mod(x, is_first_microbatch=True)
+        with fp8_autocast(enabled=True, fp8_recipe=MXFP4QATMXFP8BlockScaling()):
+            y_class = mod(x, is_first_microbatch=True)
+    assert torch.equal(y_field, y_class), "field-driven QAT differs from class-driven QAT"
+
+    assert Float8BlockScaling(mxfp4_qat_weights=True).mxfp4_qat()
+    try:
+        Float8BlockScaling(mxfp4_qat_weights=True, w_block_scaling_dim=1)
+        raise SystemExit("field-driven blockwise QAT must validate 2D weight blocks")
+    except ValueError:
+        pass
+
+
+def test_ops_api_rejected():
+    """te.ops bypasses the QAT weight hook and must reject MXFP4 QAT recipes."""
+    import transformer_engine.pytorch.ops as te_ops
+
+    op = te_ops.Linear(256, 128, bias=False, device=DEV, dtype=torch.bfloat16)
+    x = torch.randn(64, 256, device=DEV, dtype=torch.bfloat16)
+    for recipe in (MXFP4QATMXFP8BlockScaling(), MXFP4QATFloat8BlockScaling()):
+        try:
+            with fp8_autocast(enabled=True, fp8_recipe=recipe):
+                op(x)
+            raise SystemExit(f"ops API must reject {recipe.__class__.__name__}")
+        except NotImplementedError:
+            pass
+
+
+def _run_module(module_kind, recipe, override, fuse):
+    torch.manual_seed(7)
+    m_splits = [192, 320]
+    tokens, in_f, out_f = sum(m_splits), 256, 512
+    recipe.backward_override = override
+
+    if module_kind == "linear":
+        mod = te.Linear(
+            in_f,
+            out_f,
+            bias=False,
+            params_dtype=torch.bfloat16,
+            device=DEV,
+            fuse_wgrad_accumulation=fuse,
+        )
+        params = [mod.weight]
+    else:
+        mod = te.GroupedLinear(
+            2,
+            in_f,
+            out_f,
+            bias=False,
+            params_dtype=torch.bfloat16,
+            device=DEV,
+            fuse_wgrad_accumulation=fuse,
+        )
+        params = [mod.weight0, mod.weight1]
+    with torch.no_grad():
+        for p in params:
+            p.copy_((torch.randn_like(p, dtype=torch.float32) * 0.02).to(p.dtype))
+    w_hats = [mxfp4_fake_quantize(p.detach()) for p in params]
+    if fuse:
+        for p in params:
+            p.main_grad = torch.zeros_like(p, dtype=torch.float32)
+            p.grad_added_to_main_grad = False
+
+    x = (torch.randn(tokens, in_f, device=DEV, dtype=torch.float32) * 0.5).to(torch.bfloat16)
+    x.requires_grad_(True)
+    with fp8_autocast(enabled=True, fp8_recipe=recipe):
+        if module_kind == "linear":
+            out = mod(x)
+        else:
+            out = mod(x, m_splits)
+
+    offs = [0, m_splits[0], tokens] if module_kind == "grouped" else [0, tokens]
+    nsplit = len(params)
+
+    def cat_ref(weights, mat, transpose):
+        return torch.cat(
+            [
+                mat[offs[i] : offs[i + 1]].to(torch.float32)
+                @ (
+                    weights[i].detach().to(torch.float32).t()
+                    if transpose
+                    else weights[i].detach().to(torch.float32)
+                )
+                for i in range(nsplit)
+            ]
+        )
+
+    ref_out_hat = cat_ref(w_hats, x, True)
+    ref_out_master = cat_ref(params, x, True)
+    e_hat, e_master = rel_err(out, ref_out_hat), rel_err(out, ref_out_master)
+    assert e_hat < FP8_TOL_1OP, f"fwd err {e_hat}"
+    assert e_hat < e_master, "fwd closer to master than to the MXFP4-grid weight"
+
+    g = (torch.randn_like(out, dtype=torch.float32) * 0.1).to(out.dtype)
+    out.backward(g)
+
+    ref_dx_hat = cat_ref(w_hats, g, False)
+    ref_dx_master = cat_ref(params, g, False)
+    e_h, e_m = rel_err(x.grad, ref_dx_hat), rel_err(x.grad, ref_dx_master)
+    if override == "high_precision":
+        assert e_m < BF16_TOL, f"hp dgrad err vs master {e_m}"
+        assert e_m < e_h, "hp dgrad must use the unquantized master (B)"
+    elif override == "dequantized":
+        assert e_h < BF16_TOL, f"dequantized dgrad err vs w_hat {e_h}"
+        assert e_h < e_m, "dequantized dgrad must use the MXFP4-grid weight (A)"
+    else:
+        assert e_h < FP8_TOL_1OP, f"mode-None dgrad err {e_h}"
+        assert e_h < e_m, "mode-None dgrad must use the MXFP4-grid weight (C)"
+
+    if module_kind == "linear":
+        with fp8_autocast(enabled=True, fp8_recipe=recipe):
+            out_c1 = mod(x.detach(), is_first_microbatch=True)
+            out_c2 = mod(x.detach(), is_first_microbatch=False)
+        assert rel_err(out_c2, ref_out_hat) < FP8_TOL_1OP, "cached-weight fwd drifted"
+        assert rel_err(out_c1, out_c2) < 1e-6, "workspace update changed the weights"
+
+    wgrad_tol = BF16_TOL if override == "high_precision" else FP8_TOL_2OP
+    for i, p in enumerate(params):
+        ref_wg = g[offs[i] : offs[i + 1]].to(torch.float32).t() @ x[offs[i] : offs[i + 1]].to(
+            torch.float32
+        )
+        if fuse:
+            assert p.grad_added_to_main_grad is True, "fused-wgrad flag not set (G)"
+            e = rel_err(p.main_grad, ref_wg)
+            assert e < wgrad_tol, f"main_grad err {e} (G)"
+        else:
+            assert p.grad is not None, "STE did not deliver a weight grad (G)"
+            e = rel_err(p.grad, ref_wg)
+            assert e < wgrad_tol, f"wgrad err {e} (G)"
+
+
+def test_e2e_matrix():
+    for module_kind in ("linear", "grouped"):
+        for recipe_cls in (MXFP4QATMXFP8BlockScaling, MXFP4QATFloat8BlockScaling):
+            for override in (None, "dequantized", "high_precision"):
+                for fuse in (False, True):
+                    _run_module(module_kind, recipe_cls(), override, fuse)
+                    print(
+                        f"  e2e OK {module_kind} {recipe_cls.__name__} "
+                        f"override={override} fuse={fuse}"
+                    )
+
+
+TESTS = [
+    test_fake_quant_grid_scale_rtne,
+    test_D_lossless_in_bf16,
+    test_E_mxfp8_rowwise_lossless,
+    test_dequantize_e8m0_extreme_codes,
+    test_G_blockwise_lossless_and_transpose,
+    test_C_mxfp8_colwise_bound,
+    test_lossless_dequant_to_bf16,
+    test_kernel_matches_torch_reference,
+    test_edge_value_domains,
+    test_blockwise_recipe_is_128x128,
+    test_blockwise_over_ratio_is_bounded_not_lossless,
+    test_exp2f_e8m0_all_codes,
+    test_ste_identity_gradient,
+    test_misaligned_and_noncontiguous,
+    test_bf16_exhaustive_bit_patterns,
+    test_fp32_bit_fuzz,
+    test_scale_threshold_and_rtne_midpoints,
+    test_deployment_top_cap_domain,
+    test_blockwise_feasibility_enumeration,
+    test_recipe_switch_invalidates_cache,
+    test_recipe_field_controls_qat,
+    test_ops_api_rejected,
+    test_fp16_pipeline_rejected,
+    test_kernel_perf,
+    test_kernel_fast_math_immune,
+    test_e2e_matrix,
+]
+
+if __name__ == "__main__":
+    import sys
+    import traceback
+
+    failed = 0
+    for t in TESTS:
+        try:
+            t()
+            torch.cuda.synchronize()
+            print(f"PASS {t.__name__}")
+        except Exception:
+            failed += 1
+            print(f"FAIL {t.__name__}")
+            traceback.print_exc()
+    print(f"{len(TESTS) - failed}/{len(TESTS)} passed")
+    sys.exit(1 if failed else 0)
diff --git a/transformer_engine/common/cast/cast.cu b/transformer_engine/common/cast/cast.cu
index 1e3c04573b..bf44e3bffb 100644
--- a/transformer_engine/common/cast/cast.cu
+++ b/transformer_engine/common/cast/cast.cu
@@ -16,6 +16,7 @@
 #include "../utils.cuh"
 #include "dispatch/dequantize.cuh"
 #include "dispatch/quantize.cuh"
+#include "mxfp4/fake_quantize_mxfp4.cuh"
 #include "transformer_engine/transpose.h"
 
 void nvte_quantize(const NVTETensor input, NVTETensor output, cudaStream_t stream) {
@@ -47,6 +48,38 @@ void nvte_quantize_v2(const NVTETensor input, NVTETensor output,
   dispatch::quantize_fwd_helper<IS_ACT, Empty, nullptr>(input, output, quant_config, stream);
 }
 
+void nvte_mxfp4_fake_quantize(const NVTETensor input, NVTETensor output, cudaStream_t stream) {
+  NVTE_API_CALL(nvte_mxfp4_fake_quantize);
+  using namespace transformer_engine;
+  const Tensor &input_tensor = *convertNVTETensorCheck(input);
+  Tensor &output_tensor = *convertNVTETensorCheck(output);
+  NVTE_CHECK(input_tensor.dtype() == DType::kBFloat16 || input_tensor.dtype() == DType::kFloat32,
+             "MXFP4 fake-quantize supports BF16/FP32 inputs, got ",
+             to_string(input_tensor.dtype()));
+  NVTE_CHECK(output_tensor.dtype() == input_tensor.dtype(),
+             "MXFP4 fake-quantize output dtype must match the input dtype.");
+  NVTE_CHECK(input_tensor.data.shape == output_tensor.data.shape,
+             "MXFP4 fake-quantize output shape must match the input shape.");
+  const size_t numel = input_tensor.numel();
+  NVTE_CHECK(numel % 32 == 0, "MXFP4 fake-quantize needs the element count divisible by 32.");
+  NVTE_CHECK(!input_tensor.data.shape.empty() && input_tensor.data.shape.back() % 32 == 0,
+             "MXFP4 fake-quantize needs the innermost dimension divisible by 32 "
+             "(1x32 blocks must not cross rows).");
+  if (numel == 0) {
+    return;
+  }
+  if (input_tensor.dtype() == DType::kBFloat16) {
+    dispatch::mxfp4::fake_quantize_mxfp4_launch<bf16>(
+        reinterpret_cast<const bf16 *>(input_tensor.data.dptr),
+        reinterpret_cast<bf16 *>(output_tensor.data.dptr), numel, stream);
+  } else {
+    dispatch::mxfp4::fake_quantize_mxfp4_launch<fp32>(
+        reinterpret_cast<const fp32 *>(input_tensor.data.dptr),
+        reinterpret_cast<fp32 *>(output_tensor.data.dptr), numel, stream);
+  }
+  NVTE_CHECK_CUDA(cudaGetLastError());
+}
+
 void nvte_dequantize(const NVTETensor input, NVTETensor output, cudaStream_t stream) {
   NVTE_API_CALL(nvte_dequantize);
   using namespace transformer_engine;
diff --git a/transformer_engine/common/cast/mxfp4/fake_quantize_mxfp4.cuh b/transformer_engine/common/cast/mxfp4/fake_quantize_mxfp4.cuh
new file mode 100644
index 0000000000..6b4c437221
--- /dev/null
+++ b/transformer_engine/common/cast/mxfp4/fake_quantize_mxfp4.cuh
@@ -0,0 +1,127 @@
+/*************************************************************************
+ * Copyright (c) 2022-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+ *
+ * See LICENSE for license information.
+ ************************************************************************/
+
+#ifndef TRANSFORMER_ENGINE_CAST_MXFP4_FAKE_QUANTIZE_MXFP4_CUH_
+#define TRANSFORMER_ENGINE_CAST_MXFP4_FAKE_QUANTIZE_MXFP4_CUH_
+
+#include <cuda.h>
+#include <cuda_runtime.h>
+
+#include "../../utils.cuh"
+
+namespace transformer_engine {
+namespace dispatch {
+namespace mxfp4 {
+namespace fake_quantize_kernel {
+
+constexpr size_t MXFP4_BLOCK_SIZE = 32;
+constexpr size_t THREADS_PER_CHUNK = 256;
+constexpr float E2M1_MAX = 6.0f;
+constexpr int SCALE_EXP_MIN = -126;
+constexpr int SCALE_EXP_MAX = 125;
+
+__device__ __forceinline__ float round_e2m1_magnitude(const float y) {
+  if (y <= 2.0f) {
+    return rintf(y * 2.0f) * 0.5f;
+  }
+  if (y <= 4.0f) {
+    return rintf(y);
+  }
+  return rintf(y * 0.5f) * 2.0f;
+}
+
+__device__ __forceinline__ int block_scale_exponent_from_bits(const int bits) {
+  const int exp_field = bits >> 23;
+  const int mantissa = bits & 0x7FFFFF;
+  int e = exp_field - 129 + (mantissa > 0x400000 ? 1 : 0);
+  if (exp_field == 0) {
+    e = SCALE_EXP_MIN;
+  }
+  return max(SCALE_EXP_MIN, min(SCALE_EXP_MAX, e));
+}
+
+__device__ __forceinline__ float exp2_from_int(const int e) {
+  return __int_as_float((e + 127) << 23);
+}
+
+__device__ __forceinline__ float mul_rn_no_ftz(const float a, const float b) {
+  float r;
+  asm("mul.rn.f32 %0, %1, %2;" : "=f"(r) : "f"(a), "f"(b));
+  return r;
+}
+
+template <typename IType>
+__global__ void __launch_bounds__(THREADS_PER_CHUNK)
+    fake_quantize_mxfp4_kernel(const IType *__restrict__ input, IType *__restrict__ output,
+                               const size_t num_vectors) {
+  constexpr size_t VEC_ELTS = 16 / sizeof(IType);
+  constexpr int LANES_PER_BLOCK = MXFP4_BLOCK_SIZE / VEC_ELTS;
+
+  const size_t stride = static_cast<size_t>(gridDim.x) * blockDim.x;
+  for (size_t vec_idx = static_cast<size_t>(blockIdx.x) * blockDim.x + threadIdx.x;
+       vec_idx < num_vectors; vec_idx += stride) {
+    Vec<IType, VEC_ELTS> in_vec;
+    in_vec.load_from(input, vec_idx);
+
+    float elts[VEC_ELTS];
+    int amax_bits = 0;
+#pragma unroll
+    for (int i = 0; i < static_cast<int>(VEC_ELTS); ++i) {
+      elts[i] = static_cast<float>(in_vec.data.elt[i]);
+      amax_bits = max(amax_bits, __float_as_int(elts[i]) & 0x7fffffff);
+    }
+    const unsigned lane = threadIdx.x & 31u;
+    const unsigned group_mask = (LANES_PER_BLOCK == 32) ? 0xffffffffu
+                                                        : (((1u << LANES_PER_BLOCK) - 1u)
+                                                           << (lane & ~(LANES_PER_BLOCK - 1u)));
+#pragma unroll
+    for (int offset = LANES_PER_BLOCK / 2; offset > 0; offset >>= 1) {
+      amax_bits = max(amax_bits, __shfl_xor_sync(group_mask, amax_bits, offset));
+    }
+    const int nonfinite = amax_bits >= 0x7f800000;
+
+    Vec<IType, VEC_ELTS> out_vec;
+    if (nonfinite) {
+      const float qnan = __int_as_float(0x7fc00000);
+#pragma unroll
+      for (int i = 0; i < static_cast<int>(VEC_ELTS); ++i) {
+        out_vec.data.elt[i] = static_cast<IType>(qnan);
+      }
+    } else {
+      const int e = (amax_bits > 0) ? block_scale_exponent_from_bits(amax_bits) : 0;
+      const float scale = exp2_from_int(e);
+      const float inv_scale = exp2_from_int(-e);
+#pragma unroll
+      for (int i = 0; i < static_cast<int>(VEC_ELTS); ++i) {
+        const float av = __int_as_float(__float_as_int(elts[i]) & 0x7fffffff);
+        const float y = fminf(mul_rn_no_ftz(av, inv_scale), E2M1_MAX);
+        const float q = copysignf(round_e2m1_magnitude(y), elts[i]);
+        out_vec.data.elt[i] = static_cast<IType>(mul_rn_no_ftz(q, scale));
+      }
+    }
+    out_vec.store_to(output, vec_idx);
+  }
+}
+
+}  // namespace fake_quantize_kernel
+
+template <typename IType>
+void fake_quantize_mxfp4_launch(const IType *input, IType *output, const size_t numel,
+                                cudaStream_t stream) {
+  using namespace fake_quantize_kernel;
+  constexpr size_t VEC_ELTS = 16 / sizeof(IType);
+  const size_t num_vectors = numel / VEC_ELTS;
+  const size_t blocks_needed = (num_vectors + THREADS_PER_CHUNK - 1) / THREADS_PER_CHUNK;
+  const int grid = static_cast<int>(blocks_needed < 65535 ? blocks_needed : 65535);
+  fake_quantize_mxfp4_kernel<IType>
+      <<<grid, THREADS_PER_CHUNK, 0, stream>>>(input, output, num_vectors);
+}
+
+}  // namespace mxfp4
+}  // namespace dispatch
+}  // namespace transformer_engine
+
+#endif
diff --git a/transformer_engine/common/include/transformer_engine/cast.h b/transformer_engine/common/include/transformer_engine/cast.h
index 4d6d24ba65..e19f40d3fd 100644
--- a/transformer_engine/common/include/transformer_engine/cast.h
+++ b/transformer_engine/common/include/transformer_engine/cast.h
@@ -418,6 +418,15 @@ void nvte_group_quantize_dbias_dsrelu(const NVTEGroupedTensor input,
  */
 void nvte_dequantize(const NVTETensor input, NVTETensor output, cudaStream_t stream);
 
+/*! \brief Project a BF16/FP32 tensor onto the MXFP4 (E2M1) grid with 1x32
+ *         power-of-two block scales, writing back dequantized values.
+ *
+ *  \param[in]     input     Input tensor (BF16/FP32, innermost dim divisible by 32).
+ *  \param[in,out] output    Output tensor with the same dtype and shape.
+ *  \param[in]     stream    CUDA stream used for the operation.
+ */
+void nvte_mxfp4_fake_quantize(const NVTETensor input, NVTETensor output, cudaStream_t stream);
+
 /*! \brief Casts input grouped tensor from reduced to higher precision.
  *         In case of the MXFP8 dequantization, the dequantized values are stored to the rowwise
  *         data of the output tensor, regardless of whether the row- or columnwise scaling is used.
diff --git a/transformer_engine/common/recipe/__init__.py b/transformer_engine/common/recipe/__init__.py
index 8c209ace15..a7224f1e63 100644
--- a/transformer_engine/common/recipe/__init__.py
+++ b/transformer_engine/common/recipe/__init__.py
@@ -168,6 +168,12 @@ def custom(cls):
         """Whether the given recipe is custom."""
         return issubclass(cls, CustomRecipe)
 
+    def mxfp4_qat(self):
+        """Whether the recipe applies MXFP4 weight QAT on top of its base recipe."""
+        return bool(getattr(self, "mxfp4_qat_weights", False)) or isinstance(
+            self, (MXFP4QATMXFP8BlockScaling, MXFP4QATFloat8BlockScaling)
+        )
+
 
 @dataclass(repr=False)
 class DelayedScaling(Recipe):
@@ -368,6 +374,7 @@ class MXFP8BlockScaling(Recipe):
     fp8_dpa: bool = False
     fp8_mha: bool = False
     backward_override: Optional[str] = os.getenv("NVTE_BACKWARD_OVERRIDE", None)
+    mxfp4_qat_weights: bool = os.getenv("NVTE_MXFP4_QAT", "0") == "1"
 
     def __post_init__(self) -> None:
         assert self.fp8_format != Format.E5M2, "Pure E5M2 training is not supported."
@@ -437,6 +444,7 @@ class Float8BlockScaling(Recipe):
     fp8_dpa: bool = False
     fp8_mha: bool = False
     backward_override: Optional[str] = os.getenv("NVTE_BACKWARD_OVERRIDE", None)
+    mxfp4_qat_weights: bool = os.getenv("NVTE_MXFP4_QAT", "0") == "1"
 
     def __post_init__(self) -> None:
         assert self.x_block_scaling_dim in [1, 2], "Only 1D or 2D blocks supported for x"
@@ -461,6 +469,17 @@ def __post_init__(self) -> None:
         assert (
             self.backward_override in _BACKWARD_OVERRIDES
         ), "NVTE_BACKWARD_OVERRIDE must be unset or one of: 'high_precision', 'dequantized'."
+        if self.mxfp4_qat_weights:
+            if self.w_block_scaling_dim != 2:
+                raise ValueError(
+                    "MXFP4 QAT requires 128x128 (2D) weight scaling blocks, got "
+                    f"w_block_scaling_dim={self.w_block_scaling_dim}."
+                )
+            if not self.fp8_quant_fwd_weight.power_2_scale:
+                raise ValueError(
+                    "MXFP4 QAT requires power-of-two weight scales; "
+                    "NVTE_FP8_BLOCK_SCALING_FP32_SCALES=1 is incompatible."
+                )
 
     def _make_repr(self) -> str:
         return (
@@ -481,6 +500,34 @@ def _make_repr(self) -> str:
         )
 
 
+@dataclass(repr=False)
+class MXFP4QATMXFP8BlockScaling(MXFP8BlockScaling):
+    """
+    MXFP8 recipe with MXFP4 weight quantization-aware training.
+
+    Weights are projected onto the MXFP4 (E2M1, 1x32 power-of-two scale) grid
+    before the regular MXFP8 weight quantization; the rowwise encoding of the
+    projected weight is lossless. Activations, gradients and
+    ``backward_override`` behave as in the base recipe. bf16/fp32 only.
+    """
+
+    mxfp4_qat_weights: bool = True
+
+
+@dataclass(repr=False)
+class MXFP4QATFloat8BlockScaling(Float8BlockScaling):
+    """
+    Float8 block-scaling recipe with MXFP4 weight quantization-aware training.
+
+    Weights are projected onto the MXFP4 (E2M1, 1x32 power-of-two scale) grid
+    before the regular 128x128 blockwise FP8 weight quantization (lossless
+    within each tile's FP8 dynamic-range headroom). Activations, gradients and
+    ``backward_override`` behave as in the base recipe. bf16/fp32 only.
+    """
+
+    mxfp4_qat_weights: bool = True
+
+
 @dataclass(repr=False)
 class NVFP4BlockScaling(Recipe):
     """
diff --git a/transformer_engine/common/util/ptx.cuh b/transformer_engine/common/util/ptx.cuh
index 2814aa3490..94a82990e3 100644
--- a/transformer_engine/common/util/ptx.cuh
+++ b/transformer_engine/common/util/ptx.cuh
@@ -380,6 +380,8 @@ __device__ __forceinline__ bf16 exp2f_rcp<bf16>(e8m0_t biased_exp) {
 }
 
 __device__ __forceinline__ float exp2f(e8m0_t biased_exp) {
+  if (biased_exp == 0) return __int_as_float(0x00400000);
+  if (biased_exp == 255) return __int_as_float(0x7fffffff);
   return __int_as_float(biased_exp << FP32_MANTISSA_BITS);
 }
 
diff --git a/transformer_engine/pytorch/csrc/extensions.h b/transformer_engine/pytorch/csrc/extensions.h
index 6edfbdc00e..2a77895ead 100644
--- a/transformer_engine/pytorch/csrc/extensions.h
+++ b/transformer_engine/pytorch/csrc/extensions.h
@@ -347,6 +347,8 @@ py::object nvfp4_quantize_with_amax(const at::Tensor &tensor, py::handle quantiz
 
 py::object dequantize(const py::handle &input, DType otype);
 
+at::Tensor mxfp4_fake_quantize(const at::Tensor &input);
+
 py::object group_quantize(const at::Tensor &tensor, py::handle quantizer, const size_t num_tensors,
                           std::optional<at::Tensor> first_dims, std::optional<at::Tensor> last_dims,
                           std::optional<at::Tensor> tensor_offsets,
diff --git a/transformer_engine/pytorch/csrc/extensions/cast.cpp b/transformer_engine/pytorch/csrc/extensions/cast.cpp
index 8b1cd384aa..41dcd03e7f 100644
--- a/transformer_engine/pytorch/csrc/extensions/cast.cpp
+++ b/transformer_engine/pytorch/csrc/extensions/cast.cpp
@@ -556,6 +556,24 @@ py::object bgrad_group_quantize(const at::Tensor &tensor, py::handle quantizer,
                         py::cast(std::move(dbias_torch)));
 }
 
+at::Tensor mxfp4_fake_quantize(const at::Tensor &input) {
+  TORCH_CHECK(input.is_cuda(), "mxfp4_fake_quantize: input must be a CUDA tensor.");
+  TORCH_CHECK(input.dim() == 2, "mxfp4_fake_quantize: input must be 2D, got dim=", input.dim());
+  TORCH_CHECK(input.size(1) % 32 == 0,
+              "mxfp4_fake_quantize: innermost dimension must be divisible by 32, got ",
+              input.size(1));
+  const c10::cuda::CUDAGuard device_guard(input.device());
+  auto input_contiguous = input.contiguous();
+  if (reinterpret_cast<uintptr_t>(input_contiguous.data_ptr()) % 16 != 0) {
+    input_contiguous = input_contiguous.clone();
+  }
+  auto output = at::empty_like(input_contiguous);
+  auto input_cpp = makeTransformerEngineTensor(input_contiguous);
+  auto output_cpp = makeTransformerEngineTensor(output);
+  nvte_mxfp4_fake_quantize(input_cpp.data(), output_cpp.data(), at::cuda::getCurrentCUDAStream());
+  return output;
+}
+
 py::object dequantize(const py::handle &input, transformer_engine::DType otype) {
   init_extension();
 
diff --git a/transformer_engine/pytorch/csrc/extensions/pybind.cpp b/transformer_engine/pytorch/csrc/extensions/pybind.cpp
index 7e9d114be8..facbb0c489 100644
--- a/transformer_engine/pytorch/csrc/extensions/pybind.cpp
+++ b/transformer_engine/pytorch/csrc/extensions/pybind.cpp
@@ -199,6 +199,9 @@ PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
 
   m.def("quantize", transformer_engine::pytorch::quantize, py::arg("tensor"), py::arg("quantizer"),
         py::arg("output") = py::none(), py::arg("noop") = py::none());
+  m.def("mxfp4_fake_quantize", &transformer_engine::pytorch::mxfp4_fake_quantize,
+        "Project a BF16/FP32 tensor onto the MXFP4 grid (QAT fake-quantization)", py::arg("input"),
+        py::call_guard<py::gil_scoped_release>());
   m.def("dequantize", &transformer_engine::pytorch::dequantize, "Dequantize", py::arg("input"),
         py::arg("otype"));
   m.def("create_empty_quantized_tensor",
diff --git a/transformer_engine/pytorch/module/base.py b/transformer_engine/pytorch/module/base.py
index 7c11e635c6..5447aa676f 100644
--- a/transformer_engine/pytorch/module/base.py
+++ b/transformer_engine/pytorch/module/base.py
@@ -46,6 +46,7 @@
     _fsdp_gather_tensors,
 )
 from ..constants import dist_group_type
+from ..mxfp4_qat import mxfp4_fake_quantize
 from ..cpp_extensions.gemm import _NUM_MAX_UB_STREAMS
 from ..quantized_tensor import QuantizedTensor, QuantizedTensorStorage, Quantizer
 from ..tensor.float8_tensor import Float8Quantizer, Float8CurrentScalingQuantizer
@@ -800,8 +801,26 @@ def quantize_weight(
         ``_fp8_workspaces``.
     """
 
+    _mxfp4_qat_active = (
+        FP8GlobalStateManager.get_fp8_recipe().mxfp4_qat()
+        if FP8GlobalStateManager.is_fp8_enabled()
+        else False
+    )
+    if _mxfp4_qat_active and workspace_dtype == torch.float16:
+        raise NotImplementedError(
+            "MXFP4 QAT does not support fp16 as the activation/dequantize dtype: "
+            "the MXFP4 grid (values up to 6*2^125) exceeds fp16 range. Use bf16 "
+            "or fp32."
+        )
+
     # Already-quantized weight (primary FP8 parameters)
     if isinstance(tensor, QuantizedTensor):
+        if _mxfp4_qat_active:
+            raise NotImplementedError(
+                "MXFP4 QAT recipes do not support primary quantized weights: the "
+                "high-precision master weight is required to project onto the "
+                "MXFP4 grid."
+            )
         update_rowwise = True if quantizer.rowwise_usage else None
         update_columnwise = True if quantizer.columnwise_usage else None
         tensor.update_usage(
@@ -833,6 +852,8 @@ def quantize_weight(
         if update_workspace:
             if tensor is None:
                 raise ValueError("tensor kwarg must be provided to update FP8 workspace")
+            if _mxfp4_qat_active:
+                tensor = mxfp4_fake_quantize(tensor)
             if hasattr(workspace, "quantize_"):
                 workspace.quantize_(tensor, noop_flag=skip_update_flag)
             else:
@@ -842,6 +863,8 @@ def quantize_weight(
     # Cache miss — create new workspace
     if tensor is None or quantizer is None:
         raise ValueError("tensor and quantizer kwargs must be provided to construct FP8 workspace")
+    if _mxfp4_qat_active:
+        tensor = mxfp4_fake_quantize(tensor)
     if cache:
         # Ensure the tensor in the cache is an instance of torch.Tensor,
         # as it persists beyond a single forward pass.
@@ -1506,9 +1529,12 @@ def init_fp8_metadata(self, num_gemms: int = 1) -> None:
             meta["recipe"] = FP8GlobalStateManager.get_fp8_recipe()
 
         _current_recipe = meta["recipe"]
-        if _original_recipe is not None and not (
-            issubclass(_current_recipe.__class__, _original_recipe.__class__)
-            or issubclass(_original_recipe.__class__, _current_recipe.__class__)
+        if _original_recipe is not None and (
+            not (
+                issubclass(_current_recipe.__class__, _original_recipe.__class__)
+                or issubclass(_original_recipe.__class__, _current_recipe.__class__)
+            )
+            or _current_recipe.mxfp4_qat() != _original_recipe.mxfp4_qat()
         ):
             warnings.warn(
                 f"Recipe type changed from {_original_recipe.__class__.__name__} "
diff --git a/transformer_engine/pytorch/module/layernorm_linear.py b/transformer_engine/pytorch/module/layernorm_linear.py
index a588e21a0c..d17e221325 100644
--- a/transformer_engine/pytorch/module/layernorm_linear.py
+++ b/transformer_engine/pytorch/module/layernorm_linear.py
@@ -32,6 +32,7 @@
     _2X_ACC_WGRAD,
 )
 from ..quantization import FP8GlobalStateManager, QuantizerRole
+from ..mxfp4_qat import mxfp4_fake_quantize
 from ..utils import (
     assert_dim_for_fp8_exec,
     cast_if_needed,
@@ -810,6 +811,8 @@ def backward(
                     # fsdp2 quantized-tensor hooks when workspace was not saved.
                     weight = saved_weight
                 elif ctx.weight_quantizer is not None:
+                    if ctx.fp8_recipe is not None and ctx.fp8_recipe.mxfp4_qat():
+                        saved_weight = mxfp4_fake_quantize(saved_weight)
                     ctx.weight_quantizer.set_usage(rowwise=True, columnwise=True)
                     weight = ctx.weight_quantizer(saved_weight)
 
diff --git a/transformer_engine/pytorch/module/layernorm_mlp.py b/transformer_engine/pytorch/module/layernorm_mlp.py
index 19e20d775c..dfa9f7f0f8 100644
--- a/transformer_engine/pytorch/module/layernorm_mlp.py
+++ b/transformer_engine/pytorch/module/layernorm_mlp.py
@@ -33,6 +33,7 @@
     _2X_ACC_WGRAD,
 )
 from ..quantization import FP8GlobalStateManager, QuantizerRole
+from ..mxfp4_qat import mxfp4_fake_quantize
 from ..jit import (
     bias_gelu_fused,
     bgrad_dgelu_fused,
@@ -247,11 +248,12 @@ def _forward(
             backward_override = FP8GlobalStateManager.get_fp8_recipe().backward_override
         else:
             backward_override = None
-        assert backward_override is None, (
-            "NVTE_BACKWARD_OVERRIDE=high_precision/dequantized is not implemented in LayerNormMLP."
-            " Replace LayerNormMLP with LayerNormLinear + Linear to enable"
-            " high_precision/dequantized backward."
-        )
+        if backward_override is not None:
+            raise NotImplementedError(
+                "NVTE_BACKWARD_OVERRIDE=high_precision/dequantized is not implemented in"
+                " LayerNormMLP. Replace LayerNormMLP with LayerNormLinear + Linear to enable"
+                " high_precision/dequantized backward."
+            )
 
         # if grad is enabled and this is not the bwd stage, we must save this so bwd knows which path to take
         if is_grad_enabled and not recompute_for_bwd:
@@ -1228,6 +1230,8 @@ def backward(
                 if isinstance(origin_fc2_weight, QuantizedTensorStorage):
                     fc2_weight = origin_fc2_weight
                 elif ctx.fc2_weight_quantizer is not None:
+                    if ctx.fp8_recipe is not None and ctx.fp8_recipe.mxfp4_qat():
+                        origin_fc2_weight = mxfp4_fake_quantize(origin_fc2_weight)
                     ctx.fc2_weight_quantizer.set_usage(rowwise=True, columnwise=True)
                     fc2_weight = ctx.fc2_weight_quantizer(origin_fc2_weight)
 
@@ -1522,6 +1526,8 @@ def fc2_wgrad_gemm(
                 if isinstance(origin_fc1_weight, QuantizedTensorStorage):
                     fc1_weight = origin_fc1_weight
                 elif ctx.fc1_weight_quantizer is not None:
+                    if ctx.fp8_recipe is not None and ctx.fp8_recipe.mxfp4_qat():
+                        origin_fc1_weight = mxfp4_fake_quantize(origin_fc1_weight)
                     ctx.fc1_weight_quantizer.set_usage(rowwise=True, columnwise=True)
                     fc1_weight = ctx.fc1_weight_quantizer(origin_fc1_weight)
 
diff --git a/transformer_engine/pytorch/module/linear.py b/transformer_engine/pytorch/module/linear.py
index 6b0f941bc4..0fe378b11a 100644
--- a/transformer_engine/pytorch/module/linear.py
+++ b/transformer_engine/pytorch/module/linear.py
@@ -15,6 +15,7 @@
 import transformer_engine_torch as tex
 
 from transformer_engine.common.recipe import Recipe
+from ..mxfp4_qat import mxfp4_fake_quantize
 from transformer_engine.pytorch.torch_version import torch_version
 
 from .base import (
@@ -1003,6 +1004,8 @@ def _linear_backward(args: LinearBwdArgs) -> Tuple[Union[torch.Tensor, None], ..
                     # fsdp2 quantized-tensor hooks when workspace was not saved.
                     weight_fp8 = saved_weight
                 elif bwd_args.weight_quantizer is not None:
+                    if bwd_args.fp8 and FP8GlobalStateManager.get_fp8_recipe().mxfp4_qat():
+                        saved_weight = mxfp4_fake_quantize(saved_weight)
                     bwd_args.weight_quantizer.set_usage(rowwise=True, columnwise=True)
                     weight_fp8 = bwd_args.weight_quantizer(saved_weight)
             elif (
@@ -1013,6 +1016,8 @@ def _linear_backward(args: LinearBwdArgs) -> Tuple[Union[torch.Tensor, None], ..
             ):
                 # Distributed weight re-gathered a BF16 weight: quantize with the layer quantizer
                 # so the dgrad operand isn't cast by the delayed recipe.
+                if FP8GlobalStateManager.get_fp8_recipe().mxfp4_qat():
+                    weight_fp8 = mxfp4_fake_quantize(weight_fp8)
                 bwd_args.weight_quantizer.set_usage(rowwise=True, columnwise=True)
                 weight_fp8 = bwd_args.weight_quantizer(weight_fp8)
 
diff --git a/transformer_engine/pytorch/mxfp4_qat.py b/transformer_engine/pytorch/mxfp4_qat.py
new file mode 100644
index 0000000000..57ee4f854c
--- /dev/null
+++ b/transformer_engine/pytorch/mxfp4_qat.py
@@ -0,0 +1,71 @@
+# Copyright (c) 2022-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+#
+# See LICENSE for license information.
+
+"""MXFP4 weight fake-quantization for quantization-aware training.
+
+Projects weights onto the MXFP4 (E2M1) grid with 1x32 power-of-two block
+scales and returns the dequantized values in the input dtype. The output is
+exact in bf16/fp32, and the host recipe's weight quantization (MXFP8 rowwise;
+128x128 blockwise within the tile scale-spread bound) is DECODED-VALUE exact
+on this finite QAT domain. It is not raw-tuple identical to a direct
+MXFP4->MXFP8 converter: the host encoder re-canonicalizes (payload, scale)
+per block, so equal values may carry different E4M3/UE8M0 bytes, and the
+original MXFP4 scale metadata is not preserved. Scale floor 2^-126 / cap 2^125
+(deployment contract); non-finite values NaN-poison their 1x32 block; fp16 is
+rejected. The numerical specification lives in
+tests/pytorch/references/mxfp4_qat_reference.py.
+"""
+import torch
+
+__all__ = ["mxfp4_fake_quantize"]
+
+_MXFP4_BLOCK = 32
+
+
+class _MXFP4FakeQuantizeSTE(torch.autograd.Function):
+    """Straight-through estimator: identity gradient through the projection."""
+
+    @staticmethod
+    def forward(ctx, weight):  # pylint: disable=arguments-differ
+        import transformer_engine_torch as tex
+
+        if not hasattr(tex, "mxfp4_fake_quantize"):
+            raise RuntimeError(
+                "transformer_engine_torch was built without mxfp4_fake_quantize; "
+                "rebuild Transformer Engine."
+            )
+        w = weight.contiguous()
+        if w.data_ptr() % 16 != 0:
+            # logically-contiguous storage-offset views can be misaligned for
+            # the kernel's 16-byte vector accesses
+            w = w.clone()
+        return tex.mxfp4_fake_quantize(w)
+
+    @staticmethod
+    def backward(ctx, grad_output):  # pylint: disable=arguments-differ
+        return grad_output
+
+
+def mxfp4_fake_quantize(weight: torch.Tensor) -> torch.Tensor:
+    """Project ``weight`` onto the MXFP4 grid and return it in the input dtype.
+
+    Per 1x32 block: scale = 2^clamp(ceil(log2(amax / 6)), -126, 125), RTNE onto
+    E2M1, saturate at 6, multiply back. The gradient is the straight-through
+    estimator (identity).
+    """
+    if not weight.is_cuda:
+        raise ValueError("MXFP4 QAT expects a CUDA weight tensor")
+    if weight.dim() != 2:
+        raise ValueError(f"MXFP4 QAT expects a 2D weight, got {tuple(weight.shape)}")
+    if weight.dtype not in (torch.bfloat16, torch.float32):
+        raise ValueError(
+            f"MXFP4 QAT supports bf16/fp32 weights, got {weight.dtype} "
+            "(fp16 cannot represent the full MXFP4 scale range)"
+        )
+    if weight.shape[1] % _MXFP4_BLOCK != 0:
+        raise ValueError(
+            "MXFP4 QAT needs the weight inner dim divisible by "
+            f"{_MXFP4_BLOCK}, got {weight.shape[1]}"
+        )
+    return _MXFP4FakeQuantizeSTE.apply(weight)
diff --git a/transformer_engine/pytorch/ops/basic/basic_linear.py b/transformer_engine/pytorch/ops/basic/basic_linear.py
index cb429055a4..2ceec72b5e 100644
--- a/transformer_engine/pytorch/ops/basic/basic_linear.py
+++ b/transformer_engine/pytorch/ops/basic/basic_linear.py
@@ -558,6 +558,16 @@ def _functional_forward(
 
         # Check weight tensor
         w = weight
+        if (
+            with_quantized_compute
+            and FP8GlobalStateManager.is_fp8_enabled()
+            and FP8GlobalStateManager.get_fp8_recipe().mxfp4_qat()
+        ):
+            raise NotImplementedError(
+                "MXFP4 QAT recipes are not implemented for the operations API "
+                "(te.ops): weight quantization here bypasses the QAT projection. "
+                "Use the module API (te.Linear/te.GroupedLinear/...)."
+            )
         if not with_quantized_compute:
             w = maybe_dequantize(w, dtype)
         elif with_quantized_compute and not is_quantized_tensor(w):
diff --git a/transformer_engine/pytorch/ops/basic/grouped_linear.py b/transformer_engine/pytorch/ops/basic/grouped_linear.py
index d3e2dc0a13..01f245fdec 100644
--- a/transformer_engine/pytorch/ops/basic/grouped_linear.py
+++ b/transformer_engine/pytorch/ops/basic/grouped_linear.py
@@ -488,6 +488,16 @@ def _quantize_weights(
     ) -> Sequence[torch.Tensor]:
         """Construct quantized weight tensors."""
 
+        if (
+            FP8GlobalStateManager.is_fp8_enabled()
+            and FP8GlobalStateManager.get_fp8_recipe().mxfp4_qat()
+        ):
+            raise NotImplementedError(
+                "MXFP4 QAT recipes do not support quantized primary weights in the "
+                "operations API: the high-precision master is required for the QAT "
+                "projection."
+            )
+
         # Manually construct MXFP8 weights
         if isinstance(quantizers[0], MXFP8Quantizer):
             return self._quantize_weights_mxfp8(weights, quantizers)
diff --git a/transformer_engine/pytorch/ops/fused/grouped_mlp.py b/transformer_engine/pytorch/ops/fused/grouped_mlp.py
index 83954a9b3d..0ed4f0fb81 100644
--- a/transformer_engine/pytorch/ops/fused/grouped_mlp.py
+++ b/transformer_engine/pytorch/ops/fused/grouped_mlp.py
@@ -26,6 +26,7 @@
     finalize_weight_grads,
 )
 from ...module.base import _2X_ACC_WGRAD
+from ...fp8 import FP8GlobalStateManager
 from ...quantization import Recipe
 from ...tensor import NVFP4Quantizer, NVFP4Tensor, NVFP4TensorStorage, Quantizer
 from ...tensor.grouped_tensor import GroupedTensor
@@ -960,6 +961,14 @@ def fuser_forward(
         fc2_input_quantizer = fc2_op.get_quantizer("forward", 0)
         fc2_weight_quantizer = fc2_op.get_quantizer("forward", 1)
         fc2_grad_output_quantizer = fc2_op.get_quantizer("backward", 0)
+        if (
+            FP8GlobalStateManager.is_fp8_enabled()
+            and FP8GlobalStateManager.get_fp8_recipe().mxfp4_qat()
+        ):
+            raise NotImplementedError(
+                "MXFP4 QAT recipes are not implemented for fused grouped-MLP ops: "
+                "weight quantization here bypasses the QAT projection."
+            )
 
         # Extract split sizes from extra input
         fc1_split_sizes = basic_op_extra_inputs[0][0]
diff --git a/transformer_engine/pytorch/ops/fused/userbuffers_forward_linear.py b/transformer_engine/pytorch/ops/fused/userbuffers_forward_linear.py
index 3a8ff5438d..ee03a63f06 100644
--- a/transformer_engine/pytorch/ops/fused/userbuffers_forward_linear.py
+++ b/transformer_engine/pytorch/ops/fused/userbuffers_forward_linear.py
@@ -319,6 +319,10 @@ def fuser_forward(
                 raise RuntimeError(
                     f"Unsupported recipe for Userbuffers ({recipe.__class__.__name__})"
                 )
+            if recipe.mxfp4_qat():
+                raise NotImplementedError(
+                    "MXFP4 QAT recipes are not implemented for Userbuffers fused ops."
+                )
 
         # Get autocast dtype if needed
         if torch.is_autocast_enabled():