GEMM-GS: Accelerated 3D Gaussian Splatting

• GEMM-GS is a GPU rendering acceleration method that reformulates the Gaussian opacity calculation in 3DGS into a GEMM-compatible matrix multiplication while preserving exact blending semantics.
• It employs a hardware-aware algebraic transformation with intra-tile coordinate recoding to map complex quadratic evaluations to efficient Tensor Core operations.
• Empirical evaluations show up to 1.74× speedup and demonstrate composability with existing 3DGS accelerators across diverse 3D scenes.

GEMM-GS is a GPU rendering acceleration method for 3D Gaussian Splatting (3DGS) that targets the dominant bottleneck of the 3DGS pipeline—the blending stage—by reformulating Gaussian–pixel opacity evaluation into a GEMM-compatible matrix multiplication that can execute on Tensor Cores. Its defining claim is not an approximation of the rendering equation, but an equivalent algebraic transformation of the Gaussian power term, so that the original blending semantics are preserved while the implementation becomes compatible with high-throughput GEMM hardware paths on modern NVIDIA GPUs (Li et al., 2 Apr 2026).

1. Pipeline context and performance bottleneck

3DGS was introduced as an explicit-scene alternative to NeRF, whose point-sampling design can require over a million model inferences per image. In the formulation discussed by GEMM-GS, the standard 3DGS renderer proceeds through four stages: preprocessing, duplication, sorting, and blending. Preprocessing projects 3D Gaussians to the image plane, computes projected ellipse parameters, depth $d$, and color $\mathbf{c}$, and determines tile overlap. Duplication replicates each Gaussian for every tile it intersects. Sorting orders duplicated Gaussians within each tile by depth. Blending then composites the depth-sorted Gaussians for each tile in front-to-back order (Li et al., 2 Apr 2026).

The paper reports that, on an NVIDIA A100 and across scenes from Tanks & Temples, Deep Blending, and Mip-NeRF 360, the blending stage contributes nearly 70% of total rendering time. This localization of the bottleneck is central: existing acceleration methods often reduce model size or preprocessing cost, but do not change the fact that standard blending is dominated by scalar and elementwise operations rather than GEMM-like dense linear algebra. Since Tensor Core throughput has grown dramatically from Volta onward and can exceed CUDA-core throughput by over an order of magnitude, the absence of GEMM-shaped work in 3DGS leaves a large fraction of available hardware underutilized (Li et al., 2 Apr 2026).

The obstacle is the Gaussian opacity calculation itself. In vanilla 3DGS, opacity depends on a quadratic form in pixel offsets, involving $\Delta x_{ij}^2$, $\Delta y_{ij}^2$, and $\Delta x_{ij}\Delta y_{ij}$. That structure is not directly representable as a standard matrix multiply. GEMM-GS is therefore best understood as a hardware-aware algebraic recoding of the blending kernel rather than a change to scene representation, training objective, or compositing order.

2. Original blending equations and exact GEMM-compatible reformulation

For a tile with Gaussians sorted front-to-back, the rendered color of pixel $p_j$ is written as

$$\mathbf{C}_j = \sum_{i=1}^{N} T_i \alpha_{ij} \mathbf{c}_i, \quad T_i = \prod_{k=1}^{i-1} (1 - \alpha_{kj}),$$

where $\alpha_{ij}$ is the opacity of Gaussian $i$ at pixel $j$, $\mathbf{c}$0 is the Gaussian color, and $\mathbf{c}$1 is the accumulated transmittance. GEMM-GS leaves this front-to-back compositing rule unchanged. The transformation targets only the opacity exponent (Li et al., 2 Apr 2026).

The original opacity is

$\mathbf{c}$2

with

$\mathbf{c}$3

The exponent is expanded as

$\mathbf{c}$4

GEMM-GS introduces intra-tile coordinates by selecting a tile reference pixel $\mathbf{c}$5 with coordinates $\mathbf{c}$6 and writing

$\mathbf{c}$7

so that

$\mathbf{c}$8

Defining

$\mathbf{c}$9

the paper expands and regroupes the quadratic form into a 6D dot product: $\Delta x_{ij}^2$0 Thus

