Profiling, GPU utilisation and writing a fused kernel

A slow model is slow for one of three reasons: the GPU is out of arithmetic, out of memory bandwidth, or sitting idle waiting for work. This page shows how to tell which one with a profiler, how to debug a GPU that looks two-thirds idle, and how to write and fairly benchmark a small fused kernel.

Source: Vishal Singh's post of 25 inference-engineering problems · saved in sources/09-vishal-inference-25.md

Going deeper: everything below is the full version, with some maths. On a first read, skim the headings, the interview answers and the "Why it matters for robotics" box, and come back for the rest later.

The numbers this page uses

Every worked answer below uses one GPU, so the arithmetic is comparable. These figures come from NVIDIA's H100 product page, SXM column. The PCIe-class variant listed beside it differs (different memory size and bandwidth), so don't mix the two.

Kernel launch costs, CPU step times and request rates below are assumptions chosen to make the arithmetic clear. They are not measurements. Measure your own.

Problem 13 · Profile a slow model with Nsight: compute, memory bandwidth, or launch overhead?

What the interviewer is testing: whether you have a method, not a hunch. Can you pick the right tool, read a timeline, use the roofline model, and name the fix for each kind of bottleneck?

The concepts

Two tools, two zoom levels. Nsight Systems (nsys) records a timeline of the whole program: CPU threads, CUDA API calls, every kernel, memory copies, and communication between GPUs. It answers "where does wall-clock time go, and is the GPU ever idle?" Nsight Compute (ncu) replays individual kernels many times with hardware counters. It answers "why is this kernel slow?" Always start with Systems. Running Compute on everything is slow and drowns you in detail.

# 1. Timeline of a few hundred decode steps
nsys profile -t cuda,nvtx,osrt --cuda-graph-trace=node -o trace python bench.py

# 2. Deep dive on the kernels that dominate the timeline
ncu --set full -k regex:rms -c 20 -o rms_report python bench.py

Wrap the phases of your code in NVTX ranges (torch.cuda.nvtx.range_push("prefill") … range_pop()) so the timeline shows "prefill", "decode step", "sampling" as labelled bars. Without them you are reading thousands of anonymous kernel names.

The roofline model. A kernel does some floating-point operations (FLOPs) and moves some bytes to and from HBM. Its arithmetic intensity is the ratio. Its best possible speed is capped by whichever resource runs out first:

arithmetic intensity  I = FLOPs / bytes moved            (FLOP per byte)
attainable FLOP/s       = min( peak FLOP/s ,  I × bandwidth )
ridge point          I* = peak FLOP/s / bandwidth

H100 SXM, dense BF16:  I* = 989e12 / 3.35e12 ≈ 295 FLOP/byte
I < 295  → memory-bound   (more bandwidth or fewer bytes helps)
I > 295  → compute-bound  (more FLOP/s or fewer FLOPs helps)

Launch overhead is the third case, and the roofline doesn't show it. Each kernel launch costs CPU time in Python, the framework and the driver. If the GPU finishes a kernel before the CPU has queued the next one, the GPU sits idle. Nsight Systems shows this as gaps between kernels on the GPU row, while the CPU row is busy.

A worked answer with real numbers

Take one linear layer of a model with hidden size d = 4,096, weights in BF16 (2 bytes). Feed it N tokens at once: N is the batch size during decode, or the chunk length during prefill.

FLOPs = 2 × N × d²                       (multiply + add per weight per token)
bytes = 2 × (d² + 2·N·d)                 (read weights, read input, write output)
I     = 2·N·d² / (2·(d² + 2·N·d)) = N·d / (d + 2N)

Decode, N = 1:
  I = 4096 / 4098 ≈ 1.0 FLOP/byte          → far below 295, memory-bound
  time ≈ bytes / bandwidth = 33.6 MB / 3.35 TB/s ≈ 10 µs
  attainable = 1.0 × 3.35 TB/s ≈ 3.35 TFLOP/s → about 0.3% of the 989 peak

Prefill chunk, N = 2,048:
  I = 2048 × 4096 / (4096 + 4096) = 1,024 FLOP/byte → compute-bound
  compute time = 2 × 2048 × 4096² / 989e12 ≈ 69 µs
  memory time  = 67 MB / 3.35 TB/s ≈ 20 µs  → compute wins

Crossover: N·d/(d + 2N) = 295  →  N ≈ 345 tokens per step

So for this layer on this GPU, any step with fewer than about 345 tokens is memory-bound. Single-user decode is deep in the memory-bound corner. Attention during decode is worse still, because every sequence reads its own KV cache (see page 02).

The whole-model version gives a hard upper bound on single-stream decode speed:

