Efficient Differentiable Hardware Rasterization for 3D Gaussian Splatting

Kerbl et al. (2023) proposed 3D Gaussian Splatting (3DGS) as a differentiable scene representation that models 3D geometry using anisotropic 3D Gaussians. Each Gaussian primitive is parameterized by: (1) a center position $\mu_{\text{3D}}\in\mathbb{R}^{3}$, (2) opacity $o\in[0,1]$, (3) view-dependent color encoded through spherical harmonics (SHs), and (4) a covariance matrix $\mathbf{\Sigma_{\text{3D}}}$ constructed from the learnable rotation parameters $\mathbf{R}$ and scaling parameters $\mathbf{S}$ as $\mathbf{\Sigma_{\text{3D}}}=\mathbf{RS}\mathbf{S}^{\top}\mathbf{R}^{\top}$. The influence of a Gaussian (i.e. density) at point $\mathbf{x}$ follows:

Unlike implicit volumetric representations (e.g., NeRF Mildenhall et al. (2021)), 3DGS enables real-time rendering via explicit splat rasterization while preserving differentiability.

Following the formulation in Kerbl et al. (2023), 3DGS rendering first projects 3D Gaussians to screen space. Given a Gaussian splat with 3D mean $\mu_{\text{3D}}$ and covariance $\mathbf{\Sigma}_{\text{3D}}$, its 2D projection is computed using the view matrix $\mathbf{W}$ and the projection Jacobian $\mathbf{J}$:

yielding the 2D Gaussian density function $\alpha_{\text{2D}}(\mathbf{u}|o,\mu_{\text{2D}},\Sigma_{\text{2D}})$, where $\mathbf{u}$ denotes 2D screen-space coordinates and all parameters are projected to screen space.

Depth-sorted splats (indexed $1,\dots,N$) are alpha-composited to compute the final pixel color $C$. Let $c_{i}$ denote the view-dependent color of splat $i$, and $\alpha_{i}=\alpha_{\text{2D}}(\mathbf{u})$ (ignore $o,\mu_{\text{2D}},\Sigma_{\text{2D}}$ for simplicity) as the transparency at pixel $\mathbf{u}$. The blending equations are:

Where $T_{i}$ represents accumulated transmittance. To enable gradient-based optimization, the renderer computes derivatives of the loss $\mathcal{L}(C)$ with respect to splat parameters. The key gradients $\partial\mathcal{L}/\partial\alpha_{i}$ and $\partial\mathcal{L}/\partial c_{i}$ propagate to spatial parameters ($\mu_{\text{3D}}$, $\mathbf{\Sigma}_{\text{3D}}$) via chain rule differentiation. From Equation 2, we can derive:

The gradient computation supports both front-to-back($i=1,\dots,N$) and back-to-front ($i=N,\dots,1$) processing of depth-sorted splats.

Although Gaussian splats have infinite theoretical support (since $\alpha_{\text{2D}}(\mathbf{u})>0$ for all $\mathbf{u}$), practical implementations discard marginal contributions with $\alpha_{\text{2D}}(\mathbf{u})<1/255$, a threshold aligned with 8-bit color precision. This early culling applies to both color blending and gradient computation.

The differentiable tile-based rasterizer proposed by Kerbl et al. (2023) partitions the screen into $16\times 16$ tiles. For each splat, intersecting tiles are computed according to the splat boundary, and buffers are dynamically allocated to store splat-tile pairs. These pairs are sorted by tile ID and depth to generate per-tile depth-ordered splat lists. During the forward pass, pixel colors are computed in front-to-back order as Equation 2.

In the backward pass, splats are processed in reversed depth order. Initial gradient computation focuses on per-pixel $\partial\mathcal{L}/\partial c_{i}$ and $\partial\mathcal{L}/\partial\alpha_{i}$, followed by propagation to intermediate derivatives, such as $\partial\mathcal{L}/\partial\mathbf{\Sigma_{\text{2D},i}}^{-1}$, $\partial\mathcal{L}/\partial o_{i}$ and $\partial\mathcal{L}/\partial\mu_{\text{2D},i}$. These gradients are then aggregated per-splat via atomic operations. High-dimensional parameter gradients (3D positions $\mu_{\text{3D}}$, SH coefficients, etc.) are processed in a dedicated subsequent pass to maintain performance.

Kerbl et al. (2023) also implements an optimization termed T-Culling, which early terminates front-to-back splat traversal in the forward pass when the accumulated transparency product $T_{i}$ drops below $0.0001$, with its results stored and reused to accelerate the backward pass as well. Leveraging the monotonic decay of $T_{i}$ with depth order, this heuristic skips deeper splats whose contributions to pixel color and gradients become negligible. This optimization has been widely adopted in subsequent tile-based rasterizers as in gsplatYe et al. (2024). Though being straightforward for tile-based methods, T-Culling is incompatible with hardware 3DGS rasterization when employing graphics pipeline fixed-function blending, as GPU’s fixed-function blending intrinsically lacks support for such optimization.

