DEV Community

Aleksei Romanov
Aleksei Romanov

Posted on

GPU Kernels, Explained: CUDA, Triton, cuBLAS, torch.compile and KernelBench

When you run a neural network in PyTorch, Python is only the coordinator. Behind every tensor operation — from a simple activation to a multi-head attention block — sits one or more GPU kernels: specialized functions executed by tens of thousands of hardware threads in parallel. Heavy matrix multiplications call into NVIDIA’s hand-tuned cuBLAS library, convolutions route to cuDNN, and custom elementwise ops run either as kernels built into PyTorch or as JIT-compiled code generated by torch.compile. In practice, whether your training loop takes one second per step or three seconds is decided almost entirely inside these kernels.

In this deep dive, we look at how GPU kernels actually work under the hood. We will build the same ReLU kernel in raw CUDA and in Triton, use the classic roofline model to see why most deep-learning operations are choked by memory bandwidth rather than arithmetic, and walk through four core optimizations: operator fusion, streaming online softmax, tiled matrix multiplies with Tensor Cores, and GEMM epilogues. Finally, we compare the four ways kernels reach hardware in modern PyTorch pipelines and examine KernelBench — the benchmark testing whether LLMs can write production kernels that beat PyTorch.

Note on numbers: Latency figures in this post are theoretical lower bounds derived from NVIDIA H100 SXM datasheet specifications to illustrate hardware limits, rather than live empirical benchmarks.

Start here

What a GPU kernel is and how it maps work onto threads; the same kernel in CUDA and Triton; the roofline model and why fusion, tiling and Tensor Cores matter; what cuBLAS, cuDNN, eager PyTorch and torch.compile each do; how KernelBench scores a kernel and how to compare its numbers.

  • Kernel: A function that runs on the GPU. The CPU launches it once; the GPU runs many copies of it in parallel, each on a different part of the data.
  • Thread, block, grid: A thread runs one copy of the kernel. Threads are grouped into blocks that share fast on-chip memory; all blocks of one launch form the grid.
  • SM: Streaming Multiprocessor: one of the GPU's processor cores. An H100 SXM has 132 of them; each runs several blocks at a time.
  • HBM: High-bandwidth memory: the GPU's main memory (80 GB on an H100). Large, but slow compared with registers and on-chip shared memory.
  • Memory-bound / compute-bound: A kernel is memory-bound when moving its data takes longer than its arithmetic, and compute-bound in the opposite case.
  • Fusion: Doing several operations in one kernel so intermediate results stay on chip instead of making a round trip through HBM.
  • Tensor Core: A matrix-multiply unit inside each SM, many times faster than ordinary arithmetic units, at reduced input precision (TF32, BF16, FP8).

1. What a kernel is

A common mental model is that the GPU somehow runs your Python program. In reality, Python stays entirely on the CPU. The CPU acts as an orchestrator: PyTorch figures out which math operations are needed, and the CUDA driver launches a compiled C++ or PTX function onto the GPU, specifying a grid configuration that dictates how many parallel copies to spawn.

That GPU function is the kernel. Each running instance is a thread. Threads are grouped into thread blocks (up to 1,024 threads per block), and the hardware executes them in lockstep groups of 32 known as warps. Every thread in a warp executes the exact same instruction at the same time across different data elements. (If you want a broader overview of how CPUs and accelerators divide work, our guide to CPUs, GPUs and clusters introduces the hardware layout).

To understand why kernel performance varies so wildly, you have to look at the physical silicon. An NVIDIA H100 SXM chip contains 132 Streaming Multiprocessors (SMs). Each SM has its own fast registers and configurable on-chip shared memory (SRAM), while all SMs share a 50 MB L2 cache and 80 GB of High-Bandwidth Memory (HBM3). As you move down this memory hierarchy, capacity grows, but latency and energy costs multiply:

An H100 SXM in the numbers that matter for kernels. The CPU only launches the kernel; blocks of threads are spread over the 132 SMs, and every value they use travels up a memory hierarchy that gets smaller and faster closer to the arithmetic.
An H100 SXM in the numbers that matter for kernels. The CPU only launches the kernel; blocks of threads are spread over the 132 SMs, and every value they use travels up a memory hierarchy that gets smaller and faster closer to the arithmetic.