8B parameters × 2 bytes = 16 GB read per decode step
16 GB / 3.35 TB/s ≈ 4.8 ms per token  →  at most ~209 tokens/s at batch 1
(real numbers are lower: KV-cache reads, imperfect bandwidth use, overheads)

Launch overhead, with assumed numbers. Suppose a decode step launches about 15 kernels per layer across 32 layers: 480 launches. Assume each costs about 5 µs of CPU time in eager PyTorch. That's 2.4 ms of CPU work per step. For the 8B model the GPU needs about 4.8 ms, so the CPU just keeps up. Shrink the model to 1B parameters (2 GB, about 0.6 ms of GPU time per step) and the same 2.4 ms of launches now dominates: the GPU is idle about three-quarters of the time. Small models are where launch overhead bites.

Reading the profile: a decision table

Trade-offs

• Nsight Systems vs Nsight Compute: Systems is low overhead and shows the whole picture; Compute is precise per kernel but replays each kernel, so it's slow and changes timings. Use Systems to choose which kernels to send to Compute.
• CUDA graphs: remove launch cost by replaying a recorded sequence of kernels, but need fixed shapes, so engines capture one graph per batch-size bucket. That costs memory and startup time.
• Profiling a server vs a script: a script is clean and repeatable; a live server shows real scheduler gaps. Profile both.

Common mistakes

• Running Nsight Compute on the whole program first, and drowning in data.
• Using the "with sparsity" FLOP figure from a spec sheet as the peak for a dense model, which doubles the apparent headroom.
• Profiling a cold start, so the timeline is dominated by compilation, autotuning and memory allocation.
• Calling a kernel "slow" without comparing it to its roofline bound. A kernel at 90% of peak bandwidth is done; only the algorithm can help.

Problem 14 · Your GPU sits at 30% utilisation under load. Debug it

What the interviewer is testing: whether you know what the utilisation number actually measures, and whether you debug from the outside in (load, scheduler, CPU) before blaming kernels.

The concepts

What "GPU utilisation" means. The figure in nvidia-smi is the share of a sample period during which at least one kernel was running. It says nothing about how much of the GPU that kernel used. A tiny kernel on one SM counts as 100% busy. So 30% means the GPU ran nothing for about 70% of the time. That points at starvation, not slow kernels. For how well the busy time was used, DCGM's profiling metrics are better: SM activity, SM occupancy, tensor pipe activity and DRAM activity.

Little's law links load to batch size. The average number of requests in the system equals arrival rate × time each spends there:

in-flight requests = arrival rate × mean request duration
GPU busy fraction  = GPU time per step / (GPU time + non-overlapped CPU time per step)

A worked answer with real numbers

Two common stories, with assumed numbers.

Story A: the load is too light. The engine is configured for up to 64 concurrent sequences. Each request lasts about 10 s (prompt plus a few hundred output tokens). The load tester sends 1 request per second.

in-flight = 1 req/s × 10 s = 10 sequences   (the engine allows 64)
to fill 64 slots: 64 / 10 s = 6.4 req/s needed

The GPU is not broken. The test is too gentle, or a gateway in front is limiting concurrency. Check the engine's running and waiting counts. If nothing is waiting, the GPU isn't the bottleneck.

Story B: the CPU is in the critical path. Each decode step needs 12 ms of GPU work. Between steps, the CPU spends 28 ms on scheduling, sampling, detokenising and streaming, and none of it overlaps with the GPU.

busy fraction = 12 / (12 + 28) = 30%
step time     = 40 ms  →  25 steps/s

If CPU work overlaps the next GPU step:
step time = max(12, 28) = 28 ms → 43% busy, 36 steps/s
If the CPU work is also cut to 8 ms:
step time = max(12, 8) = 12 ms → ~100% busy, 83 steps/s

Overlap alone isn't enough here, because the CPU part is longer than the GPU part. You also have to shrink it: move detokenisation off the main loop, sample on the GPU, batch the streaming writes.

The debugging ladder

1. Is there enough work? Look at the waiting queue and running batch size. Empty queue plus small batch means a load or routing problem.
2. Is the scheduler capping the batch? Limits on concurrent sequences or tokens per step, or a KV cache that's full and causing preemptions, keep batches small even when requests are waiting.
3. Is the CPU in the path? Nsight Systems: GPU gaps lined up with Python, tokenisation, sampling or network work on the CPU threads. Also check the CPU itself is not throttled or shared.
4. Is it waiting on transfers or other GPUs? Host-to-device copies, synchronising calls such as .item() in the loop, or tensor-parallel all-reduces waiting on a slow peer.
5. Only then, the kernels. Many small kernels with gaps: launch overhead (Problem 13).