$\Delta x_{ij}^2$1

This identity is the mathematical core of GEMM-GS. Pixel-side vectors $\Delta x_{ij}^2$2 depend only on intra-tile coordinates, so for a fixed $\Delta x_{ij}^2$3 tile they are reusable across all tiles and can be precomputed offline once. Gaussian-side vectors $\Delta x_{ij}^2$4 depend on projected Gaussian parameters and on Gaussian position relative to the tile reference, so they are built online per tile. Stacking the Gaussian vectors into $\Delta x_{ij}^2$5 and pixel vectors into $\Delta x_{ij}^2$6 yields the intended matrix product

$\Delta x_{ij}^2$7

whose entries are the original $\Delta x_{ij}^2$8, followed by

$\Delta x_{ij}^2$9

The paper explicitly treats this as an equivalent transformation, not a lossy approximation. It also notes that some printed matrix dimensions are inconsistent, but the intended algebraic construction is clear (Li et al., 2 Apr 2026).

3. Tensor Core kernel architecture and runtime system

GEMM-GS preserves the tile structure of vanilla 3DGS: each $\Delta y_{ij}^2$0 image tile is rendered by one thread block of 256 threads, with one thread per pixel, and Gaussians overlapping the tile are processed in batches of 256. The CUDA kernel is organized as a three-stage pipeline. First, Gaussian indices for the current tile batch are loaded into shared memory. Second, each thread fetches one Gaussian’s attributes—projected coordinates $\Delta y_{ij}^2$1, covariance parameters $\Delta y_{ij}^2$2, opacity $\Delta y_{ij}^2$3, and color $\Delta y_{ij}^2$4—and cooperatively constructs one row of $\Delta y_{ij}^2$5. Third, the block performs the GEMM, reconstructs $\Delta y_{ij}^2$6, and executes standard front-to-back volume rendering for the tile (Li et al., 2 Apr 2026).

The implementation is explicitly double-buffered across Gaussian indices, Gaussian features, and $\Delta y_{ij}^2$7, so that loading the next batch overlaps with current computation. The paper states that asynchronous memory copy, introduced in Ampere, is used to overlap data movement and compute. This is important because GEMM-GS does not accelerate the entire blending stage uniformly; it accelerates the power-term evaluation while attempting to hide the memory-system cost of preparing Tensor Core operands.

A further implementation constraint is shared-memory capacity. The conceptual $\Delta y_{ij}^2$8 multiply is therefore decomposed into mini-batches. The paper states that each iteration multiplies a $\Delta y_{ij}^2$9 submatrix of $\Delta x_{ij}\Delta y_{ij}$0 with $\Delta x_{ij}\Delta y_{ij}$1 using warp-level PTX mma(m16n8k8) instructions. Each thread block contains 8 warps, and the implementation description states that $\Delta x_{ij}\Delta y_{ij}$2 warp-level MMA calls complete an effective $\Delta x_{ij}\Delta y_{ij}$3 multiplication. The notation is not fully consistent with the high-level matrix dimensions printed earlier in the paper, but the hardware strategy is clear: GEMM-GS lowers the transformed blending kernel to warp-tiled Tensor Core fragments and keeps synchronization largely at warp scope rather than full-thread-block scope (Li et al., 2 Apr 2026).

Data organization is correspondingly split. The pixel matrix $\Delta x_{ij}\Delta y_{ij}$4 is offline-precomputed because every tile shares the same intra-tile coordinate stencil. The Gaussian matrix $\Delta x_{ij}\Delta y_{ij}$5 is built online from per-tile Gaussian state. Shared memory stages sorted Gaussian lists and Gaussian attributes, while Tensor Cores consume the packed fragment layout generated from $\Delta x_{ij}\Delta y_{ij}$6 and $\Delta x_{ij}\Delta y_{ij}$7.

4. Empirical evaluation and measured speedups