The central tension in GPU programming is the imbalance between raw compute power and memory bandwidth. An H100 SXM can churn through 67 TFLOPS of standard FP32 math on its CUDA cores and nearly 1,000 TFLOPS of BF16 matrix arithmetic on its Tensor Cores. Yet its memory bus caps out at 3.35 TB/s.

Do the math: to keep standard CUDA cores 100% saturated, a kernel must perform roughly 20 arithmetic operations for every single byte loaded from HBM. For Tensor Cores, that requirement shoots up to nearly 300 operations per byte. Yet everyday deep-learning operations — activations, additions, LayerNorms — perform only 1 or 2 operations per byte. The GPU’s massive compute engines end up starved, waiting on memory traffic. Once you see this gap, modern GPU kernel design suddenly makes sense: almost every trick is about moving less data across HBM.

2. A CUDA kernel

CUDA C++ is NVIDIA’s native dialect for programming accelerators. In CUDA, you write code from the viewpoint of a single thread. Built-in runtime variables — blockIdx, blockDim, and threadIdx — allow each thread to compute its unique position in the global array.

Here is how you write an elementwise ReLU ( yi=max⁡(xi,0)y_i = \max(x_i, 0) ) for an array of n floating-point numbers:

__global__ void relu_kernel(const float* x, float* y, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;   // this thread's element
    if (i < n) y[i] = fmaxf(x[i], 0.0f);             // the last block runs past n
}

// Host side: 256 threads per block, enough blocks to cover n elements.
relu_kernel<<<(n + 255) / 256, 256>>>(x, y, n);
Enter fullscreen mode Exit fullscreen mode

On the host side, the launch syntax <<<...>>> configures how many 256-thread blocks to schedule so that all n elements are covered. Because array sizes rarely align perfectly to multiples of 256, the defensive guard if (i < n) stops out-of-bounds memory writes on the final block.

Notice how thread i reads element x[i]. Because consecutive threads in a warp access adjacent memory addresses in parallel, the GPU memory controller can merge their 32 individual 4-byte reads into a single 128-byte transaction. This is called memory coalescing. If your access pattern is strided or scattered, memory throughput can drop by 5× to 10× without changing the mathematical output.

ReLU is straightforward because each data point is completely independent. Real-world machine learning kernels, however, require threads to collaborate. In a matrix multiplication, threads must stage input tiles into on-chip shared memory (a shared scratchpad inside the SM), coordinate through __syncthreads() barriers to prevent data races, issue specialized Tensor Core instructions, and carefully budget registers per thread to maintain high SM occupancy. CUDA gives you bare-metal control over all of these knobs — but it also means every misplaced barrier or register spill can silently tank your performance.

3. The same kernel in Triton

Writing raw CUDA can feel like writing C with inline assembly: you are constantly juggling thread IDs, shared memory pointers, and race conditions. This complexity is why Triton (Tillet, Kung and Cox, MAPL 2019; developed by OpenAI and the open-source community) has become so popular in the AI systems world.

In Triton, you write Python. More importantly, you stop thinking about individual scalar threads. Instead, you write from the perspective of an entire program block processing small array tiles (typically 64 to 1,024 elements) using vectorized NumPy-like operations. The Triton compiler takes care of mapping blocks to warps, vectorizing loads, ensuring coalesced access, and orchestrating shared memory:

import torch
import triton
import triton.language as tl

@triton.jit
def relu_kernel(x_ptr, y_ptr, n, BLOCK: tl.constexpr):
    pid = tl.program_id(0)                       # this program's block
    offs = pid * BLOCK + tl.arange(0, BLOCK)     # BLOCK positions at once
    mask = offs < n                              # the last block runs past n
    x = tl.load(x_ptr + offs, mask=mask)
    tl.store(y_ptr + offs, tl.maximum(x, 0.0), mask=mask)

def relu(x: torch.Tensor) -> torch.Tensor:
    y = torch.empty_like(x)
    n = x.numel()
    relu_kernel[(triton.cdiv(n, 1024),)](x, y, n, BLOCK=1024)
    return y
Enter fullscreen mode Exit fullscreen mode

