Five Self-Contained CUDA Optimizations for 3D Gaussian Splatting

George Kopanas·Google DeepMind, San Francisco

TL;DR. 7.3× faster forward and backward kernels at 30k (8.4× at 7k) over Inria [1] and gsplat [4], broken down in 6 isolated CUDA optimization steps you can port ad hoc into any splat rasterizer and training harness.

Since its introduction by Kerbl & Kopanas et al. (SIGGRAPH 2023) [1], 3D Gaussian Splatting (3DGS) has become one of the most heavily optimized GPU pipelines in computer graphics and vision. I still remember, even before we started working on the original 3DGS paper, writing my first differentiable software rasterizers: a single forward pass took 500 ms and could barely handle 1.5 million points. I was doing my best, and in every weekly meeting my supervisor, George Drettakis, would insist: “I don't know how to do it, that's for you to figure out, but there is no way this is the fastest a GPU can render 1.5 million points.” His intuition gave him the conviction that kept us going, and I only realized how invaluable that was later, when that early code was transformed with the help of Bernhard Kerbl into the published kernels of the Inria implementation. When we released those kernels, it was hard to imagine that within three years they could be made another ~7× faster, as we show in this post. Yet what surprised me even more is that even after combining all the optimizations proposed since then by students, academics, and engineers across the community, we are still not hitting the hardware ceiling—another 1.78× speedup on the full forward-backward training step remains on the table. In this post, we explore where that remaining headroom comes from and how to take advantage of it. I wrote this for graphics and machine-learning practitioners, whether you are working on 3DGS itself or want a concrete, end-to-end case study of how to optimize a non-trivial GPU pipeline from first principles. I assume basic prior knowledge of 3DGS and some familiarity with the technical details of the original CUDA implementation, while explaining the GPU hardware concepts we need right before we use them.

To establish a rigorous literature baseline, we consolidate five families of optimizations from recent publications and open-source repositories into a single package, s01_literature, drawing on DISTWAR [2], Taming 3DGS [3], gsplat [4], StopThePop [5], FlashGS [6], Speedy-Splat [7], Splatshop [8], Faster-GS [9], HiGS [10], Yang et al. [11], and LichtFeld Studio (fastgs) [12]. Together, they deliver a 4.09× speedup at 30k (4.53× at 7k) over Inria across the 9 Mip-NeRF 360 scenes. This consolidated baseline is already faster than any existing open-source implementation we could find, outperforming the standalone LichtFeld Studio (fastgs) [12] rasterizer by 1.67× at 30k (1.63× at 7k).

Where does the remaining 1.78× speedup come from? Once those literature techniques are in place, profiling reveals a second tier of bottlenecks rooted in how the kernels map onto the GPU's warp execution model, register file, L2 cache sectors, and arithmetic execution pipes. In the rest of this post, we will first recap and explain all the optimizations we incorporated from the literature (s01_literature), and then walk through five self-contained optimizations (s02–s06) that address these hardware bottlenecks one by one. Each step lives in its own standalone package with a direct code diff against the previous step, reaching our 7.30× end-to-end speedup at 30k (8.38× at 7k) over Inria [1] and gsplat [4] (and 2.97× over LichtFeld Studio [12]) with identical image quality.

Collapsible section (click to expandcollapse) — SM sub-partitions, warp schedulers, FMA/ALU/MUFU/LSU execution pipes, L1 and shared memory, 32-byte L2 cache sectors, write-through stores, and ROP atomics.

40 GB HBM2 Off-Chip Memory (1,555 GB/s Peak · ~400–600 Cycles Latency · 32-Byte Sector Granularity) 5.83M Gaussians × 236 Bytes = 1.38 GB of scene parameters (~0.89 ms minimum just to read once at peak HBM2 bandwidth)

40 MB Shared L2 Cache (Split into 80 Parallel Slices · 128-Byte Cache Lines = 4 × 32-Byte Sectors) Global writes skip L1 and go straight to L2 in 32-byte sectors; atomicAdds to the same address queue up at one slice

Inside One of the 108 Streaming Multiprocessors (SMs) · Up to 64 Resident Warps (2,048 Threads)

Sub-Partition 0 1 warp scheduler issues 1 warp instruction / cycle

16,384 × 32-bit registers (64 KB) 16 FP32 · 16 INT32 · 8 LSU · 4 MUFU units

Sub-Partition 1 1 warp scheduler issues 1 warp instruction / cycle

16,384 × 32-bit registers (64 KB) 16 FP32 · 16 INT32 · 8 LSU · 4 MUFU units

Sub-Partition 2 1 warp scheduler issues 1 warp instruction / cycle

16,384 × 32-bit registers (64 KB) 16 FP32 · 16 INT32 · 8 LSU · 4 MUFU units

Sub-Partition 3 1 warp scheduler issues 1 warp instruction / cycle

16,384 × 32-bit registers (64 KB) 16 FP32 · 16 INT32 · 8 LSU · 4 MUFU units

SM-Wide Pipe Throughput (4 Sub-Partitions Combined)

FMA Pipe FP32 add/mul/fma 64 threads / cycle = 2 warp instr / cycle

ALU Pipe INT32, cmp, sel 64 threads / cycle = 2 warp instr / cycle

MUFU Pipe ex2, rcp, rsqrt 16 threads / cycle = 0.5 warp instr / cycle

LSU & Shuffle ld/st, _shfl* 32 threads / cycle = 1 warp instr / cycle

192 KB On-Chip SRAM (~20–30 Cycles)

L1 data cache hardware-managed caches global loads

Shared memory programmer-managed shared, up to 164 KB one SRAM, split configurably · 32 banks × 4 B / cycle