The evaluation uses 13 scenes: Train and Truck from Tanks & Temples, Playroom and Drjohnson from Deep Blending, and bicycle, bonsai, counter, flowers, garden, kitchen, room, stump, and treehill from Mip-NeRF 360. Models are trained for 30K iterations using official 3DGS. Rendering latency is measured on NVIDIA A100 and H100 GPUs by rendering the full test set 10 times and averaging. Reported scene resolutions range roughly from $\Delta x_{ij}\Delta y_{ij}$8 to $\Delta x_{ij}\Delta y_{ij}$9 for Tanks & Temples and Deep Blending, and around $p_j$0 for Mip-NeRF 360, with Gaussian counts from roughly 1.04M to 4.74M (Li et al., 2 Apr 2026).

The headline result is a mean speedup of $p_j$1 over vanilla 3DGS on A100 and $p_j$2 on H100. On A100, per-scene latency reductions include truck from 4.28 ms to 2.78 ms ($p_j$3), train from 4.51 ms to 3.01 ms ($p_j$4), drjohnson from 6.16 ms to 4.16 ms ($p_j$5), playroom from 4.37 ms to 2.86 ms ($p_j$6), bicycle from 9.64 ms to 6.87 ms ($p_j$7), and treehill from 5.02 ms to 3.63 ms ($p_j$8) (Li et al., 2 Apr 2026).

Setting	Result	Platform
GEMM-GS vs vanilla 3DGS	$p_j$9 average speedup	A100
GEMM-GS vs vanilla 3DGS	$$\mathbf{C}_j = \sum_{i=1}^{N} T_i \alpha_{ij} \mathbf{c}_i, \quad T_i = \prod_{k=1}^{i-1} (1 - \alpha_{kj}),$$0 average speedup	H100
GEMM-GS at $$\mathbf{C}_j = \sum_{i=1}^{N} T_i \alpha_{ij} \mathbf{c}_i, \quad T_i = \prod_{k=1}^{i-1} (1 - \alpha_{kj}),$$1 resolution	$$\mathbf{C}_j = \sum_{i=1}^{N} T_i \alpha_{ij} \mathbf{c}_i, \quad T_i = \prod_{k=1}^{i-1} (1 - \alpha_{kj}),$$2 speedup	relative to vanilla
GEMM-GS at $$\mathbf{C}_j = \sum_{i=1}^{N} T_i \alpha_{ij} \mathbf{c}_i, \quad T_i = \prod_{k=1}^{i-1} (1 - \alpha_{kj}),$$3 resolution	$$\mathbf{C}_j = \sum_{i=1}^{N} T_i \alpha_{ij} \mathbf{c}_i, \quad T_i = \prod_{k=1}^{i-1} (1 - \alpha_{kj}),$$4 speedup	relative to vanilla

The resolution study is particularly revealing. The paper reports $$\mathbf{C}_j = \sum_{i=1}^{N} T_i \alpha_{ij} \mathbf{c}_i, \quad T_i = \prod_{k=1}^{i-1} (1 - \alpha_{kj}),$$5 speedup at $$\mathbf{C}_j = \sum_{i=1}^{N} T_i \alpha_{ij} \mathbf{c}_i, \quad T_i = \prod_{k=1}^{i-1} (1 - \alpha_{kj}),$$6 resolution and $$\mathbf{C}_j = \sum_{i=1}^{N} T_i \alpha_{ij} \mathbf{c}_i, \quad T_i = \prod_{k=1}^{i-1} (1 - \alpha_{kj}),$$7 at $$\mathbf{C}_j = \sum_{i=1}^{N} T_i \alpha_{ij} \mathbf{c}_i, \quad T_i = \prod_{k=1}^{i-1} (1 - \alpha_{kj}),$$8 resolution relative to vanilla 3DGS, indicating that the transformed Tensor Core path becomes more beneficial as the blending workload grows. The batch-size study also shows that smaller Gaussian batches degrade performance, because a $$\mathbf{C}_j = \sum_{i=1}^{N} T_i \alpha_{ij} \mathbf{c}_i, \quad T_i = \prod_{k=1}^{i-1} (1 - \alpha_{kj}),$$9 tile maps to 256 threads and a batch size $\alpha_{ij}$0 prevents even distribution of $\alpha_{ij}$1 construction work across the block, thereby lowering parallel efficiency and Tensor Core utilization (Li et al., 2 Apr 2026).