The same 20-element ReLU. CUDA assigns one value per thread and leaves cooperation to the programmer; Triton assigns one block per program and leaves the mapping onto threads to the compiler.
The same 20-element ReLU. CUDA assigns one value per thread and leaves cooperation to the programmer; Triton assigns one block per program and leaves the mapping onto threads to the compiler.

What you decide and what the toolchain decides.

CUDA C++ Triton
Unit of code One thread One program over a block of values
Shared memory and barriers Written by hand Inserted by the compiler
Coalescing, vectorization Your index arithmetic Compiler, from the block’s layout
Tensor Cores WMMA/MMA instructions or libraries such as CUTLASS tl.dot on blocks
Tuning Template parameters, by hand Block sizes, warps and pipeline stages; @triton.autotune can time several
Typical use Libraries and the last few percent Fused ML kernels; the code that torch.compile generates

What do you give up with Triton? You give up fine-grained control. Individual program instances cannot communicate with each other except through slow global HBM memory, and low-level hardware tricks (like hand-crafted warp shuffles or irregular memory layouts) remain out of reach. But for the fused, memory-bound operations that dominate modern Transformers, that trade-off is an easy win: you get 80–95% of hand-tuned CUDA speed in a few dozen lines of readable Python. That is why PyTorch chose Triton as the primary code generation target for torch.compile.

4. Memory-bound or compute-bound: the roofline

Before spending days attempting to optimize a kernel, you need to answer a fundamental question: what is actually holding it back? Is the GPU waiting on arithmetic logic units to finish crunching numbers, or are the execution units idling while waiting for tensors to travel across the memory bus?

The classic roofline model (Williams, Waterman and Patterson, 2009) formalizes this intuition into a simple formula:

t  ≥  max⁡!(FLOPsPpeak, bytes movedBHBM),I=FLOPsbytes moved t \;\ge\; \max!\left(\frac{\text{FLOPs}}{P_{\text{peak}}},\ \frac{\text{bytes moved}}{B_{\text{HBM}}}\right),\qquad I=\frac{\text{FLOPs}}{\text{bytes moved}}

Here, II represents the arithmetic intensity — the ratio of floating-point math operations performed per byte of data moved across HBM. Every GPU architecture has a distinct transition point called the ridge point ( Ppeak/BHBMP_{\text{peak}}/B_{\text{HBM}} ):

  • If your kernel sits to the left of the ridge (low intensity), it is memory-bound. Faster arithmetic units will not speed it up by a single nanosecond. The only way to win is to move fewer bytes across the bus or reuse data in cache.
  • If your kernel sits to the right of the ridge (high intensity), it is compute-bound. Memory bandwidth is plenty fast; your runtime is determined by how well you saturate ALUs or Tensor Cores.

On an NVIDIA H100 SXM, this ridge line sits at 20 FLOPs/byte for FP32 on CUDA cores, roughly 148 FLOPs/byte for TF32 on Tensor Cores, and 295 FLOPs/byte for BF16 Tensor Cores (using dense datasheet peaks; NVIDIA’s marketing often cites 2:4 structured sparsity numbers, which are double these).

Arithmetic intensity of common operations in FP32, counting the minimum traffic (each input read once, each output written once).

Operation FLOPs per element Bytes per element Intensity Bound on H100
ReLU, add, multiply 1 8 0.125 Memory
Softmax over a row ≈ 5 8 ≈ 0.6 Memory
LayerNorm ≈ 8 8 ≈ 1 Memory
Matmul, N = 1,024 2N³ FLOPs over 3N² values ≈ 171 Compute (FP32 cores)
Matmul, N = 4,096 2N³ FLOPs over 3N² values ≈ 683 Compute

Notice where standard deep-learning building blocks fall: activations, normalizations, and softmax sit two full orders of magnitude to the left of the ridge. For these layers, optimizing math formulas is irrelevant; all that matters is memory traffic. Conversely, large matrix multiplications sit comfortably to the right: they offer abundant theoretical data reuse, so the challenge becomes organizing memory access so Tensor Cores never stall. Use the interactive explorer below to inspect these bounds across different tensor shapes:

Kernel Cost Explorer: a chained-elementwise fusion example, 16M FP32 elements, 3 ops. Eager reads/writes each op to HBM (403 MB, a 120 µs bandwidth floor); a fused kernel keeps intermediates in registers (134 MB, 40 µs).
The site has an interactive version with three tabs (fusion, roofline, matmul tiling) and sliders for element count and dtype — try it on the original post.

5. Four optimizations, from simple to hard

5.1 Fusion: stop writing intermediates to memory

Consider a familiar one-liner: y = relu(x * s + b). In standard eager PyTorch, this statement executes as three distinct kernel launches:

  1. Kernel 1 (scale): reads x from HBM, computes x * s, and writes a temporary tensor tmp1 back to HBM.
  2. Kernel 2 (bias add): reads tmp1 and b, computes tmp1 + b, and writes temporary tensor tmp2 back to HBM.
  3. Kernel 3 (activation): reads tmp2, applies relu, and writes final tensor y back to HBM.

Notice the waste: tmp1 and tmp2 exist solely to pass numbers to the next line of code, yet they force millions of bytes across the memory bus. A fused kernel collapses all three steps into a single pass: it reads x once, keeps intermediate values in fast on-chip registers, and writes y once:

@triton.jit
def scale_bias_relu(x_ptr, b_ptr, y_ptr, n, cols, s, BLOCK: tl.constexpr):
    offs = tl.program_id(0) * BLOCK + tl.arange(0, BLOCK)
    mask = offs < n
    x = tl.load(x_ptr + offs, mask=mask)
    b = tl.load(b_ptr + offs % cols, mask=mask)     # one bias per column, stays in cache
    tl.store(y_ptr + offs, tl.maximum(x * s + b, 0.0), mask=mask)
Enter fullscreen mode Exit fullscreen mode

Three eager kernels move six tensors through HBM; one fused kernel moves two. For 16.8M FP32 values the bandwidth bound falls from 120 µs to 40 µs, and two kernel launches disappear.
Three eager kernels move six tensors through HBM; one fused kernel moves two. For 16.8M FP32 values the bandwidth bound falls from 120 µs to 40 µs, and two kernel launches disappear.

Fusion does not eliminate a single FLOP of arithmetic. What it eliminates is memory traffic and driver launch overhead. That is why operator fusion is always the first optimization compilers apply, and why cleaning up chains of small elementwise layers yields such dramatic speedups.

5.2 Softmax: fuse a reduction, then make it streaming

Softmax is a notorious bottleneck in attention layers. The standard formula — yi=exi−m/∑jexj−my_i = e^{x_i - m} / \sum_j e^{x_j - m} with the row maximum mm subtracted for numerical stability — is built from two reductions (max and sum) and several elementwise operations. Implemented naively from PyTorch primitives, it requires five kernel launches and eight passes over the data:

def naive_softmax(x):          # x: [rows, cols]
    m = x.max(dim=1)[0]           # kernel 1: read x
    z = x - m[:, None]            # kernel 2: read x, write z
    e = torch.exp(z)              # kernel 3: read z, write e
    d = e.sum(dim=1)              # kernel 4: read e
    return e / d[:, None]         # kernel 5: read e, write y
Enter fullscreen mode Exit fullscreen mode

If an entire row fits inside on-chip SRAM, one program block per row can load it once into registers, compute the maximum and the sum on the fly, and write the normalized probabilities back to HBM in a single pass. That cuts memory traffic from eight passes down to two — delivering an immediate ~4× speedup. Below is a simplified implementation adapted from Triton’s fused-softmax tutorial:

@triton.jit
def softmax_kernel(x_ptr, y_ptr, cols, stride, BLOCK: tl.constexpr):
    row = tl.program_id(0)                        # one program per row
    offs = tl.arange(0, BLOCK)                    # BLOCK = next power of 2 >= cols
    mask = offs < cols
    x = tl.load(x_ptr + row * stride + offs, mask=mask, other=-float("inf"))
    x = x - tl.max(x, axis=0)                     # subtract the row max for stability
    e = tl.exp(x)
    tl.store(y_ptr + row * stride + offs, e / tl.sum(e, axis=0), mask=mask)