Programmable Blending is a hardware rasterization technique that overcomes the predefined constraints of traditional fixed-function blending, enabling custom color mixing logic to be implemented in shader code. This capability, supported through graphics API extensions on modern GPUs, provides two essential features: (1) race-condition-free operations on per-pixel resources in fragment shaders, and (2) execution order guarantees which match the rasterization sequence for resource access.

Tile-based architecture GPUs prevalent in mobile platforms natively support programmable blendingArm Holdings plc (2025); Apple Inc. (2024). These GPUs guarantee that fragment shaders within a tile execute in strict rasterization order, requiring only framebuffer synchronization to make programmable blending available. For example, within the Vulkan graphics API Bailey (2019), this capability is formally specified through the \seqsplitVK_EXT_shader_tile_image and \seqsplitVK_EXT_rasterization_order_attachment_access extensions, which permit safe read-modify-write operations on the same attachment in fragment shaders while eliminating explicit memory barrier requirements.

Desktop GPUs programs can implement programmable blending via the \seqsplitVK_EXT_fragment_shader_interlock Vulkan extension, which defines a critical section through \seqsplitbeginInvocationInterlock() and \seqsplitendInvocationInterlock() calls, enabling safe access to pixel-local data structures in fragment shaders. Execution order can be further configured within these critical sections as rasterization-ordered. However, unlike tile-based mobile GPUs that guarantee rasterization-order execution, most desktop GPUs employ out-of-order fragment shader scheduling for improved hardware utilization. This architectural approach requires fragment shader interlocks to enforce synchronization via pipeline stalls or execution order constraints, incurring measurable performance overhead.

Programmable blending enables persistent per-pixel state management, making per-pixel gradient computation within fragment shaders feasible. While both front-to-back ($1\to N$) and back-to-front ($N\to 1$) traversal orders are valid for gradient computation, we empirically find that front-to-back traversal better aligns with T-Culling optimizations. This stems from hardware blending’s limitations in forward passes (not able to apply T-Culling), thereby constraining the optimization exclusively to backward passes. A back-to-front traversal in backward pass would require dividing out culled $1-\alpha_{i}$ terms for $T_{i}$ values needed by T-Culling, leading to extra memory accesses. Front-to-back traversal avoids this by directly computing $T_{i}$ through incremental updates.

The front-to-back gradient computation proceeds as follows. Firstly, both Equation 3 and Equation 4 depend on the following recursively updated transmittance:

Additionally, the summation term in Equation 4 requires recursive computation. We define an auxiliary sequence $C^{\prime}_{i}$:

Substituting Equation 6 into Equation 4 yields the reformed gradient expression:

Using these relations, we compute $\partial\mathcal{L}/\partial c_{i}$ in Equation 3 and $\partial\mathcal{L}/\partial\alpha_{i}$ in Equation 7 through front-to-back recurrence of $T_{i}$ and $C^{\prime}_{i}$, serving as the foundation for deriving any other gradients.

Our backward pass draws splats in front-to-back order, and Procedure 1 illustrates our gradient computation approach within the fragment shader, enabled by programmable blending through \seqsplitVK_EXT_fragment_shader_interlock (targeting desktop GPUs). During the rasterization-ordered critical section enforced by the interlocks, we maintain and retrieve values of $C^{\prime}_{i}$ and $T_{i}$ by reading/writing a dedicated texture $\text{tex}_{C^{\prime},T}$. The shader then loads $\partial\mathcal{L}/\partial C$ from a read-only texture $\text{tex}_{\partial\mathcal{L}/\partial C}$, then computes $\partial\mathcal{L}/\partial c_{i}$ and $\partial\mathcal{L}/\partial\alpha_{i}$ using analytical derivatives derived from Equation 3 and 7. Similar to Kerbl et al. (2023), these gradients are subsequently propagated to intermediate splat parameters $\theta_{i}$ (e.g., 2D positions, colors) via the chain rule, denoted as $\partial\mathcal{L}/\partial\theta_{i}$.

Additionally, we maintain a boolean \seqsplitcullingFlag to track fragment rejection status, which combines two culling criteria: splat boundary culling ($\alpha_{i}<1/255$) and T-Culling ($T_{i}<0.0001$). Fragments marked by \seqsplitcullingFlag skip all gradient computations, reducing both texture accesses and arithmetic operations.

Prior to detailing our method, we briefly clarify the architectural foundations of subgroups and quads in modern GPUs. Modern GPU architectures employ two key execution models for fragment shaders:

Subgroups (warps/wavefronts): Minimal schedulable thread blocks (32/64 threads on NVIDIA/AMD GPUs respectively) supporting intra-group data exchange via operations like \seqsplitsubgroupShuffle.

Quads: $2\times 2$ fragment blocks enabling derivative calculations, with helper invocations (detectable via \seqsplitgl_HelperInvocation) handling primitive boundary conditions.