The paper does not provide PSNR, SSIM, or LPIPS tables in the provided text. Its quality argument is algebraic rather than metric-based: since the reformulation preserves $\alpha_{ij}$2, then $\alpha_{ij}$3 and the front-to-back compositing semantics are preserved as well.

5. Composition with existing 3DGS accelerators

A significant property of GEMM-GS is that it replaces only the blending-stage implementation and is therefore orthogonal to methods that attack other parts of the 3DGS stack. The paper explicitly evaluates composition with FlashGS, StopThePop, Speedy-Splat, c3dgs, and LightGaussian. Averaged over those baselines on A100, GEMM-GS provides an additional $\alpha_{ij}$4 speedup, with per-method averages of $\alpha_{ij}$5 on FlashGS, $\alpha_{ij}$6 on StopThePop, $\alpha_{ij}$7 on Speedy-Splat, $\alpha_{ij}$8 on c3dgs, and $\alpha_{ij}$9 on LightGaussian (Li et al., 2 Apr 2026).

The per-scene numbers illustrate the same pattern. For FlashGS + GEMM-GS, truck improves from 2.53 ms to 2.29 ms ($i$0), train from 2.60 ms to 2.15 ms ($i$1), and drjohnson from 2.08 ms to 1.58 ms ($i$2). For StopThePop + GEMM-GS, truck improves from 3.52 ms to 2.56 ms ($i$3) and playroom from 3.53 ms to 2.23 ms ($i$4). For Speedy-Splat + GEMM-GS, truck improves from 2.83 ms to 1.86 ms ($i$5), train from 2.93 ms to 1.88 ms ($i$6), drjohnson from 2.85 ms to 1.79 ms ($i$7), and playroom from 2.01 ms to 1.26 ms ($i$8) (Li et al., 2 Apr 2026).

This composability is central to the method’s systems significance. GEMM-GS does not reduce the number of Gaussians, alter tile overlap tests, or compress the scene representation. It accelerates the arithmetic core of blending. As a result, it composes naturally with preprocessing-oriented methods and with compression or pruning methods. The smaller gains on top of FlashGS than on top of vanilla 3DGS or c3dgs also align with the paper’s general limitation statement that benefits diminish when the workload is less blending-dominated.

6. Exactness claims, limitations, and broader significance

A persistent misconception would be to interpret GEMM-GS as turning all of blending into GEMM. The paper does not make that claim. What is transformed into GEMM is the Gaussian power computation. Exponentiation, $i$9-skipping when $j$0, early termination when transmittance is exhausted, and recursive front-to-back compositing remain non-GEMM control flow. Sorting is still required, because front-to-back compositing remains order dependent. The method is therefore best described as a partial algebraic lowering of blending to Tensor Core-compatible linear algebra rather than a wholesale replacement of the blending stage by a dense-matrix kernel (Li et al., 2 Apr 2026).

The implementation is also explicitly architecture-sensitive. It is designed around NVIDIA Tensor Core execution, uses PTX MMA instructions, and relies on asynchronous copy as introduced in Ampere. The kernel is substantially more complex than vanilla tile blending because it must manage packed fragments, shared-memory staging, mini-batching, and a three-stage double-buffered pipeline. The paper also notes that mini-batching is necessitated by shared-memory limits, that the precision discussion is limited in the provided text, and that some matrix-dimension notation in the printed equations is inconsistent (Li et al., 2 Apr 2026).

Within those limits, GEMM-GS is significant for two reasons. First, it demonstrates that a nominally non-GEMM neural rendering kernel can be recast exactly into a small polynomial basis whose evaluation is Tensor Core compatible. Second, it shows that this algebraic recoding is sufficiently practical to deliver consistent end-to-end rendering gains—$j$1 over vanilla 3DGS on A100, with larger gains at higher resolution and further gains when composed with other 3DGS accelerators—without changing the 3DGS scene model or the blending equation itself. The released code underscores that the contribution is intended as a deployable renderer-level systems optimization rather than only a symbolic reformulation (Li et al., 2 Apr 2026).