Figure 1: The memory hierarchy and compute layout of the NVIDIA GA100 GPU in the A100-SXM4-40GB. Top: 40 GB of HBM2 feed a 40 MB L2 cache, split into 80 slices, that all 108 SMs share; global atomics execute inside the L2, and global writes travel to it in 32-byte sectors. Bottom: the inside of a single SM. Each of its four sub-partitions has its own warp scheduler, which issues one warp instruction per cycle, and its own 64 KB register file. The lower row adds the four sub-partitions together into the throughput of each pipe for the whole SM, next to the 192 KB of on-chip SRAM that the SM splits between L1 cache and shared memory. Grey boxes hold data, white boxes compute; the slate line marks the atomic serialization that dominates the baseline.

Every performance cliff we encounter in 3D Gaussian Splatting, and every speedup we achieve across Steps 1 through 6, traces directly back to how the NVIDIA GA100 (A100-SXM4-40GB) schedules instructions inside its 108 Streaming Multiprocessors (SMs) and moves data through its two-level cache hierarchy.

An A100 GPU contains 108 SMs clocked at up to 1.41 GHz. Each SM contains four independent processing blocks (sub-partitions), each with its own warp scheduler and register file (65,536 32-bit registers per SM in total). The hardware partitions the threads of a thread block into 32-thread warps that execute a single instruction stream in Single-Instruction, Multiple-Thread (SIMT) lockstep. Throughout this post we denote the 32 threads of a warp by T0..T31. Because each SM can hold up to 64 resident warps (2,048 threads) but only 65,536 registers, the number of 32-bit registers allocated per thread directly governs how many warps can reside on an SM at once to hide memory and instruction latency: only a kernel using at most 32 registers per thread can sustain full occupancy (64 warps/SM), whereas holding large local arrays per thread forces the compiler (ptxas) either to allocate more registers, throttling warp residency, or to spill excess live values into slow local memory in L1/L2.

Each of the four sub-partitions has its own warp scheduler, and every cycle each scheduler can issue one warp instruction, that is, one operation for all 32 threads T0..T31. The SM as a whole therefore issues up to four warp instructions per cycle. Each instruction goes to one of several separate functional units, or pipes: the FMA pipe for floating-point math, the ALU pipe for integer arithmetic, comparisons, and predicate logic, the MUFU pipe for transcendentals such as exponentials and reciprocals, and the LSU pipe for memory accesses and cross-thread shuffles. How fast a pipe drains those instructions depends on how many threads it can execute per cycle. A sub-partition has 16 FP32 units, so one 32-thread FMA instruction occupies its FMA pipe for two cycles; summed over the four sub-partitions, the SM completes 64 FP32 thread operations, or two warp instructions, per cycle, and the ALU pipe has the same width. The LSU pipe is half as wide, one warp instruction per cycle per SM, and the MUFU pipe a quarter, half a warp instruction per cycle. Four instructions per cycle is thus a ceiling the SM reaches only when it spreads work across several pipes. Because the pipes are independent, a loop made mostly of floating-point math bottlenecks on the FMA pipe while the ALU pipe idles, and a loop made mostly of integer indexing, bitmask manipulation, and predicate tests bottlenecks on the ALU pipe no matter how much floating-point work we remove.

Branches interact with this lockstep execution in two ways, and Step 5 relies on both. First, all 32 threads T0..T31 of a warp share a single program counter, so a branch only saves work when all 32 threads take the same side. When 8 pixels of a 32-pixel warp fall inside a small Gaussian ellipse and 24 fall outside, an early-exit branch such as if (!hit) continue; does not make the warp any faster: the 8 active threads still step through every instruction of the loop body while the other 24 wait. Second, the branch costs time even when no thread diverges. Its mere presence forces the compiler to wrap the loop body in extra control-flow instructions: BSSY before the branch to record where diverged threads will later rejoin, BRA for the conditional jump itself, and BSYNC at the end of the loop body to re-synchronize the warp. Every warp executes these instructions on every iteration, whether its threads split or not. They also cut one straight-line loop body into several short basic blocks, and the compiler cannot move independent loads and math across those block boundaries to hide their latency.

Next to its registers, each SM has 192 KB of on-chip SRAM (~20–30 cycles away) split between the hardware-managed L1 data cache and up to 164 KB of programmer-managed shared memory (shared, organized into 32 four-byte banks delivering 128 B/cycle). Shared memory is what makes tile-based rasterization fast: instead of every pixel thread fetching each Gaussian from L2, the 256 threads of a 16 × 16 tile cooperatively stage a batch of Gaussians into shared memory once, and when all 32 threads of a warp read the same Gaussian, the hardware broadcasts the shared-memory word in a single conflict-free access.

Outside the SMs sits a 40 MB L2 cache (~200 cycles away, backed by 1,555 GB/s HBM2) partitioned into 80 parallel hardware slices. Each slice stores 128-byte cache lines subdivided into four 32-byte sectors (8 floats)—the minimum unit transferred between an SM and L2. While global loads cache in L1, global stores are write-through: every warp store immediately sends each modified 32-byte sector to L2. When T0..T31 write 32 adjacent floats (128 B), they fill four complete sectors at 100% efficiency; when they write one float per Gaussian at a 192-byte stride, each 4-byte store touches a separate 32-byte sector, sending 8× more L2 traffic than payload.

Finally, global floating-point atomicAdd instructions execute inside the Raster Operation (ROP) unit of each L2 slice. Throughput depends on whether atomics hit the same 32-bit address: atomics to different cache lines spread across the 80 slices, and atomics to distinct floats within the same cache line pipeline without dependency stalls (as in Step 2, where 18 threads update 18 distinct gradient floats at once). Only atomics to the exact same 32-bit address serialize strictly one after another—which is why the baseline rasterizer stalls for hundreds of cycle