What the 3DGS Literature Left on the Table

Step 1: Consolidating What the Literature Already Knowss01_literature·diff vs. s00_baseline ↗

Before claiming any new optimization, we want to establish an uncompromising baseline that consolidates every major algorithmic and systems improvement already known in the 3DGS literature into a single, cohesive CUDA implementation. While this step does not introduce any new technique of its own, it is the first open-source implementation to combine every trick from the literature under one roof, achieving state-of-the-art performance that is already 1.67× faster than LichtFeld Studio (fastgs) [12] at 30k (and 1.63× faster at 7k).

(a) Warp-butterfly backward gradient reduction (DISTWAR [2], Taming 3DGS [3], gsplat [4], LichtFeld Studio [12]). Because all 32 threads T0..T31 of a warp step through the same shared-memory splat simultaneously and hold their 9 per-pixel screen-space derivatives gr[0..8] in registers (2D center μ′_x, μ′_y, conic a, b, c, opacity o, and RGB color r, g, b), we sum each splat's gradients across the 32 threads of the warp in registers before touching global memory. In CUDA, threads in the same warp can exchange a 32-bit register directly with a partner thread in a single hardware instruction—a warp shuffle (__shfl_xor_sync)—without touching shared or global memory. Using a 5-step butterfly reductionWarp Butterfly Reduction (__shfl_xor_sync)Sums a value across all 32 threads of a warp in log₂(32) = 5 register steps. At each step offset ∈ {16, 8, 4, 2, 1}, every thread swaps its running sum with partner thread ^ offset and adds the two:Named after the crossing "butterfly" edges of an FFT network. (5 shuffles per float with offset = 16, 8, 4, 2, 1, or 9 × 5 = 45 shuffles per splat), thread 0 accumulates the warp's 32-pixel subtotal and issues 9 global atomicAdds, reducing L2 atomic traffic by up to 32×.

(b) Exact 2D conic-tile intersection & opacity-aware cutoffs (StopThePop [5], FlashGS [6], Speedy-Splat [7], LichtFeld Studio [12]). Kerbl & Kopanas et al.'s 3𝜎 circle radius ignores both the anisotropy and the opacity o of a splat, emitting every tile in a conservative square around the Gaussian (Figure 3, left). Since a splat's contribution at pixel offset d is 𝛼 = o e^−q(d)/2 with q(d) = d^⊤Λd, it can only exceed the 1/255 blending cutoff inside the opacity-aware ellipse q(d) ≤ 2 ln(255 o). A first improvement (StopThePop [5], FlashGS [6]) shrinks the search region to the tight axis-aligned bounding box (AABB) of that ellipse and tests each tile in the AABB in a 2D loop, though tilted ellipses still waste loop iterations evaluating and rejecting empty corner tiles inside the AABB (Figure 3, middle). Following Speedy-Splat [7] and LichtFeld Studio [12], s01_literature instead solves the quadratic boundary Δx(Δy) in closed form for each covered tile row ty to obtain that row's exact horizontal tile span [tx0, tx1) directly, evaluating zero rejected tiles (Figure 3, right). On room 30k view 0, this cuts emitted tile-Gaussian pairs by 73.8% (3.82× fewer pairs, 7,004,474 → 1,832,356).

