Differentiable rendering: the backward pass

A differentiable renderer computes not just an image but how the image would change if any scene parameter changed: ∂pixel/∂θ for every position, scale, opacity and colour. With that, "reconstruct a scene from photos" becomes "minimise the difference between rendered and real photos by gradient descent". That's how 3D Gaussian Splatting is trained, and how NeRF, inverse rendering and most of neural graphics work.

This lesson derives the backward pass of the splatting rasterizer from gfx-gaussian-splatting — the part of 3DGS that is hand-written CUDA rather than PyTorch autograd — shows how to verify a gradient, and how the optimiser and the densification heuristics use it. It is the core of "implementing research techniques": a paper gives you the forward model; you derive, implement and check the backward.

1. Reverse mode in one paragraph

A renderer is a long chain of small functions: parameters → activations (exp, sigmoid) → covariance → conic → α per pixel → composited colour → loss L, a single number. Reverse-mode differentiation (backpropagation) runs the chain backwards carrying the adjoint v̄ = ∂L/∂v of every intermediate v: if u = f(v, w), then v̄ += ū·∂f/∂v and w̄ += ū·∂f/∂w. One backward pass costs a small constant times the forward and yields ∂L/∂θ for all parameters — millions of them — which is why it beats finite differences (one forward per parameter) or forward mode.

The bookkeeping rules: every use of a variable adds to its adjoint (+=, which is why the GPU kernels use atomicAdd); the backward needs the forward's intermediate values, either stored or recomputed; and every branch the forward took (skip, clamp, early exit) must be replayed exactly.

One Gaussian at one pixel: forward and backward