Trade-offs

Common mistakes

• Treating nvidia-smi utilisation as "how much of the GPU is used". 100% can still mean one busy SM.
• Jumping straight to kernel optimisation when the queue is empty.
• Testing with a closed-loop client at low concurrency and concluding the server can't scale.
• Adding a hidden .item() or .cpu() in the decode loop, which forces the CPU to wait for the GPU every step.

Problem 15 · Write a fused RMSNorm kernel and benchmark it against PyTorch

What the interviewer is testing: can you write a correct GPU kernel (indexing, masking, reductions, precision), reason about why fusion helps, and benchmark honestly?

The concepts

RMSNorm (Zhang and Sennrich, 2019) is LayerNorm without the mean subtraction and bias. Llama-style models apply it twice per layer. For a row x of length d:

rms(x) = sqrt( (1/d) × Σ xᵢ²  + ε )
yᵢ     = xᵢ / rms(x) × wᵢ

It does about 4 FLOPs per element (square, add, two multiplies). In BF16 a fused version reads 2 bytes and writes 2 bytes per element, so its intensity is about 1 FLOP per byte. That's far below the H100 ridge of about 295. RMSNorm is memory-bound, so its speed is set by bytes moved. That's why fusion helps.

Written naively in eager PyTorch, each operation is its own kernel and writes its result back to HBM:

x.pow(2)        read x, write tmp1          4 bytes/elem in BF16
.mean(-1)       read tmp1                   2
x * rstd        read x, write tmp2          4
tmp2 * w        read tmp2, write y          4
                                            ── ~14 bytes/elem
fused kernel    read x once, write y once   4 bytes/elem

That simple count says unfused moves about 3.5× more bytes. Real eager code often upcasts to FP32 for the sum, which adds more traffic; PyTorch's built-in torch.nn.functional.rms_norm and torch.compile may already fuse it. So the fair baselines are eager, built-in and compiled, not eager alone.

A Triton kernel sketch

One program instance handles one row. The whole row fits in one block for typical hidden sizes (up to a few tens of thousands), so the reduction is a single tl.sum. Accumulate in FP32 even when inputs are BF16.

import torch, triton, triton.language as tl

@triton.jit
def rmsnorm_kernel(X, W, Y, stride, N, eps, BLOCK: tl.constexpr):
    row = tl.program_id(0)
    cols = tl.arange(0, BLOCK)
    mask = cols < N
    x = tl.load(X + row * stride + cols, mask=mask, other=0.0).to(tl.float32)
    w = tl.load(W + cols, mask=mask, other=0.0).to(tl.float32)
    rstd = tl.rsqrt(tl.sum(x * x, axis=0) / N + eps)   # masked lanes add 0
    y = x * rstd * w
    tl.store(Y + row * stride + cols, y.to(Y.dtype.element_ty), mask=mask)

def rmsnorm(x, w, eps=1e-6):
    N = x.shape[-1]
    x2 = x.reshape(-1, N).contiguous()
    y = torch.empty_like(x2)
    BLOCK = triton.next_power_of_2(N)
    warps = 8 if BLOCK >= 4096 else 4
    rmsnorm_kernel[(x2.shape[0],)](x2, w, y, x2.stride(0), N, eps,
                                   BLOCK=BLOCK, num_warps=warps)
    return y.view_as(x)

Points worth saying aloud: the mask handles hidden sizes that aren't a power of two; other=0.0 keeps masked lanes out of the sum; dividing by N (not BLOCK) gives the true mean; the FP32 accumulator avoids BF16 rounding in the sum. For very wide rows you would loop over chunks of the row, as the Triton layer-norm tutorial does. In a real model, the bigger win is fusing the residual add into the same kernel (read x and residual, write the sum and the normalised output), which serving engines ship as a standard op.

The benchmarking method

from triton.testing import do_bench
import torch.nn.functional as F

def eager(x, w, eps=1e-6):
    xf = x.float()
    return (xf * torch.rsqrt(xf.pow(2).mean(-1, keepdim=True) + eps) * w.float()).to(x.dtype)

builtin  = lambda x, w: F.rms_norm(x, (x.shape[-1],), w, 1e-6)
compiled = torch.compile(eager)