(c) Per-warp 8×4 sub-tile culling bitmask & compact 48-byte splat layout (HiGS [10], Yang et al. [11], FlashGS [6]). Each 16 × 16 pixel tile is processed by 8 warps (W0..W7), where each warp is responsible for an 8×4 rectangular sub-tile of 32 pixels (Figure 4, left). Many splats that intersect the 16 × 16 tile only touch one or two of those eight sub-tiles; if all 8 warps blindly loop through every splat in the tile's list, warps whose 8×4 pixels lie completely outside a splat waste cycles evaluating it. To avoid this, whenever a batch of splats is staged into shared memory, we apply the same analytic 1D row-span algorithm from (b) at 8×4 sub-tile granularity (WarpMask<1>) to compute an 8-bit mask (batch_mask) for each splat (one bit per warp) recording which of the eight 8×4 sub-tiles the splat's ellipse overlaps. Before blending a 32-splat chunk, each warp uses a single hardware warp-vote instruction (__ballot_sync) to gather its sub-tile bit across all 32 staged splats into a 32-bit integer worklist, and then uses find-first-set (__ffs) to jump directly from one overlapping splat to the next (Figure 4, right), skipping splats that miss its sub-tile in zero loop iterations. In addition, all 2D attributes of each visible Gaussian (center, conic, RGB color, and opacity) are packed contiguously into a single 16-byte-aligned 48-byte struct holding three float4 fields, allowing threads to stage each splat into shared memory with three wide 16-byte float4 reads instead of scattered scalar loads.

(d) Two-stage visibility-indexed 32-bit / 16-bit radix sorting (Splatshop [8], Faster-GS [9]). Instead of running an 8-pass 64-bit CUB radix sort over all I (tile | depth) keys, s01_literature implements the two-stage sorting pipeline of Splatshop [8] and Faster-GS [9] (Figure 5): first, filter the M visible Gaussians (M ≪ N) and sort them front-to-back once by their 32-bit float depth key; second, emit 16-bit tile keys in that depth-sorted order and run a narrow 2-pass stable radix sort over the 16-bit tile IDs. On room 30k view 0, sorting time drops by 79% (0.97 ms → 0.20 ms across the two radix sorts).