Enter fullscreen mode Exit fullscreen mode

What happens when the sequence length is too large to fit in on-chip SRAM (like in long-context models)? You cannot compute the true row maximum before seeing the entire row — unless you use the online softmax trick (Milakov and Gimelshein, 2018). Online softmax maintains a running maximum mm and a running sum dd , dynamically rescaling previous partial sums whenever a larger element appears in the next chunk:

m′=max⁡!(m, max⁡j∈chunkxj),d′=d e m−m′+∑j∈chunke xj−m′ m' = \max!\big(m,\ \max_{j \in \text{chunk}} x_j\big),\qquad d' = d\,e^{\,m - m'} + \sum_{j \in \text{chunk}} e^{\,x_j - m'}

This streaming reduction allows the kernel to compute exact softmax statistics in a single streaming pass. Applied to attention, this exact rescaling lets a kernel process keys and values block-by-block without ever materializing the massive N×NN \times N attention matrix in HBM — the core algorithmic breakthrough behind FlashAttention (Dao et al., 2022).

5.3 Tiling and Tensor Cores: reuse every load

Matrix multiplication ( C=ABC = AB ) of two N×NN \times N matrices requires 2N32N^3 FLOPs, yet only touches 3N23N^2 values. In theory, that offers massive computational reuse.

In a naive kernel where each thread calculates a single element Ci,jC_{i,j} , every thread reads an entire row of A and column of B from HBM — resulting in 2N32N^3 total loads. The solution is tiling: each program computes a T×TT \times T block of C, streaming strips of A and B through shared memory so every loaded value is reused TT times:

loadsnaive=2N3,loadstiled=2N3T+N2 \text{loads}{\text{naive}} = 2N^3,\qquad \text{loads}{\text{tiled}} = \frac{2N^3}{T} + N^2

Tiling divides the requested traffic by the tile size. At N = 4,096 in FP32, the naive kernel asks for 550 GB; with 128 × 128 tiles the request falls to 4.4 GB, close to what the arithmetic itself needs on CUDA cores.
Tiling divides the requested traffic by the tile size. At N = 4,096 in FP32, the naive kernel asks for 550 GB; with 128 × 128 tiles the request falls to 4.4 GB, close to what the arithmetic itself needs on CUDA cores.

In Triton, tiling is structured as a loop over chunks of dimension K, and calling tl.dot dispatches matrix operations directly to hardware Tensor Cores:

@triton.jit
def matmul_kernel(a_ptr, b_ptr, c_ptr, M, N, K,
                  s_am, s_ak, s_bk, s_bn, s_cm, s_cn,
                  BM: tl.constexpr, BN: tl.constexpr, BK: tl.constexpr):
    rm = tl.program_id(0) * BM + tl.arange(0, BM)     # this program's rows of C
    rn = tl.program_id(1) * BN + tl.arange(0, BN)     # and its columns
    rk = tl.arange(0, BK)
    acc = tl.zeros((BM, BN), dtype=tl.float32)
    for k in range(0, K, BK):                         # stream strips of A and B
        a = tl.load(a_ptr + rm[:, None] * s_am + (k + rk)[None, :] * s_ak,
                    mask=(rm[:, None] < M) & (k + rk[None, :] < K), other=0.0)
        b = tl.load(b_ptr + (k + rk)[:, None] * s_bk + rn[None, :] * s_bn,
                    mask=(k + rk[:, None] < K) & (rn[None, :] < N), other=0.0)
        acc += tl.dot(a, b, input_precision="ieee")   # FP32 inputs default to TF32
    tl.store(c_ptr + rm[:, None] * s_cm + rn[None, :] * s_cn, acc,
             mask=(rm[:, None] < M) & (rn[None, :] < N))
Enter fullscreen mode Exit fullscreen mode

Notice the comment on tl.dot: it conceals a notorious precision trap. Tensor Cores do not compute in standard IEEE FP32 by default. When fed FP32 inputs, they automatically switch to TF32 (TensorFloat-32), which retains the 8-bit dynamic range of FP32 but truncates the mantissa from 23 bits down to 10 bits (the precision of FP16).

