Graphics processing unit (GPU) generations scale matrix, special-function, and memory pipelines at different rates, so kernel bottlenecks move as hardware evolves. FlashAttention-4 exposed this imbalance inside attention on NVIDIA Blackwell. We test whether short polynomial programs can accelerate other special-function-unit (SFU) operations in large language models (LLMs). We first compare native PyTorch evaluation with packed fused multiply--add (FMA) programs in an isolated IEEE binary16 (FP16) sweep spanning L2-resident and high-bandwidth-memory (HBM)-resident working sets. We then replace native sigmoid, tanh, and sigmoid linear unit (SiLU) with degree-3 or degree-4 bfloat16 (BF16) programs in four GB200 integration tasks: dense SiLU, tanh-softcapped attention, sigmoid attention, and routed-expert Swish-gated linear unit (SwiGLU). The programs combine analytical symmetry, target-format rounding, and packed arithmetic inside consuming kernels. The isolated paths improve by 1.19--2.19x in L2 and 1.00--1.70x in HBM. The dense-SiLU, tanh-softcapped-attention, and routed-expert substitutions improve complete training-step throughput by 2.7%, 2.9%, and 8.0%, respectively. The sigmoid-attention substitution improves complete-attention forward by 7.4% and the complete GPU step by 0.3%. Same-checkpoint open-weight ablations and one paired pre-training comparison per task extend the evaluation to model behavior. At common horizons near 100 billion tokens, the final smoothed training-loss differences (polynomial minus native) range from −0.107 to +0.079 across the four tasks.
Figures & tables
Case
Baseline
Polynomial intervention
Question
B1
Native SiLU in dense MLPs
Packed D3 SiLU
Does packed SiLU increase dense-model training-step throughput?
B2
Tanh-softcapped attention
D4 tanh inside FA4
Does packed tanh reduce complete softcapped-attention runtime?
B3
Native sigmoid attention
Direct D3 sigmoid and D4 derivative
Does polynomial sigmoid reduce complete sigmoid-attention runtime?
B4
Native SiLU in routed experts
Fused D3 polynomial Swish-gated linear unit (SwiGLU)
Does packed SwiGLU increase routed-expert training-step throughput?
Table 1: The four principal interventions. D3 and D4 denote degree-3 and degree-4 polynomials. B3 replaces sigmoid and its derivative; softmax exponential is outside its scope.
Figure 1: Analytic reduction before fitting. The lightly shaded positive half-domain is the only side that needs a polynomial. (a) Centering sigmoid exposes an odd residual, so the negative half is recovered by a sign-bit operation. (b) Subtracting x/2 from SiLU exposes an even residual despite SiLU itself being asymmetric. (c) Centering the SiLU derivative exposes an odd residual and enforces s′(x)+s′(−x)=1 by construction.
Figure 2: Constrained fitting and target-program replay, illustrated with the selected D4 centered-sigmoid candidate. (a) In normalized coordinates, every free basis term uj−u vanishes at both endpoints, so optimization cannot violate the two fixed values. (b) Fit interval Li and runtime clamp Lc are searched independently; gray cells violate the Lc≤Li+1 search constraint and the star marks the selected (5.0,5.28125) configuration. (c) Absolute-error curves separate the FP64 fit, BF16 coefficient rounding, the fitter’s split multiply/add replay, and a packed-FMA-style replay that rounds once per multiply-add. The execution schedule changes both error magnitude and where the worst error occurs.
Maximum absolute error
Direction
Function
Reduced program
Endpoint-constrained least squares
Sollya minimax
Forward
Sigmoid
D3, 21+ odd
5.44
5.47
Forward
Tanh
D4, odd
9.61
7.22
Forward
SiLU
D3, x/2+ even
21.74
19.95
Forward
GELU
D5, x/2+ even
1.06
0.50
Backward
Sigmoid
D4, even
4.46
1.97
Table 2: Maximum absolute error ( ×10−3 ) for representative constrained fits and same-form Sollya fpminimax controls. Both use BF16-rounded coefficients under the same real-arithmetic replay; intermediate target-precision rounding is evaluated in appendix A .
Direction
Function
Polynomial program
4M elements (L2 cache)
268M elements (HBM)
Forward
Sigmoid
D3 polynomial
1.702x
1.249x
Forward
Tanh
D4 polynomial
1.187x
1.143x
Forward
SiLU
D3 polynomial
1.536x
1.267x
Forward
GELU
D5 polynomial
1.673x
1.702x
Backward
Sigmoid
D4 direct derivative
2.135x
1.004x
Backward
Tanh
D4 direct derivative
2.184x
1.003x
Table 3: Isolated FP16 function speed on GB200. Speedup is native PyTorch time divided by packed polynomial CUDA time, rounded to three decimals. The two columns report 4,194,304 elements (L2-cache resident) and 268,435,456 elements (HBM resident). Backward polynomial rows are direct derivative fits evaluated from x ; native sigmoid and tanh backward use precomputed forward outputs.
Case
Timed component
Program
Forward speedup
Backward speedup
B1
One-kernel SiLU activation
D3 direct-fit probe
1.055x
1.007x
B2
Complete causal GQA attention
D4 tanh softcap
2.053x
1.046x
B3
Complete causal GQA attention
D3/D4 sigmoid
1.074x
1.001x
B4
One-kernel expert activation
D3 SwiGLU
1.003x
1.005x
Table 4: B1–B4 component performance. Speedup is native time divided by polynomial time; ratios are rounded to three decimals. Each ratio uses the medians of five alternating, 100-iteration round means. B1 evaluates 58,720,256 values; B2 and B3 use causal grouped-query attention (GQA) with batch 1, sequence length 4096, 32 query heads, 8 key/value heads, and head dimension 128; B4 evaluates 24,576 rows with hidden size 1280. The B1 row is a unary direct-derivative probe; the complete-step B1 measurements use the deployed packed SwiGLU path described below.
Case
Model / intervention
Forward speedup
Backward speedup
GPU-step speedup
Training-throughput ratio
B1
Llama 3 8B, dense D3 SiLU
1.005x
1.038x
1.027x
1.027x
B2
Llama 3 8B, D4 tanh softcap
1.107x
1.004x
1.029x
1.029x
B3
Llama 3 8B, D3/D4 sigmoid attention
1.011x
0.997x
1.003x
1.002x
B4
DeepSeek-V3 27A4B, expert D3 SwiGLU
1.087x
1.077x
1.080x
1.080x
Table 5: Complete-model performance. B1–B3 use matched 100-step Llama 3 8B probes on one GB200. B4 uses the complete DeepSeek-V3 27A4B architecture on four GB200s. Timing speedup is native time divided by polynomial time; training-throughput ratio is polynomial tokens/s divided by native tokens/s. Ratios use medians over steps 20–100, are rounded to three decimals, and favor the polynomial above one.
Model
Operator
Fit method
Mean accuracy change (percentage points)
WikiText-2 word perplexity
Prefill throughput ratio
Decode throughput ratio
Qwen2.5 7B
Dense SwiGLU
Endpoint LS
-0.51
8.734→197.677
1.035x
1.013x
Sollya minimax
-39.80
8.734→298,700.735
1.036x
1.015x
Qwen3 30B-A3B
Expert SwiGLU
Endpoint LS
-0.10
7.648→7.655
1.056x
1.025x
Sollya minimax
+0.04
7.648→7.651
1.051x
1.025x
GLM-4.7 Flash
Expert SwiGLU + router sigmoid
Endpoint LS
-0.30
11.406→11.400
0.969x
1.004x
Sollya minimax
-0.18
11.406→11.614
0.966x
1.003x
Table 6: Same-checkpoint degree-4 ablations. “Endpoint LS” is this paper’s endpoint-constrained least-squares fit; “Sollya minimax” is the same-form control. Mean accuracy is the unweighted mean of seven tasks; its change is variant minus native and is rounded to 0.01 percentage points. WikiText-2 word perplexity is native → variant and is rounded to three decimals. Throughput ratios are variant divided by native, are rounded to three decimals, and favor the polynomial above one.
Case
Intervention
Shared horizon (billion tokens)
Native EWM training loss
Polynomial EWM training loss
Training-loss change
B1
Dense D3 SiLU
100.2701
2.6264
2.7050
+0.0786
B2
D4 tanh softcap
100.2701
2.5415
2.5439
+0.0025
B3
D3/D4 sigmoid attention
100.2701
4.9471
4.8397
−0.1073
B4
Routed-expert D3 SwiGLU
100.2439
2.1532
2.1560
+0.0028
Table 7: Paired pre-training at the common observed horizon. Loss is the exponentially weighted mean (EWM) of sampled training loss with a half-life of one billion tokens. Training-loss change is polynomial minus native. Losses and changes are rounded to four decimals. Each row is one paired trajectory.
Figure 3: Paired B1–B4 pre-training loss. Gray is the native/SFU control and coral is the polynomial candidate. Thin lines are uniformly thinned sampled training loss; thick lines are token-domain exponentially weighted means (EWMs) with a half-life of one billion tokens, interpolated after smoothing onto a shared grid with 50 million tokens between points. Each panel ends at the common observed horizon. Insets show the final 25 billion tokens.
Appendix figures & tables5 assets
Supplementary material from the paper’s appendix.
Appendix
Figure 4: Forward numerical accuracy for D3–D6 BF16 target-precision programs. Absolute error is measured against the FP32 reference at BF16-representable inputs.
Figure 5: Backward numerical accuracy for D3–D6 BF16 target-precision programs. Absolute error is measured against the FP32 reference at BF16-representable inputs.
Family
Objective
Fractional core
Maximum relative error
Root-mean-square relative error
D1
Endpoint constrained
1 packed FMA
6.148×10−2
4.466×10−2
PWL2
Least squares
2 packed FMAs + hinge
1.526×10−2
6.754×10−3
PWL2
Sampled minimax
2 packed FMAs + hinge
1.074×10−2
7.174×10−3
D2
Least squares
2 packed FMAs
3.204×10−3
1.987×10−3
D2
Sampled minimax
2 packed FMAs
2.679×10−3
1.921×10−3
FA4 D3
Sollya/minimax
3 packed FMAs; p(1)=1.99983
8.763×10−5
6.029×10−5
Appendix
Table 8: Fractional 2f fits on f∈[0,1] . Errors are measured in FP64 on a dense grid before coefficient quantization. “Core” excludes common range reduction and exponent reconstruction. D1, PWL2, and D2 preserve both endpoints exactly.
Table 9: The exp2 performance ladder on GB200. Entries are speedups over the matching native/all-SFU path. Full-exp rows include unchanged FP32 reduction and reconstruction; FA4 rows use identical causal routing. Values are rounded to three decimals, except the 8B optimizer-step probes, which are rounded to four decimals.
Measured quantity
FA4 D3
PWL2 FP16
D2 FP32
Maximum over output and input-gradient tensors
8.641×10−4
3.458×10−3
2.019×10−3
Appendix
Table 10: Maximum tensor-relative L2 error against all-SFU FA4 on GB200.
Blackwell's 4-bit floating-point (FP4) tensor cores do not automatically make attention faster because softmax conversion and on-chip dependencies dominate once its matrix products shrink. We address this with Direct-P for noncausal inference and a causal path that passes the forward quantization directly into backward. Direct-P maps scores directly to FP4 probabilities and reaches up to 2.13× the bfloat16 (BF16) forward throughput on an NVIDIA GB200. The causal path reconstructs probabilities from saved quantized queries and keys and uses 8-bit floating-point (FP8) gradient operands, accelerating a complete single-GPU 8-billion-parameter update by up to 1.14×. Matched distributed training retains FP8 probabilities and values; every tested MXFP4 probability/value training trajectory diverges.
While most attention logits can be computed in low precision without degrading numerical stability, current attention kernels fail to exploit this phenomenon. We introduce a novel hardware-algorithm co-design in the form of mixed-precision FlashAttention. Our method accumulates key-query products and evaluates their exponentials in 8-bit formats, then adaptively identifies sensitive sub-blocks and recomputes them in 16-bit formats. We propose the specifications for a dedicated accelerator capable of executing this pipeline efficiently. Simulated experiments with Qwen3 and Gemma 3 show that rerouting a selective minority of sub-blocks to high precision is sufficient to recover the baseline model performance.
Stanislav Budzinskiy, Marian Gloser, Tolunay Yilmaz +5
Faculty of Mathematics University of Vienna, Austria · Huawei Heisenberg Research Center, Munich, Germany · Huawei Technologies Co. Ltd
SwiGLU is the dominant MLP activation in modern large language models, yet its intermediate tensor materialization costs 9-37% of MLP execution time. We present two complementary CUTLASS-based SM90 kernels that fuse SwiGLU into GeMM at the tile level. Kernel-1 overlaps Swish computation on the Gate accumulator with Up-tile loading using the Pingpong warp-specialized schedule; Kernel-2 interleaves SwiGLU with tile stores via a custom Epilogue Visitor Tree. Evaluated on Qwen-2.5 models (0.5B-72B) on NVIDIA H100, our kernels achieve up to 2.47x speedup over PyTorch, shifting workloads from memory-bound to compute-bound and reaching 79.5% peak BF16 utilization. We demonstrate that torch.compile cannot replicate this fusion (3-7x slower than our kernels), validating the need for hand-crafted tile-level design. Our fused kernels are also numerically superior, achieving zero mismatches compared to 4.5-11% for cuBLAS.
Abhinav Jangda, Tyler Sorensen, Sebastian Burckhardt +3