Forward, left to right. The backward visits the same nodes right to left, multiplying the adjoint by each local derivative (shown under the edges' targets).

2. The derivative of compositing

Front-to-back compositing of n samples (nearest first) over a background b:

C = Σᵢ cᵢ αᵢ Tᵢ + T_N b, Tᵢ = Π_{j<i} (1 − αⱼ).

Colour. cᵢ appears once: ∂C/∂cᵢ = αᵢ Tᵢ.

Opacity. αᵢ appears in term i (as αᵢ) and in every term behind it (through T_j = Tᵢ(1 − αᵢ)·Π_{i<k<j}(1 − α_k) for j > i, and in T_N). Differentiating,

∂C/∂αᵢ = cᵢ Tᵢ − (1/(1 − αᵢ)) · [Σ_{j>i} c_j α_j T_j + T_N b].

Factor out Tᵢ using Tᵢ₊₁ = Tᵢ(1 − αᵢ), and name the colour accumulated behind i, seen from just behind it:

Sᵢ = Σ_{j>i} c_j α_j T_j / Tᵢ₊₁ (background excluded).

Then ∂C/∂αᵢ = Tᵢ (cᵢ − Sᵢ) − T_N b / (1 − αᵢ). Read it: making sample i more opaque shows more of its own colour and hides what's behind it, so the gradient is the difference — weighted by how much light reaches i at all. The background term is negative for b > 0: more opacity, less background.

Background. ∂C/∂b = T_N.

The suffix Sᵢ obeys a one-step recurrence from the back: Sᵢ = αᵢ₊₁ cᵢ₊₁ + (1 − αᵢ₊₁) Sᵢ₊₁, S_{n−1} = 0. And Tᵢ can be rebuilt from the back too: Tᵢ = Tᵢ₊₁ / (1 − αᵢ), starting from T_N. That is exactly the reference backward kernel (backward.cu, renderCUDA): it stores only the final T and the index of the last contributing splat per pixel, then walks each tile's list in reverse.

C++

// the reference, per pixel, walking splats back to front
T = T / (1.f - alpha);                                  // rebuild T_i
const float dchannel_dcolor = alpha * T;                // ∂C/∂c_i
for (int ch = 0; ch < C; ch++) {
    accum_rec[ch] = last_alpha * last_color[ch] + (1.f - last_alpha) * accum_rec[ch];   // S_i
    last_color[ch] = c;
    dL_dalpha += (c - accum_rec[ch]) * dL_dchannel;
    atomicAdd(&(dL_dcolors[global_id * C + ch]), dchannel_dcolor * dL_dchannel);
}
dL_dalpha *= T;
last_alpha = alpha;
dL_dalpha += (-T_final / (1.f - alpha)) * bg_dot_dpixel;   // the background term

Rebuilding T backwards

start / 5

The forward multiplies T by (1 − αᵢ) and keeps only T_N. The backward divides it back out — possible because α ≤ 0.99, so 1 − α ≥ 0.01.

C++

#include <cstddef>
#include <print>
#include <vector>

double composite(const std::vector<double>& a, const std::vector<double>& c, double bg) {
    double T = 1, C = 0;
    for (std::size_t i = 0; i < a.size(); ++i) { C += c[i] * a[i] * T; T *= 1 - a[i]; }
    return C + T * bg;
}

int main() {
    const std::vector<double> a{0.5, 0.25, 0.2, 0.6}, c{0.9, 0.1, 0.5, 0.3};
    const double bg = 1.0;
    double T_final = 1;
    for (double x : a) T_final *= 1 - x;

    double T = T_final, accum = 0, last_alpha = 0, last_color = 0;
    std::vector<double> d_alpha(a.size());
    for (std::size_t i = a.size(); i-- > 0;) {
        T /= 1 - a[i];
        accum = last_alpha * last_color + (1 - last_alpha) * accum;
        last_alpha = a[i];
        last_color = c[i];
        d_alpha[i] = T * (c[i] - accum) - T_final / (1 - a[i]) * bg;
    }
    for (std::size_t i = 0; i < a.size(); ++i) {
        auto ap = a, am = a;
        const double h = 1e-6;
        ap[i] += h;
        am[i] -= h;
        const double fd = (composite(ap, c, bg) - composite(am, c, bg)) / (2 * h);
        std::println("dC/dalpha[{}] = {:+.8f}   finite difference {:+.8f}", i, d_alpha[i], fd);
    }
}

3. Into the parameters: the chain rule

Each α is α = min(0.99, a), a = o·G, o = sigmoid(logit), G = exp(power), and for an isotropic 2D Gaussian power = −½|p − x|²/σ² with σ = e^{log σ}. Given ∂L/∂α:

The clamp. If a > 0.99 the forward used a constant, so ∂α/∂a = 0 and the pixel passes no gradient on (gsplat: if (opac * vis <= 0.999f) { … }). The original graphdeco kernel computes dL_dG = con_o.w * dL_dalpha regardless — a "straight-through" estimate that trains fine but disagrees with finite differences. Whichever you choose, the T rebuild must divide by the clamped α: that's what the forward multiplied in.
Opacity. ∂L/∂logit = ∂L/∂a · G · o(1 − o).
Shape and position. ∂L/∂power = ∂L/∂a · o · G; ∂power/∂x = −(x − p_x)/σ²; ∂power/∂log σ = |p − x|²/σ² = −2·power.

In full 3DGS the chain continues: ∂L/∂power → the 2D mean and the conic (∂power/∂conic = −½(dx², 2dx·dy, dy²)); the conic → the 2D covariance (∂Σ′⁻¹ = −Σ′⁻¹ ∂Σ′ Σ′⁻¹); Σ′ = J W Σ Wᵀ Jᵀ → Σ and the mean (through J); Σ = M Mᵀ → M = R S → scale and quaternion; and the SH colour → SH coefficients and the mean (the view direction depends on it). Each step is a small local Jacobian; the Mathematical Supplement for the gsplat Library writes them all down.

C++

#include <algorithm>
#include <cmath>
#include <print>

// α of one Gaussian at pixel (0, 0) as a function of its parameters.
double alpha(double x, double y, double log_sigma, double logit) {
    const double power = -0.5 * (x * x + y * y) * std::exp(-2 * log_sigma);
    const double a = std::exp(power) / (1 + std::exp(-logit));
    return std::min(0.99, a);
}

int main() {
    const double x = 0.8, y = -0.5, ls = std::log(1.3), logit = 1.2;
    const double inv_var = std::exp(-2 * ls);
    const double power = -0.5 * (x * x + y * y) * inv_var;
    const double G = std::exp(power), o = 1 / (1 + std::exp(-logit)), a = o * G;
    const double dL_da = a > 0.99 ? 0.0 : 1.0;    // L = α here, so dL/dα = 1
    const double dL_dpower = dL_da * o * G;
    const double analytic[4] = {dL_dpower * -x * inv_var, dL_dpower * -y * inv_var,
                                dL_dpower * -2 * power, dL_da * G * o * (1 - o)};
    const double h = 1e-6;
    const double numeric[4] = {
        (alpha(x + h, y, ls, logit) - alpha(x - h, y, ls, logit)) / (2 * h),
        (alpha(x, y + h, ls, logit) - alpha(x, y - h, ls, logit)) / (2 * h),
        (alpha(x, y, ls + h, logit) - alpha(x, y, ls - h, logit)) / (2 * h),
        (alpha(x, y, ls, logit + h) - alpha(x, y, ls, logit - h)) / (2 * h)};
    const char* names[4] = {"x", "y", "log sigma", "logit"};
    for (int k = 0; k < 4; ++k)
        std::println("{:>9}: {:+.8f}  vs  {:+.8f}", names[k], analytic[k], numeric[k]);
}

4. Gradient checking

Never trust a hand-derived backward until it matches central differences:

f′(x) ≈ (f(x + h) − f(x − h)) / 2h, error ≈ h²·|f‴|/6 (truncation) + ε·|f|/h (rounding).

Minimising over h gives h ≈ ∛ε·max(1, |x|): about 5·10⁻³ in float and 6·10⁻⁶ in double, with best-case errors around 10⁻⁵ and 10⁻¹¹. Too small an h is worse, not better — in float, 0.7f + 10⁻⁸ == 0.7f and the "derivative" is 0. Hence the rules: check in double (PyTorch's gradcheck insists on float64), divide by the step actually taken, (x + h) − (x − h), compare by relative error with a floor for gradients that should be zero, and test away from kinks (the 1/255 skip, the 0.99 clamp, early termination), where the function isn't differentiable.

Central-difference error vs step size

d/dx sin x, worst case over 400 points x ∈ [0.6, 0.8], measured. Left of each minimum, rounding error ∝ ε/h dominates; right of it, truncation ∝ h². Float bottoms out near 10⁻⁵ at h ≈ 5·10⁻³; below ~10⁻⁸ x ± h rounds to x and the quotient is 0/0.

C++

#include <cmath>
#include <limits>
#include <print>

template <class Real>
Real central(Real x, Real h) {
    const Real xp = x + h, xm = x - h;
    return (std::sin(xp) - std::sin(xm)) / (xp - xm);
}

int main() {
    const double exact = std::cos(0.7);
    for (double h : {1e-1, 1e-2, 5e-3, 1e-4, 1e-6, 1e-8, 1e-10}) {
        const float f = central(0.7f, static_cast<float>(h));
        const double d = central(0.7, h);
        std::println("h = {:7.0e}   float err {:9.2e}   double err {:9.2e}", h,
                     std::fabs(f - exact), std::fabs(d - exact));
    }
    std::println("cbrt(eps): float {:.2e}, double {:.2e}",
                 std::cbrt(std::numeric_limits<float>::epsilon()),
                 std::cbrt(std::numeric_limits<double>::epsilon()));
}

(At h = 10⁻⁸ in float, x + h rounds to x, so xp − xm is 0 and the quotient is NaN — a check that prints NaN is telling you the step is too small.)

5. Losses and optimisers

Loss. 3DGS minimises L = (1 − λ)·L1 + λ·D-SSIM with λ = 0.2. L1 = mean |I − I*| (gradient sign(I − I*)/N) is robust to outliers; D-SSIM = 1 − SSIM, where SSIM compares local means, variances and covariance over 11×11 Gaussian windows, rewards structure and contrast. Quality is reported as PSNR = 10·log₁₀(1/MSE) dB (for values in [0, 1]; +6 dB ≈ half the RMS error), plus SSIM and LPIPS. The build task uses MSE because its gradient, 2(I − I*)/N, is the simplest.

SGD steps θ −= lr·g: the step scales with the gradient, and positions, opacities and colours have gradients that differ by orders of magnitude. Adam (Kingma & Ba 2015) normalises per parameter:

m ← β₁m + (1 − β₁)g, v ← β₂v + (1 − β₂)g², m̂ = m/(1 − β₁ᵗ), v̂ = v/(1 − β₂ᵗ), θ ← θ − lr·m̂/(√v̂ + ε).

The bias correction matters: m and v start at 0, so without it the first step is lr·0.1/√0.001 ≈ 3.16·lr. With it, the step is ≈ lr regardless of g's scale.

3DGS uses one learning rate per parameter group: positions 1.6·10⁻⁴ × scene extent, decayed exponentially to 1.6·10⁻⁶ over 30k iterations; SH DC 2.5·10⁻³ and the higher bands 20× less; opacity 0.05 (0.025 in later code); log-scale 5·10⁻³; rotation 10⁻³; Adam ε = 10⁻¹⁵. In PyTorch that's one torch.optim.Adam with a list of {"params": …, "lr": …, "name": …} groups.

C++

#include <cmath>
#include <print>

// Minimise f(x, y) = 100 x² + 0.01 y² (gradients differ by 10⁴) from (1, 1).
int main() {
    double sx = 1, sy = 1;                       // SGD
    double ax = 1, ay = 1, m[2] = {0, 0}, v[2] = {0, 0};   // Adam
    const double b1 = 0.9, b2 = 0.999, eps = 1e-8;
    for (int t = 1; t <= 200; ++t) {
        sx -= 0.009 * 200 * sx;                  // the largest stable SGD step for x
        sy -= 0.009 * 0.02 * sy;
        const double g[2] = {200 * ax, 0.02 * ay};
        double* p[2] = {&ax, &ay};
        for (int i = 0; i < 2; ++i) {
            m[i] = b1 * m[i] + (1 - b1) * g[i];
            v[i] = b2 * v[i] + (1 - b2) * g[i] * g[i];
            const double mh = m[i] / (1 - std::pow(b1, t)), vh = v[i] / (1 - std::pow(b2, t));
            *p[i] -= 0.01 * mh / (std::sqrt(vh) + eps);
        }
    }
    std::println("SGD:  x = {:+.5f}  y = {:+.5f}", sx, sy);
    std::println("Adam: x = {:+.5f}  y = {:+.5f}", ax, ay);
}

SGD's learning rate is capped by the stiff x direction, so y has barely moved after 200 steps; Adam's normalised steps (≈ lr each) bring both coordinates down together.

6. Adaptive density control

Gradient descent can move, resize and recolour Gaussians but not create or delete them. 3DGS adds a heuristic loop every 100 iterations (from iteration 500 to 15k), driven by the average magnitude of each Gaussian's screen-space mean gradient — a by-product of the same backward pass:

Densify and prune

Large view-space gradients mean the region is poorly reconstructed. Small Gaussians there are cloned, large ones split; transparent or huge ones are removed. Every 3000 iterations all opacities are reset near 0 so floaters can die.

Later work replaces the heuristics with something principled: MCMC-style relocation (3DGS as Markov Chain Monte Carlo, Kheradmand et al. 2024), absolute gradients (AbsGS), or budgeted densification (Taming 3DGS) — gsplat ships several as "strategies".

7. Where custom backward passes live: PyTorch

PyTorch autograd differentiates anything built from its ops, but a per-pixel blending loop over millions of splats would be hopelessly slow as tensor ops. So the rasterizer is a custom autograd function: a CUDA forward and a hand-written CUDA backward, glued in with torch.autograd.Function.

In Python (PyTorch)

Text

import torch
import my_rasterizer_cuda as _C            # pybind11 extension: forward + backward kernels

class RasterizeGaussians(torch.autograd.Function):
    @staticmethod
    def forward(ctx, means2d, conics, opacities, colors, background, width, height):
        image, final_T, n_contrib = _C.forward(means2d, conics, opacities, colors,
                                               background, width, height)
        # save what the backward must replay; recompute the rest
        ctx.save_for_backward(means2d, conics, opacities, colors, background,
                              final_T, n_contrib)
        ctx.size = (width, height)
        return image

    @staticmethod
    def backward(ctx, dL_dimage):
        means2d, conics, opacities, colors, background, final_T, n_contrib = ctx.saved_tensors
        g_means2d, g_conics, g_opac, g_colors = _C.backward(
            means2d, conics, opacities, colors, background, final_T, n_contrib,
            dL_dimage.contiguous(), *ctx.size)
        # one gradient per forward input; None for non-tensor / non-differentiable ones
        return g_means2d, g_conics, g_opac, g_colors, None, None, None

image = RasterizeGaussians.apply(means2d, conics, opacities, colors, bg, W, H)
loss = 0.8 * (image - gt).abs().mean() + 0.2 * (1 - ssim(image, gt))
loss.backward()                                       # calls our backward, then autograd
torch.autograd.gradcheck(RasterizeGaussians.apply,   # in float64, on a tiny scene
                         (m.double().requires_grad_(), ...), eps=1e-6, atol=1e-5)

Everything before the rasterizer (activations, Σ = R S Sᵀ Rᵀ, SH evaluation, even the projection in some codebases) can stay in autograd; only the hot loop is custom.

In an interview

Expect "derive ∂C/∂αᵢ for front-to-back compositing" (answer: Tᵢ(cᵢ − Sᵢ) minus the background term, and why the backward walks back to front), "how do you know your gradient is right?" (central differences in double, h ≈ ∛ε, relative error, avoid kinks), "why Adam and not SGD here?", and "what happens at α = 0.99?".

Color theme