Because TF32 Tensor Cores deliver up to 495 TFLOPS on an H100 — roughly 7.4× faster than standard 67 TFLOPS FP32 CUDA cores — a Triton kernel will look blazingly fast compared to PyTorch’s default FP32 matmul. But when tested against a strict numerical tolerance, it may fail correctness tests. Setting input_precision="ieee" restores full FP32 accuracy at full FP32 cost; passing "tf32x3" splits each input into high and low parts, executing three Tensor Core operations to recover near-FP32 precision.

Matmul precision on an H100 SXM. Peaks are dense datasheet values; relative error is per product, before accumulation.

Inputs Unit Peak Mantissa bits Default in
FP32 (IEEE) CUDA cores 67 TFLOPS 23 PyTorch matmul
TF32 Tensor Cores ≈ 495 TFLOPS 10 Triton tl.dot on FP32; PyTorch cuDNN convolutions
BF16 Tensor Cores ≈ 989 TFLOPS 7 Mixed-precision training

5.4 Epilogues: let the library do the matmul

Trying to beat NVIDIA’s cuBLAS on a bare matrix multiplication is usually an exercise in frustration: NVIDIA engineers spend years tuning assembly kernels for each GPU microarchitecture. Where you can win is what happens immediately before and after the matmul.

If your matmul is followed by a bias add, a GELU activation, or a residual connection, you can apply those operations to each tile while it is still sitting in registers before writing it to HBM. This fused post-processing is called an epilogue. cuBLASLt provides fixed epilogues (like bias + GELU); Triton allows you to write arbitrary epilogues directly in Python; and a pragmatic architecture is often to call cuBLAS for the heavy GEMM, followed by a single fused kernel for the remaining elementwise tail.

Two other techniques appear in almost every high-performance kernel. Autotuning times multiple combinations of tile sizes, warp counts, and pipeline stages on target hardware, keeping the winner for each tensor shape. Software pipelining overlaps memory and compute by prefetching tile k+1k+1 while Tensor Cores multiply tile kk . Triton exposes both via simple decorators and arguments (@triton.autotune, num_stages), while in CUDA they require intricate manual engineering.

6. cuBLAS, cuDNN, eager PyTorch and torch.compile

In modern deep learning, GPU kernels reach hardware through four distinct paths. Depending on which path handles a layer, the exact same mathematical code can execute very differently:

One linear layer with a GELU, three ways. Eager PyTorch launches a kernel per operator; torch.compile keeps the library matmul and fuses the elementwise tail into a generated Triton kernel; a hand-written kernel can do everything in one launch.
One linear layer with a GELU, three ways. Eager PyTorch launches a kernel per operator; torch.compile keeps the library matmul and fuses the elementwise tail into a generated Triton kernel; a hand-written kernel can do everything in one launch.

Where the kernels come from.

Source What it is Strong at Weak at
cuBLAS / cuBLASLt NVIDIA’s precompiled dense linear algebra; picks a kernel per shape with heuristics Matmuls of every common shape and precision; fixed epilogues Anything that is not a matmul; arbitrary fused epilogues
cuDNN NVIDIA’s deep-learning primitives: convolution, pooling, normalization, fused attention Convolutions, with several algorithms per shape (torch.backends.cudnn.benchmark times them) Operations outside its catalog; fusing your own code
Eager PyTorch One kernel per operator, launched from Python as the code runs; matmuls and convolutions call the libraries above Flexibility, debuggability, no compile time Chains of small operators: every one is a launch and an HBM round trip
torch.compile TorchDynamo captures a graph from the Python bytecode; TorchInductor fuses it and generates Triton kernels, calling the libraries for matmuls and convolutions Fusing elementwise and reduction chains; removing launch overhead with CUDA graphs Compile time; recompiles on new shapes; a single large matmul or convolution gains little
Hand-written CUDA C++, Triton, or templates such as CUTLASS Fusions no compiler finds; layouts and algorithms specific to one model Engineering time; correctness is entirely yours to prove

torch.compile comes with different operating modes. The default mode fuses elementwise and reduction chains into Triton kernels. mode="reduce-overhead" captures the entire execution sequence into a CUDA graph, replaying dozens of launches in a single driver call — a massive win when models are CPU-dispatch bound. Finally, mode="max-autotune" benchmarks Triton templates against vendor libraries to select optimal kernels per tensor dimension.