Fragment shaders rely exclusively on these units for thread communication due to lack of shared memory, constraining gradient aggregation to intra-quad and intra-subgroup operations.

To mitigate performance bottlenecks from frequent \seqsplitatomicAdd operations, we propose aggregating gradients computed across multiple fragment shader threads using quad and subgroup intrinsics, followed by sparse \seqsplitatomicAdd calls from selected threads to update global memory buffers.

A straightforward approach involves quad reduction: Since hardware rasterization processes splat primitives at quad granularity, all threads within a quad inherently belong to the same splat. By leveraging quad intrinsics (e.g., Vulkan’s \seqsplitsubgroupQuadSwapVertical) to accumulate gradients within the quad, we perform partial sums and delegate a single thread per quad to execute \seqsplitatomicAdd. This strategy reduces \seqsplitatomicAdd calls to below $30\%$ of the baseline approach and achieves 3.05$\times$ speedup (see Table I).

To further reduce \seqsplitatomicAdd invocations, we explore subgroup gradient reduction. While there exists no guarantee that fragments within a subgroup originate from the same splat, we propose the Splat-Subgroup Cohesion Hypothesis: GPUs exhibit a statistical tendency to schedule fragments from the same splat (rendered as two adjacent triangle primitives) within a single subgroup, leveraging spatial locality for memory access efficiency during rasterization. Although not universally guaranteed, this hypothesis facilitates effective subgroup-level gradient reduction in most practical scenarios.

We validate this hypothesis through preliminary experiments on an RTX 4080 GPU. Qualitative analysis of splat-subgroup coherency (Figure 1) shows that splat-incoherent fragments (processed in subgroups containing fragments from multiple splats) mainly occur with small splats on screen, while larger splats predominantly produce splat-coherent fragments. This indicates high splat-coherency under typical rendering conditions. Quantitative measurements (Figure 2) further confirm these findings.

These results confirm that our Splat-Subgroup Cohesion Hypothesis holds robustly on RTX 4080 GPUs. While architecture-specific, the observed pixel locality patterns suggest potential applicability to other modern GPUs, requiring further cross-platform validation.

Building on this hypothesis, our subgroup reduction pipeline first assesses splat cohesion using Vulkan’s \seqsplitsubgroupAllEqual intrinsic to verify uniform splat IDs across active threads. For cohesive subgroups, gradients are aggregated via \seqsplitsubgroupAdd, with a single \seqsplitatomicAdd operation performed by a selected active lane.

Standalone subgroup reduction reduces \seqsplitatomicAdd invocations to 6% of baseline levels, and combining it with quad reduction further decreases the ratio to approximately 5%.

While minimizing \seqsplitatomicAdd invocations improves performance, excessive aggregation via subgroup reduction risks overloading streaming multiprocessors (SMs) with computation, creating an imbalance between SM workloads and L2 atomic unit (ROP) utilization that incurs performance penaltiesDurvasula et al. (2023, 2025). Inspired by Durvasula et al. (2023, 2025), we introduce a balancing threshold $X$ and restrict subgroup reduction to cases where the number of active threads requiring gradient aggregation exceeds the threshold. The optimal balancing threshold $X$ is hardware-and-scene-dependent and need to be selected through performance profiling.

Our hybrid gradient reduction method, combining quad-based and subgroup-based reduction strategies, is detailed in Procedure 2. This approach employs Vulkan GLSL intrinsics for subgroup and quad operations. The algorithm initiates by filtering out noneffective threads through a \seqsplitreduceFlag, which combines \seqsplitcullingFlag (from Procedure 1) and hardware-generated \seqsplitgl_HelperInvocation. The \seqsplitreduceFlag excludes helper invocations (which cannot perform memory writes or framebuffer updates) and fragments previously discarded by splat boundaries or T-Culling. A bitmask \seqsplitreduceMask, constructed via \seqsplitsubgroupBallot, identifies active threads eligible for gradient aggregation.

Then, the method dynamically selects between subgroup and quad operations based on runtime conditions. The algorithm first evaluates subgroup reduction feasibility by checking splat ID coherence across active threads and verifying that the valid fragment count, measured via \seqsplitsubgroupBallotBitCount, exceeds the balancing threshold $X$. When these criteria are met, gradients are aggregated within the subgroup. If unsatisfied, the method falls back to quad reduction, aggregating gradients across $2\times 2$ pixel blocks, and then restricting atomic updates to quad-local threads by zeroing out \seqsplitreduceMask bits for threads outside the current quad. The quad reduction fallback acts as a safety net for cases where subgroup coherence cannot be guaranteed, preserving performance across diverse scene complexities.

Finally, the least-significant thread in the current quad or subgroup reduction scope (identified via \seqsplitsubgroupBallotFindLSB) performs a single \seqsplitatomicAdd, accumulating the reduced gradients to global memory.

This dual-strategy approach is derived through rigorous profiling across varying 3DGS workloads, with empirical validation of its efficacy provided in Table I.