(e) Visibility-gated SH evaluation, wide float4 SH loads & base-2 MUFU exponential (Taming 3DGS [3], gsplat [4], FlashGS [6], Faster-GS [9], LichtFeld Studio [12]). Finally, s01_literature gates forward Spherical Harmonics evaluation behind visibility so off-screen Gaussians never touch their 192-byte SH table, replaces scalar 4-byte SH reads with wide 16-byte float4 vector loads, pre-scales conic coefficients by −½ log₂(e) so the inner rasterization loop evaluates 1-cycle hardware base-2 exponentials (ex2.approx.ftz.f32 on the SM's MUFU pipe), and replaces the 3-component running color accumulator in raster_bwd with a single scalar projection accum = c_behind · ∇I onto the pixel's upstream color gradient (saving 4 FMAs per pixel–splat evaluation) [3, 4, 12].

Comparison against LichtFeld Studio (fastgs) [12]. Combining all five literature pillars cuts the stage total on room 30k view 0 from 25.06 ms down to 4.86 ms (5.2×), and the end-to-end synchronous step from 25.18 ms to 4.77 ms, a 5.5× speedup over Inria's 26.09 ms. Across all 9 Mip-NeRF 360 scenes at 30k, s01_literature averages 5.87 ms (4.09× vs. Inria). Notably, s01_literature is already 1.67× faster at 30k and 1.63× faster at 7k than the standalone open-source LichtFeld Studio (fastgs) rasterizer [12], which shares much of the same algorithmic foundations ([3, 5, 7, 9]) but leaves per-tile batch_mask bitmask filtering in forward from HiGS [10] and Yang et al. [11] (c), the single packed 48-byte Splat struct from FlashGS [6] (c) (it stages conic/opacity and color from separate float4 arrays), and the two-stage depth / 16-bit tile sort from Splatshop [8] and Faster-GS [9] (d) unexploited, sorting a single packed (tile | depth) key instead. Everything from Step 2 (s02) through Step 6 (s06) is therefore strictly beyond the published literature.

Step 2: Matching the Reduction to the Gradient Width — Warp Reduce-Scatter & Two-Splat Foldings02_reduce_scatter·diff vs. s01_literature ↗

After Step 1, raster_bwd is still the largest stage on every scene (2.90 ms on room 30k, or 60% of the step; 3.98 ms on garden 30k, or 45% of the step). Although warp-level aggregation in s01_literature (following DISTWAR [2], Taming 3DGS [3], and LichtFeld Studio [12]) eliminated per-pixel L2 atomic collisions within each warp, inspect how every published warp reduction computes that sum (Figure 6, left): each thread holds 9 gradient floats for the current splat, and the kernel reduces them one component at a time using 9 full 32-thread butterfly reductionsWarp Butterfly Reduction (__shfl_xor_sync)Sums a value across all 32 threads of a warp in log₂(32) = 5 register steps. At each step offset ∈ {16, 8, 4, 2, 1}, every thread swaps its running sum with partner thread ^ offset and adds the two:Named after the crossing "butterfly" edges of an FFT network.. Each step of a butterfly reduction is a warp shuffle (__shfl_xor_sync)—a hardware instruction in which every thread exchanges one 32-bit float register directly with a partner thread in the same warp—so reducing 9 floats across 5 steps takes 9 × 5 = 45 shuffles per splat (or 30 shuffles when 6 intermediate moments are reduced before reconstructing the 9 outputs).

That naive loop (Figure 6, left) wastes work in two ways. First, at the end of a 5-step butterfly, all 32 threads hold a duplicate copy of the sum even though 31 of the 32 threads immediately throw their copy away, so repeating that across all 9 components costs 9 × 5 = 45 shuffle instructions per splat. Second, a warp only sums 32 pixels, whereas a 16 × 16 tile contains 8 warps and a splat can span multiple tiles running on different SMs, so each warp must still atomically add its 32-pixel subtotal into the splat's global gradient array in L2. Because s01_literature leaves all 9 warp subtotals on thread T0, T0 has to issue 9 global atomicAdd instructions one after another while threads T1..T31 sit idle.

Halving Warp Reduce-Scatter (45 → 12 Shuffles). As Figure 6 (right) illustrates, instead of keeping all 9 gradient components alive on every thread through all 5 steps, we can split the 9 components across the threads as we sum them. At each step, paired threads divide their remaining components in half—one partner receives and adds the lower half from its partner while the other receives and adds the upper half. Because the number of active components per thread halves at every step (9 → 5 → 3 → 2 → 1 → 1), each thread issues only 5 + 3 + 2 + 1 + 1 = 12 register shuffles (__shfl_xor_sync exchanges) instead of 45. Moreover, at the end of Step 5, 9 distinct threads in the warp each hold the 32-pixel warp sum for one of the 9 gradient components. Instead of thread T0 looping over 9 sequential global atomicAdds, all 9 threads write their warp subtotal into the splat's global gradient array in L2 at the same time in a single warp instruction.

Reducing Two Splats Together in Half-Warps (12 → 10 Shuffles per Splat). In Figure 6 (right), our toy example had 4 gradient components on 4 threads, so the active count halved cleanly all the way down (4 → 2 → 1) and every thread ended up with one sum. When the warp has more threads than gradient components, however—as in 3DGS, where a warp has 32 threads but a splat has only 9 gradients (so even a 16-thread half-warp is wide enough to hold all 9 sums)—the active count per thread hits 1 before the last step (9 → 5 → 3 → 2 → 1 → 1) and cannot be split any further. Figure 6b (left) shows this on our 4-thread toy warp with a 2-component splat (a0, a1): Step 1 splits 2 → 1 (1 shuffle), and then Step 2 still has to spend 1 more shuffle to sum that single component across pairs (1 → 1), leaving half the threads idle (2 shuffles per splat, or 4 shuffles for two splats). Instead, we can evaluate two consecutive splats (A and B) before reducing and use Step 1 to fold splat A into the lower half-warp and splat B into the upper half-warp, so both half-warps then run the halving reduce-scatter in parallel without wasted tail steps (Figure 6b, right).

Inspecting the compiled SM80 SASS of RasterBackwardKernel with cuobjdump --dump-sass confirms the structural change: static SHFL.BFLY instructions drop from 366 to 230 (−37%) and global atomic reduction instructions (RED.E.ADD.F32.FTZ.RN) drop from 72 to 16 (−78%). Because the halving reduction keeps more partial sums live per thread, we also give RasterBackwardKernel a spill-free register budget (__launch_bounds__(256, 4), up to 64 registers instead of 40), and, as in s01_literature, still skip the reduction entirely when no pixel of the warp blends the splat pair. On room 30k, raster_bwd drops from 2.90 ms to 2.08 ms (−28%), and on garden 30k from 3.98 ms to 3.46 ms (−13%).

Step 3: Eliminating SH Register Pressure & Bulk Memsets — 4-Group Register-Window Streaming & Cooperative Warp Zeroings03_sh_streaming·diff vs. s02_reduce_scatter ↗

The runtime of preprocess_bwd scales with the total number of Gaussians N in the scene rather than with the number of pixels: on room, with 1.59M Gaussians, it takes 0.54 ms, whereas on garden, with 5.83M Gaussians, it takes 2.18 ms—about a quarter of the entire step and the second-largest stage after raster_bwd. For this reason Step 3 reports garden 30k numbers. At first sight, 2.18 ms is surprising: preprocess_bwd launches only one thread per Gaussian and performs simple closed-form polynomial math. Why is such a simple kernel so slow?

The bottleneck is memory bandwidth and register pressure, not arithmetic. In 3DGS, each Gaussian stores 16 RGB Spherical Harmonics (SH) coefficients—16 × 3 = 48 contiguous floats (192 bytes per Gaussian)—laid out back-to-back in global memory so that Gaussian G0 occupies bytes 0..191, Gaussian G1 occupies bytes 192..383, and so on. Even after s01_literature gated forward SH evaluation behind visibility and widened per-thread SH loads to 16-byte float4 vectors, two bottlenecks still remain in preprocess_bwd.

96 live SH registers per thread. In preprocess_bwd, each active thread loads all 48 forward SH floats into a local array coeffs[48] and accumulates all 48 SH gradients into a second array dsh[48] before writing any of them out—tying up 96 registers per thread for SH alone, which totals 152 registers per thread and caps SM occupancy on A100 at 18.75% (12 warps per SM with 128-thread blocks).

To cut register pressure, we exploit the fact that the 16 SH basis functions Y_k(d) are additive and independent across basis index k. In the forward pass, the view-dependent RGB color ∑_k Y_k(d) s_k only requires a 3-float running sum. In the backward pass, each coefficient's gradient ∂L/∂s_k = Y_k(d) (∂L/∂c) depends only on the incoming RGB gradient ∂L/∂c and the viewing direction d, while the viewing-direction derivative ∂L/∂d = ∑_k (∂L/∂c · s_k) ∇_dY_k(d) is likewise a 3-float running sum. We therefore never need all 48 forward coefficients or all 48 gradients in registers at the same time. Instead, we rearrange the SH evaluation and VJP into four passes of 4 SH basis functions (4 × 3 = 12 floats) each, streaming one 12-float group at a time through a single register buffer v[12]—the smallest multiple of 3 RGB channels that packs evenly into 4-wide float4 (16-byte) loads, the widest load instruction a single GPU thread can issue. Streaming SH through a single 12-float register window cuts PreprocessBackwardKernel register usage by 32% (from 152 to 104 registers/thread in float32, and down to 94 when combined with bfloat16 SH in Step 6), raising resident warp occupancy on A100 from 18.75% to 25% (12 → 16 warps per SM).

Redundant 1.38 GB bulk cudaMemsetAsync & Cooperative Warp-Coalesced Zeroing of Untouched SH Rows. Before the backward pass runs, a bulk cudaMemsetAsync call zeroes all 59 gradient floats of all N = 5.83M Gaussians (1.38 GB on garden 30k) so that culled or untouched Gaussians end up with the correct zero gradients. That writes memory twice for every visible Gaussian: first zeroed by cudaMemsetAsync, then overwritten by preprocess_bwd. Yet preprocess_bwd already launches a thread for every Gaussian 0..N−1: if the kernel itself writes zeros for the inactive Gaussians while writing the computed gradients for the active ones, the separate 1.38 GB cudaMemsetAsync pass can be removed entirely.

Because one thread handles one Gaussian, if each inactive thread simply zeroed its own 192-byte SH row in a per-thread loop, adjacent threads in a warp would write 192 bytes apart in global memory (Figure 7, top). Recall from the hardware background section that the A100 L1/L2 caches transfer memory in 32-byte sectors: because adjacent threads are 192 bytes apart instead of 16 bytes apart, no two threads share a sector—every thread touches a separate 32-byte L2 sector. Instead, because preprocess_bwd assigns consecutive threads T0..T31 of a warp to 32 consecutive Gaussians G0..G31, and the SH tensor stores those Gaussians back-to-back in global memory, the warp's 32 SH rows form one contiguous 6,144-byte block. The warp first votes on activity with a single __ballot_sync instruction so every thread knows which of the warp's 32 Gaussians were untouched, and then the 32 threads march side-by-side across that contiguous block 32 chunks at a time (Figure 7, bottom), writing 16 bytes of zeros whenever a chunk belongs to an untouched Gaussian and skipping active Gaussians.

Because adjacent threads T0..T31 always address adjacent 16-byte chunks, every zeroing step coalesces cleanly into complete 32-byte L2 sectors without using a single byte of shared memory. Together, these two changes cut preprocess_bwd on garden 30k from 2.18 ms to 1.81 ms, and because removing the 5 bulk host-launched cudaMemsetAsync passes also eliminates inter-kernel launch overhead, the end-to-end synchronous fwd-bwd step on garden 30k drops from 7.62 ms to 7.38 ms.

Step 4: Pruning the Backward Pass — In-Place Contributor Compaction, Upper-Byte cmask Packing & splat_touched Gatings04_raster_pruning·diff vs. s03_sh_streaming ↗

Even with exact ellipse-tile intersection and 8 × 4 sub-tile bitmasks (Figure 4) in s01–s03, a large fraction of the splats assigned to a tile never actually blend into a single pixel because front-to-back transmittance T saturates below 10⁻⁴ on the pixels underneath the splat before the forward pass reaches its depth. Yet in s01–s03, as in prior tile rasterizers [1, 2, 3, 4, 6, 12], raster_bwd walks the tile's unpruned splat list in reverse from the last contributing splat back to the front, staging every occluded splat in that prefix into shared memory and recomputing its sub-tile mask, while preprocess_bwd runs on every frustum-visible Gaussian even if it was completely occluded across all tiles. Step 4 eliminates this wasted backward work through three co-designed mechanisms illustrated in Figure 8:

1. In-Place Forward Contributor Compaction. As each warp steps through a 32-splat chunk during raster_fwd (Figure 8, top row), it records a 32-bit bitmask of the splats that actually blended with 𝛼 ≥ 1/255 into at least one pixel of that warp and writes the mask to shared memory. At the end of each 256-splat batch, taking the bitwise OR of those masks across all warps identifies every splat that contributed anywhere in the tile (such as S0, S2, and S3 in Figure 8, while dropping occluded splats like S1, S4, and S5). A short prefix sum over the population counts of those masks gives each contributing splat its dense destination index, allowing the thread block to overwrite the front of the tile's sorted splat list in place with only the true contributors (Figure 8, middle row) and remap each pixel's stopping index into the compacted stream. On room 30k view 0, this prunes 20.7% of all tile-Gaussian pairs before the backward pass begins, shrinking the tile splat lists from 1.83M to 1.45M entries.

2. Upper-Byte cmask Packing. Because Mip-NeRF 360 Gaussian indices fit comfortably within 24 bits for scenes up to 16.7M Gaussians, forward compaction also gathers the 8 per-warp contribution bits for each surviving splat into an 8-bit contributor mask cmask and packs it into the upper byte of the compacted 32-bit splat-list entry above the 24-bit Gaussian ID (Figure 8, bottom row). When raster_bwd walks the compacted prefix in reverse, each staging thread unpacks the 24-bit Gaussian ID and the 8-bit warp mask directly from that 32-bit word—eliminating sub-tile bounding-box evaluation from the backward pass altogether and skipping warps where a splat overlapped the sub-tile ellipse in geometry but was occluded on all pixels of that warp.

3. 1-Byte splat_touched Backward Preprocess Gating. As raster_bwd stages each compacted contributing entry (Figure 8, middle row, touched[Si]=1), it writes a single byte flag marking that Gaussian ID as having contributed to the rendered image. In preprocess_bwd, threads check this flag instead of frustum visibility, skipping both the 9-float intermediate gradient load and the 48-coefficient Spherical Harmonics backward pass for every frustum-visible Gaussian that was occluded across all tiles. On room 30k view 0, spending an extra 0.16 ms in raster_fwd to compact the splat list cuts raster_bwd by 18.2%, from 2.09 ms to 1.71 ms, and preprocess_bwd from 0.46 ms to 0.42 ms, reducing total stage time from 3.96 ms to 3.71 ms and bringing the end-to-end synchronous step to 3.70 ms.

Step 5: Eliminating SIMT Warp Divergence — Branchless Predication & Exact Quadratic Roundings05_branchless_predication·diff vs. s04_raster_pruning ↗

After Step 4's pruning, the inner per-splat loops of raster_fwd (0.98 ms) and raster_bwd (1.71 ms) account for 73% of the total step time (3.71 ms) on room 30k. Crucially, the per-warp cmask sub-tile bitmask and Step 4's contributor compaction guarantee an important structural invariant: every splat that a warp's inner loop visits already intersects at least one pixel in that warp's sub-tile. Yet in s01–s04, the per-splat loop inside both rasterizers still contains divergent per-thread early-exit branches that test whether the pixel is still active and whether the evaluated opacity satisfies 𝛼 ≥ 1/255.

Why are these branches counterproductive once sub-tile culling is in place? Because all 32 threads of a warp share a single hardware instruction pointer (SIMT), a branch can only skip the rest of the loop body if all 32 threads take the early exit together. Because the warp's sub-tile bitmask already skipped every splat that missed the whole warp, at least one thread in the warp hits the current splat and stays in the loop body. Consequently, the warp must execute every blending and gradient instruction of the loop body anyway, while threads that took the early exit sit idle with their predicate masks disabled. Worse, as shown in Figure 9 (left), the compiler surrounds those divergent branches with hardware branch-reconvergence instructions (BSSY before the branch and BSYNC where threads rejoin) that chop the inner loop into multiple basic blocks and stall instruction issue.

Branchless 𝛼 = 0 Predication. Instead of branching out of the loop when a pixel misses a splat, s05_branchless_predication keeps all 32 threads converged through the entire inner loop and uses a single-cycle hardware select instruction (SEL) to set 𝛼 = 0 on missed pixels (Figure 9, right). Because a zero-opacity splat adds zero color (𝛼Tc = 0) and multiplies transmittance by 1 − 𝛼 = 1, missed threads can execute the exact same forward blending and backward gradient instructions as active threads without changing their pixel's running color or transmittance, and their zero derivatives contribute nothing to the warp reduction.

Exact Round-to-Nearest Quadratic Evaluation. Running the backward update unconditionally on every thread exposes a subtle numerical requirement: the forward and backward kernels must agree bit-for-bit on whether a pixel crossed the 𝛼 ≥ 1/255 threshold. If the compiler fuses floating-point multiply-adds differently between the two kernels, a pixel sitting right on the boundary can evaluate as a hit in the forward pass and a miss in the backward pass, leaving its transmittance T un-divided by 1 − 𝛼. To prevent this mismatch, we evaluate the 2D quadratic form in both kernels using explicit IEEE round-to-nearest-even instructions, guaranteeing bit-identical exponents and opacity values in forward and backward. Inspecting the compiled SM80 assembly confirms that branch-reconvergence instruction pairs in raster_bwd drop from 66 to 42. On room 30k, raster_fwd drops from 0.98 ms to 0.78 ms and raster_bwd drops from 1.71 ms to 1.35 ms, cutting the stage total from 3.71 ms to 3.15 ms.

Step 6: Balancing SM Execution Pipes — Multi-Row Warp Mapping & Pipe-Balanced Maths06_multirow_and_pipe_math·diff vs. s05_branchless_predication ↗

With branchless predication in place, both raster_fwd (0.78 ms) and raster_bwd (1.35 ms) execute straight-line inner loops without L2 atomic stalls or SIMT reconvergence penalties. Profiling s05_branchless_predication with Nsight Compute (on a synthetic 6M-Gaussian scene calibrated to garden 30k; all hardware counters below come from this scene) reveals the final two hardware bottlenecks inside the SM sub-partitions. First, in raster_bwd, warps spend 14.9% of their cycles stalled at block-wide __syncthreads() barriers, waiting for the slowest warp in the 16 × 16 tile to finish its batch before the block can stage the next shared-memory batch, while inside the loop body itself, the FMA pipe (floating-point multiply-add) is the busiest execution unit, running at 68.9% of its peak throughput (compared to 40.0% for the integer/comparison ALU pipe and 38.9% for the load/store LSU pipe). Second, in raster_fwd, no single execution pipe is saturated, yet the warp scheduler is: how can the scheduler be saturated when every execution pipe is well below 100%? Recall from Figure 1 that each of the four SM sub-partitions has a single warp scheduler that can issue at most 1 instruction per cycle in total across all of its functional units combined (FMA, ALU, MUFU, and LSU). Because the forward loop alternates almost evenly between integer/predicate logic (ALU, 201.9M instructions/frame), floating-point math (FMA, 205.4M instructions/frame), and fast transcendentals (XU/MUFU ex2/rcp, 54.4M instructions/frame), the warp scheduler hits its 1-instruction-per-cycle issue ceiling first—keeping issue slots busy 81.2% of active cycles—while the ALU, FMA, and XU pipes each sit at 60–65% of their individual peak throughputs. The forward loop is therefore instruction-issue bound: removing floating-point work alone cannot speed it up; the total number of issued instructions per pixel–splat evaluation, across ALU, FMA, and XU, has to come down. Step 6 (s06_multirow_and_pipe_math) resolves both bottlenecks through five co-designed transformations.

1. Square 8 × 8 Multi-Row Thread Mapping (2 Vertical Pixels per Thread). In s00–s05, a 16 × 16 tile block launched 256 threads—8 warps of 8 × 4 pixels (Figure 10, left)—so every shared-memory load, sub-tile bitmask scan, horizontal Δx setup, and backward warp reduction ran 8 times per tile. In s06_multirow_and_pipe_math, each thread instead processes two vertically adjacent pixels in the same column, so a 32-thread warp covers a square 8 × 8 quadrant (64 pixels) and the entire 16 × 16 tile needs only 128 threads = 4 warps (Figure 10, right). Because both pixels in a thread share the same horizontal coordinate, the horizontal offset Δx = p_x − μ′_x and the horizontal quadratic term depend only on the column and are computed once per thread rather than twice. Furthermore, in raster_bwd, each thread sums the partial derivatives of its two vertical pixels in registers before the warp reduction begins. Summing the two rows in registers first reduces the 256 pixels of the tile down to 128 thread totals across 4 warps instead of 8, cutting the number of warp reductions and global atomicAdd writes per tile in half (4 warps per tile instead of 8).

2. Barrier-Free Per-Warp Backward Traversal (__syncwarp()). What about the 14.9% of cycles that raster_bwd spent waiting at block-wide __syncthreads() barriers? In s00–s05, all warps in the 16 × 16 tile read splats from a single shared-memory buffer, so no warp could load the next batch until every warp in the tile had finished the current batch; fast warps whose quadrant of the tile was mostly empty had to sit idle waiting for the busiest warp. Instead, we partition shared memory in raster_bwd into 4 independent 32-splat staging slices—one private slice for each of the 4 warps—and synchronize each warp only with its own 32 threads via __syncwarp(). Using the 8-bit cmask packed into the upper byte of each entry in Step 4, each warp loads and processes only the splats that contributed to its own 8 × 8 quadrant, marching backward from its last contributor to the front completely independently of the other three warps in the tile and reducing barrier stalls from 14.9% to 0.0%.

3. Folding log₂(opacity) into the Column Constant in Forward & Backward. Because both vertical pixels in a thread share Δx, we can also fold the splat's opacity o directly into the base-2 exponential using the identity o · 2^p = 2^{p + log₂ o}. We precompute log₂ o once per visible Gaussian during preprocess; inside the tile rasterizers, each thread adds log₂ o into the horizontal A Δx² column constant once per column so that a single MUFU.EX2 instruction on each row yields the unclamped alpha 𝛼_raw = o · 2^p directly—eliminating both the per-pixel o · g multiplication and the per-row power <= 0 guard (replaced by a single per-column Δx² ≤ Δx_max² check in forward). In raster_bwd, folding o into 𝛼_raw has a second payoff: each row's log-falloff gradient g_k = 𝛼_raw T (c · ∇I − accum) already carries the factor of o required by the conic and mean derivatives, and only the single lane that outputs ∂L/∂o at the end of the warp reduction divides o back out once via 2^{−log₂ o}.

4. Zero-Multiply Row Derivative Combination & Deferred Δx Multiplication Across Vertical Shuffles. Even more FMA instructions disappear when a thread combines the conic and mean derivatives across its two vertical pixels, which sit at vertical offsets Δy₀ and Δy₀ − 1 from the splat center. Let g₀ and g₁ be the scalar log-alpha derivatives of the two rows. Because the row index shift k ∈ {0, 1} is zero on the top row and satisfies 1 = 1² = 1 on the bottom row, summing the Δy-weighted and Δy²-weighted derivatives across both rows requires zero extra multiplications by the row index: m₀ = g₀ + g₁,
S_y = g₀ Δy₀ + g₁(Δy₀ − 1) = m₀ Δy₀ − g₁,
S_yy = m₀ Δy₀² − g₁(2Δy₀ − 1). Better still, inside an 8 × 8 quadrant warp, threads that exchange values during the first two vertical shuffle steps (stride 16 and stride 8) lie in the exact same pixel column and therefore share the exact same horizontal offset Δx. Instead of multiplying m₀ and S_y by Δx and Δx² on every thread before the reduction, we sum the vertical terms down the column across the first two shuffle steps first and multiply by Δx only once per column before the horizontal shuffle steps.

5. Relieving the ALU Pipe in Forward via Sign-Bit Transmittance Termination. Finally, to relieve the integer/comparison ALU pipe in raster_fwd, we eliminate the separate boolean flag that tracked whether each pixel row had saturated. Instead, we encode saturation directly in the sign bit of the floating-point transmittance register T: when a pixel's updated transmittance crosses below the 10⁻⁴ threshold, we negate T so that a negative value marks the pixel as finished while preserving its magnitude for the final output image. A single floating-point sign check then tests whether the pixel is still active, and each thread accumulates its blended-splat bits in a register across the 32-splat chunk before combining them with a single hardware warp-wide OR reduction at the end of the chunk.

Together, these transformations cut raster_fwd on room 30k from 0.78 ms to 0.62 ms and raster_bwd from 1.35 ms to 0.93 ms (−31.1%), bringing the stage total from 3.15 ms down to 2.57 ms (9.8× below s00; 2.39 ms and 10.5× with bfloat16 SH; end-to-end, the room 30k step is 9.8× faster than Inria). Across all 9 Mip-NeRF 360 scenes, the mean training step reaches 3.29 ms at 30k (7.30× vs. Inria) and 2.89 ms at 7k (8.38× vs. Inria).