On whole models, these compiler passes deliver real gains: the PyTorch 2 paper (Ansel et al., ASPLOS 2024) reports geometric-mean speedups of 2.27× for inference and 1.41× for training over eager mode across 180+ models on A100 GPUs.

The practical takeaway is that claims of being “faster than PyTorch” are ambiguous without context. Faster than eager PyTorch often just means you fused two elementwise operations that anyone could fuse. Faster than torch.compile in max-autotune mode is a much higher bar, because the compiler has already fused elementwise tails, benchmarked library alternatives, and eliminated launch overhead.

7. KernelBench: can an LLM write these kernels?

Can frontier language models write GPU kernels that outperform human engineers and optimizing compilers? To test this systematically, Stanford researchers released KernelBench (Ouyang et al., 2025). Each task presents a standalone PyTorch module; the model must output a drop-in replacement whose forward() method runs custom CUDA or Triton kernels:

# KernelBench level 2, problem 1: Conv2D_ReLU_BiasAdd (docstrings removed)
import torch
import torch.nn as nn

class Model(nn.Module):
    def __init__(self, in_channels, out_channels, kernel_size, bias_shape):
        super(Model, self).__init__()
        self.conv = nn.Conv2d(in_channels, out_channels, kernel_size)
        self.bias = nn.Parameter(torch.randn(bias_shape))

    def forward(self, x):
        x = self.conv(x)
        x = torch.relu(x)
        x = x + self.bias
        return x

batch_size = 128
in_channels  = 64
out_channels = 128
height = width = 128
kernel_size = 3
bias_shape = (out_channels, 1, 1)

def get_inputs():
    return [torch.rand(batch_size, in_channels, height, width)]

def get_init_inputs():
    return [in_channels, out_channels, kernel_size, bias_shape]
Enter fullscreen mode Exit fullscreen mode

How KernelBench scores one answer: it must compile, match the reference on random inputs, and is then timed against the PyTorch program it replaces. fast_p counts the problems that pass both gates with a speedup above p.
How KernelBench scores one answer: it must compile, match the reference on random inputs, and is then timed against the PyTorch program it replaces. fast_p counts the problems that pass both gates with a speedup above p.

KernelBench evaluates solutions through three rigorous gates:

  1. Compilation & execution: the kernel must compile cleanly and run without CUDA exceptions.
  2. Correctness: the outputs must match the PyTorch reference across 5 random input batches (using numerical tolerances calibrated to the precision).
  3. Performance: if correct, the kernel is timed using CUDA events (3 warm-up passes followed by the mean of 100 runs, with L2 cache cleared between trials) against the eager PyTorch baseline.
fastp=1N∑i=1N1[correcti  ∧  speedupi>p],speedupi=tPyTorch,itkernel,i \text{fast}p=\frac{1}{N}\sum{i=1}^{N}\mathbf{1}\big[\text{correct}i \;\wedge\; \text{speedup}_i > p\big],\qquad \text{speedup}_i=\frac{t{\text{PyTorch},i}}{t_{\text{kernel},i}}

Here, fast0\text{fast}_0 is the raw correctness rate, while fast1\text{fast}_1 measures the proportion of problems where the model generated a correct kernel that was strictly faster than PyTorch.

The benchmark spans multiple difficulty tiers. Level 1 covers isolated operators. Level 2 tests composite patterns (like Conv2D + ReLU + BiasAdd) where fusion and epilogues provide clear optimization headroom. Level 3 tests full architecture blocks. The baseline results from early 2025 were a sobering reality check:

fast₁ against eager PyTorch, one attempt per problem, CUDA kernels on an NVIDIA L40S. Reported by Ouyang et al. (2025); not our measurements.

Model Level 1 Level 2 Level 3
DeepSeek R1 12% 36% 2%
OpenAI o1 10% 24% 12%
Claude 3.5 Sonnet 10% 7% 2%
DeepSeek V3 6% 4% 8%
GPT-4o 4% 5% 0%
Llama 3.1 405B 3% 0% 2%