d = 4096
w = torch.randn(d, device="cuda", dtype=torch.bfloat16)
for M in [1, 16, 256, 4096, 32768]:                    # decode-like → prefill-like
    x = torch.randn(M, d, device="cuda", dtype=torch.bfloat16)
    ref = eager(x, w)
    for name, f in [("eager", eager), ("builtin", builtin),
                    ("compiled", compiled), ("triton", rmsnorm)]:
        torch.testing.assert_close(f(x, w), ref, atol=2e-2, rtol=2e-2)
        ms = do_bench(lambda: f(x, w))                 # warm-up + repeats, median
        gbs = 2 * x.numel() * x.element_size() / ms / 1e6   # read x + write y
        print(f"{M:6d} {name:9s} {ms*1e3:8.1f} us  {gbs:7.0f} GB/s")

The rules that make this fair:

• Check correctness first, against an FP32 reference, with a tolerance that fits BF16.
• Report effective bandwidth, not just time. Count only the bytes a perfect kernel must move (read x, write y). Then compare with 3.35 TB/s on an H100 SXM. A kernel near 80–90% of peak has little left to gain.
• Sweep shapes. At M = 1 the kernel takes microseconds and launch overhead dominates every version. At large M bandwidth dominates. Report both regimes.
• Beat the cache. Small inputs can sit in L2 cache between repeats, which inflates bandwidth above HBM speed. Use sizes well beyond L2 for the bandwidth claim, or flush between runs.
• Warm up. Exclude Triton and torch.compile compilation from timings. Use CUDA events or do_bench, never Python time.time() without synchronising.
• Confirm with Nsight Compute. The fused kernel should show memory throughput near peak and one kernel instead of several.

What to expect, as a roofline estimate rather than a measurement: at M = 4,096 and d = 4,096 in BF16, a perfect fused kernel moves 67 MB, about 20 µs at 3.35 TB/s. The 14-bytes-per-element eager version would need about 70 µs at full bandwidth, plus four launches. Whether your kernel beats the built-in or compiled version is an empirical question. Run it.

Trade-offs

Common mistakes

• Accumulating the sum of squares in BF16, which loses precision on wide rows.
• Dividing by the block size instead of the true row length when the row is padded.
• Benchmarking only against eager PyTorch and calling it a win, when torch.compile may match it.
• Reporting a speedup from one shape, or from inputs small enough to live in L2 cache.

Practise it

On a free Colab GPU (a T4 or better), in an afternoon:

1. Write the kernel. Paste the Triton sketch above, run the benchmark loop, and record time and effective GB/s per shape for all four versions. Look up your GPU's memory bandwidth on NVIDIA's product page and compute the fraction of peak.
2. Add the residual. Extend the kernel to take a residual tensor, write x + residual back, and normalise the sum. Count bytes before and after, and check the measured change against the count.
3. Profile. Run nsys profile on a loop of 100 calls to each version. Count kernels per call and look at the gaps at M = 1. Then run ncu on the Triton kernel and read the memory throughput figure.
4. Create the 30% problem on purpose. Serve a small model (about 0.5–1B parameters) with vLLM, drive it with a load tool at concurrency 1, 8 and 64, and log the engine's running-request count, nvidia-smi utilisation and tokens per second. Then add a deliberate .item() in a custom sampling hook (or simulate CPU work) and see utilisation fall.

What to measure: effective bandwidth vs peak, kernels per call, GPU idle gaps in microseconds, and how throughput and utilisation change with concurrency.

In this course

• Module 08 · Tensor maths and autograd: Triton programming basics and FlashAttention's tiling, the bigger cousin of this kernel.
• Module 08 · Inference-aware design: KV-cache maths and why decode reads so many bytes.
• Page 06 · Quantisation and choosing GPUs: fewer bytes per weight moves decode up the memory slope.
• Page 05 · Latency: inter-token latency and tail latency, which profiling explains.
• Module 04 · Latency as a safety number: measuring wait-after-speech on the device.
• Module 08 · Course map: GPU programming resources.

Why it matters for robotics

Self-check

Resources

• Nsight Systems user guide (NVIDIA): timelines, NVTX, CUDA graph tracing
• Nsight Compute profiling guide (NVIDIA): speed-of-light sections, roofline charts, warp stall reasons
• DCGM feature overview (NVIDIA): profiling metrics beyond nvidia-smi
• CUDA C++ programming guide (NVIDIA): CUDA graphs, streams, the memory hierarchy
• NVIDIA H100 product page: the specs used on this page
• Triton tutorial: layer normalisation, the closest official example to RMSNorm
• triton.testing.do_bench reference
• torch.nn.RMSNorm and the PyTorch profiler
• Root Mean Square Layer Normalization (Zhang and Sennrich, 2019)
• FlashAttention (Dao et al., 2022), FlashAttention-2 (Dao, 2023), FlashAttention-3 (Shah et al., 2024): fusion and IO-awareness at full scale
• vLLM documentation: the engine for the utilisation exercise