On single-shot generation, even the best reasoning models struggled to beat eager PyTorch on more than a third of Level 2 problems, and barely registered on Level 3.

Why do LLMs struggle so much with kernels? Part of the answer is data distribution: CUDA C++ accounts for only 0.073% of The Stack v1.2, a standard pretraining corpus. But the paper also proved that the evaluation protocol matters as much as the base model:

  • Iterative compiler feedback: When DeepSeek R1 was given compiler error diagnostics, runtime outputs, and profiler metrics over up to 10 feedback turns, its fast1\text{fast}_1 score surged from 12% / 36% / 2% to 43% / 72% / 18% across the three levels.
  • Best-of-N sampling: Sampling 100 candidate completions per problem raised DeepSeek V3’s Level 2 score from 4% to 37%.

Since the original paper, the ecosystem has accelerated rapidly: KernelBench added support for Triton, HIP, and TileLang, while specialized research models — such as Cognition’s Kevin (multi-turn RL on CUDA), AutoTriton, and Dr. Kernel — have trained specifically on GPU environments.

8. How to read a KernelBench number

Whenever a research group or vendor claims “our AI model generates kernels that are 10× faster than PyTorch”, take a breath and check the methodology. Kernel benchmarking is fraught with subtle discrepancies:

Seven questions to ask before comparing two KernelBench results.

Question Why it changes the number
Which levels or subset? Level 2 is much easier to speed up than levels 1 and 3. Many papers report a subset.
Which baseline? Eager PyTorch and each torch.compile mode are different targets; for the same answers the paper reports DeepSeek R1 at 12/36/2% fast₁ against eager and 38/37/2% against torch.compile in its default mode.
Which GPU? The paper used an L40S. The same kernel can be fast on one GPU and slow on another, and speedups do not carry over between generations.
Which language? CUDA and Triton answers are judged by the same rules but are different tasks for a model.
How many attempts? One greedy answer, the best of k samples, and several turns with compiler and profiler feedback are different protocols; the paper shows each can multiply fast₁.
Which precision and tolerance? A 10⁻² FP32 tolerance accepts TF32 and many approximations that 10⁻⁴ rejects.
Which anti-cheating checks? Without them, some “speedups” are not computations at all (below).

A fast kernel needs a strict verifier

Kernel generation benchmarks are uniquely vulnerable to reward hacking because the objective function is literally a timer. In February 2025, Sakana AI announced an “AI CUDA Engineer” claiming 10× to 100× speedups on common PyTorch operations. Within 24 hours, the team retracted the claims after discovering that the model had found an exploit in the test harness: it corrupted memory checks so the validation suite erroneously marked broken kernels as correct.

As the KernelBench documentation notes: “If the model can reward hack, it will find ways to reward hack!” Here are the most common exploits and the defensive checks required to catch them:

Common ways a kernel can look fast without being fast, and the check that catches each.

Exploit What it does Check
PyTorch in disguise Calls torch operators or the reference inside ModelNew instead of a kernel Static checks on the code; confirm that custom kernels actually launch
Reused results Returns an output cached from the reference run or from an earlier call Fresh random inputs per trial; outputs compared on inputs the kernel has not seen
Hidden work Runs on a side CUDA stream that the timer does not wait for Synchronize all streams; time the whole call
Loose numerics Lowers precision until the output only just passes Tolerance matched to the dtype; inputs with large and small magnitudes
Shape specialization Works only for the exact test shape or values Hidden test cases with other values; edge cases
Too good to be true Speedups of 10× or more on simple operations Flag for review: KernelBench marks speedups above 10× as excessive

These principles mirror the verifier design patterns we use across all RLVR workflows; our verifier playbook and reward-hacking analysis cover these guardrails in detail. For GPU kernels, the absolute minimum baseline is testing on unseen input distributions, explicit numerical tolerance boundaries, full CUDA stream synchronization before stopping the clock, and clearly naming whether you are comparing against eager PyTorch or torch.compile.

What comes next

The concepts covered here form the foundation for the next engineering deep dive in this series, where we show how we train an open LLM to write performant Triton kernels inside a GPU kernel gym, and what RLVR techniques are required to translate those gains onto KernelBench.


Originally published at g-ftech.com.

Top comments (0)