Foundations of Physically Based Differentiable Rendering

Introduction

×

Differentiable rendering asks a simple question: if a renderer maps scene parameters to an image, can we differentiate that map? If yes, then geometry, materials, lights, and cameras can be optimized from image-space losses such as reconstruction error, perceptual losses, or task-specific objectives.

The difficulty is that rendering is not just a smooth program. It is an integral over paths, visibility changes discontinuously, and Monte Carlo estimators have sampling choices that may themselves depend on the parameters. These notes has the following flow: automatic differentiation, why naive AD fails for visibility, then boundary-aware Monte Carlo estimators, and finally the physics-based formulations used in modern differentiable renderers (Mitsuba).

The notes assume basic familiarity with physically based rendering and the rendering equation. The goal here is to make the differentiable part clear enough that the papers become much easier to read. I focused on surface transport; participating media and null-collision estimators are treated as a separate advanced topic (will probably make separate notes in future).

Note: Many of the diagrams and visualizations in this post are adapted from the respective original research papers and Delio Vicini's PhD thesis [2].

I will be using the following custom differentiable renderer to do simplified implementations throughout this post. (The full code is available on GitHub):

Path

import torch
from scene import Scene
from camera import Camera
from ray import Ray


class PathTracer:
    """Monte Carlo path tracer."""

    def __init__(self, max_depth=5, num_samples=128):
        self.max_depth = max_depth
        self.num_samples = num_samples

    def sample(self, scene: Scene, camera: Camera):
        """Render via path tracing."""
        accum_L = torch.zeros_like(camera.origins)

        for _ in range(self.num_samples):
            ray = camera.sample()
            β = torch.ones_like(ray.origins)  # Path throughput
            L = torch.zeros_like(ray.origins)  # Accumulated radiance

            for depth in range(self.max_depth):
                si = scene.intersect(ray)

                # Direct emission from light sources
                Le = β * si.emission
                L = L + torch.where(si.is_valid(), Le, 0.0)

                # Flip normal for rays hitting the back face
                shading_n = torch.where((si.n * ray.dirs).sum(-1, keepdim=True) > 0.0, -si.n, si.n)

                # Sample BSDF to get new ray direction and update throughput
                bsdf_wi, bsdf_value, bsdf_pdf = si.bsdf.sample(-ray.dirs, shading_n)

                # Update ray for next bounce
                ray = Ray(si.p + shading_n * 1e-3, bsdf_wi)
                β = torch.where(si.is_valid(), β * bsdf_value / bsdf_pdf, 0.0)

            accum_L += L

        return accum_L / self.num_samples

Differentiation Methods

Finite Differences and Simultaneous Perturbation

The simplest approach to gradient computation is the finite difference (FD) method. For a scalar function $f: \mathbb{R} \to \mathbb{R}$, the forward difference estimator approximates the derivative at $x$ using a small step size $h > 0$:

A commonly used variant is the central difference method, which provides a significantly better approximation with quadratic truncation error $\mathcal{O}(h^2)$ instead of linear $\mathcal{O}(h)$:

Finite differences are inherently biased because evaluating $f$ at a non-zero step $h$ returns a local spatial average of the true derivative rather than its point value at $x$:

where $K_h(u) = \frac{1}{2h}\mathbf{1}_{[-h,h]}(u)$ is a rectangular boxcar kernel of width $2h$. Thus, a finite difference returns a local average of $f'$ smoothed over $[-h, h]$.

This blurring can produce inaccurate gradients when $f$ contains high-frequency features or discontinuities. The bias vanishes as $h \to 0$, but infinitely small steps are numerically fragile: in floating-point arithmetic, catastrophic cancellation degrades precision, and in Monte Carlo rendering, small differences become overwhelmed by stochastic sampling noise.

It is straightforward to apply finite differences to a renderer by generating the image once with the original parameter and once with the perturbed parameter. With a Monte Carlo renderer, though, the evaluation of $f$ is noisy, and if $f(x + h)$ and $f(x)$ are evaluated independently, the FD estimator needs an enormous number of samples to converge. Using Common Random Numbers (CRN) (seeding both evaluations with identical random number generator streams) resolves this because the Monte Carlo noise in $f(x+h)$ and $f(x)$ becomes strongly correlated, a significant portion of the variance cancels out.

Fundamentally, finite differences do not scale to functions with many input parameters due to the curse of dimensionality. For inverse rendering with a scene parameter vector $\mathbf{x} = (x_1, \dots, x_n)^T \in \mathbb{R}^n$ (such as meshes, textures, and volumes with millions of degrees of freedom), central differences would require rendering the image $2n$ times per gradient step ($f(x_1, \dots, x_i \pm h, \dots, x_n)$ for every parameter $i$). This is computationally impractical. An alternative is Simultaneous Perturbation Stochastic Approximation (SPSA), which estimates high-dimensional gradients by randomly offsetting all parameters at once.

For $f: \mathbb{R}^n \to \mathbb{R}$, SPSA estimates the gradient vector using only two function evaluations, regardless of the dimension $n$:

where $\boldsymbol{\Delta}^{-1} = (\Delta_1^{-1}, \dots, \Delta_n^{-1})^T$ denotes component-wise inversion, so each gradient component is estimated as:

The random perturbation vector $\boldsymbol{\Delta}$ has entries drawn independently from a mean-zero, symmetric distribution with bounded inverse moments, in practice, almost always a Rademacher distribution (each entry $\Delta_i = \pm 1$ with equal probability). A Gaussian $\boldsymbol{\Delta}$ cannot be used here: its probability density is non-zero at $0$, so $\Delta_i^{-1}$ has infinite variance and the estimator blows up.

While SPSA requires only two function evaluations per step regardless of input dimensionality, it introduces additional stochastic direction variance into the gradient estimates. This requires careful tuning of the step size $h$ to achieve good convergence. Consequently, derivative-free methods cannot compete with gradient descent using true infinitesimal gradients computed via automatic differentiation.

Code: Generating cbox translation gradient via Finite Difference

scene_path = 'scenes/cbox/cbox_bunny.xml'
if not os.path.exists(scene_path): scene_path = 'scenes/cbox/scene.xml'
scene, cam, _ = load_scene_from_xml(scene_path, device=device, override_res=512)
img_base = render_crn(scene, cam, integrator_fd).cpu().numpy()
mesh = scene.get_mesh("bunny")

h = 0.01
mesh.translate([h, 0.0, 0.0])
img_pert = render_crn(scene, cam, integrator_fd).cpu().numpy()
mesh.translate([-h, 0.0, 0.0])

fd_grad = np.mean((img_pert - img_base) / h, axis=-1)
plot_fd(img_base, img_pert, fd_grad, "X-shift", h, vmin=-25, vmax=25)

Code: Generating cbox albedo gradient via Finite Difference

scene_path = 'scenes/cbox/cbox_bunny.xml'
if not os.path.exists(scene_path): scene_path = 'scenes/cbox/scene.xml'
scene, cam, _ = load_scene_from_xml(scene_path, device=device, override_res=512)
img_base = render_crn(scene, cam, integrator_fd).cpu().numpy()
mesh = scene.get_mesh("bunny")

h = 0.1
mesh.albedo += torch.tensor([h, 0.0, 0.0], device=device)
img_pert = render_crn(scene, cam, integrator_fd).cpu().numpy()
mesh.albedo -= torch.tensor([h, 0.0, 0.0], device=device)

fd_grad = np.mean((img_pert - img_base) / h, axis=-1)
plot_fd(img_base, img_pert, fd_grad, "Albedo", h, vmin=-0.5, vmax=0.5)

Code: Generating roughness gradient

scene_path = 'scenes/teapot/teapot.xml'
if not os.path.exists(scene_path): scene_path = 'scenes/teapot/scene.xml'
scene, cam, _ = load_scene_from_xml(scene_path, device=device, override_res=512)
img_base = render_crn(scene, cam, integrator_fd).cpu().numpy()
mesh = scene.get_mesh("teapot")

h = 0.1
mesh.roughness += torch.tensor([h], device=device)
img_pert = render_crn(scene, cam, integrator_fd).cpu().numpy()
mesh.roughness -= torch.tensor([h], device=device)

fd_grad = np.mean((img_pert - img_base) / h, axis=-1)
plot_fd(img_base, img_pert, fd_grad, "Roughness", h, vmin=-25, vmax=25)

Automatic Differentiation

Instead of approximating derivatives with finite differences or deriving a large formula by hand, we can use automatic differentiation (AD). AD runs the original computation as a sequence of simple operations and applies the chain rule to each operation. For two scalar functions, the chain rule is:

For vector-valued functions, the same rule becomes a product of Jacobian matrices.

AD evaluates the required Jacobian products without constructing the full matrices.

This evaluates the chain rule without finite-difference truncation error, though the computation remains subject to ordinary floating-point roundoff. Automatic differentiation was introduced between the 1950s and 1970s and later became widely used for neural network training. The following overview focuses on the parts that matter for inverse rendering. A good resource on the implementation of this concept by Andrej Karpathy can be found on YouTube here.

Computation graphs. The central idea is to view a computation as a graph of operations. The individual operations are nodes, and the derivatives of individual steps are assigned to the graph’s edges. For example, consider:

In a computer program, the evaluation of this expression could be implemented as a sequence of steps:

a = 2 * x
b = a * y
c = sin(b)
d = x * x
e = d * c

The corresponding computation graph is shown in Fig. . Each edge stores a local derivative. For example, because $b=ay$, the edge from $a$ to $b$ has derivative $\partial b/\partial a=y$. AD combines these local derivatives to obtain the derivative of the final output. The ordinary evaluation of the function is called the primal computation, and the saved operations are often called a tape or Wengert tape.

Forward-mode Differentiation

A key choice in AD algorithms is the directionality of the gradient computation. Forward mode starts at one input and carries its derivative toward the outputs. If the computation has a single differentiable input variable, but many outputs, it is efficient to evaluate gradients from the variable to the output in the forward direction.

Mathematically, the differentiation turns into a series of Jacobian-vector products (JVP). For a function $\mathbf{y}=f(\mathbf{x})$, forward-mode AD computes the output gradient $\delta_{\mathbf{y}}$ as the product of the Jacobian $\mathbf{J}_f$ with the input gradient $\delta_{\mathbf{x}}$:

Here and in the following, we use $\delta$ to denote vectors and scalars that are inputs and outputs of Jacobian products. To differentiate with respect to $x$, we first initialize a variable $\delta x = 1$ (and $\delta y = 0$) and then traverse the graph from left to right, in each step multiplying the derivative value by the stored edge weights. Every later $\delta v$ is computed alongside its ordinary primal value $v$.

The interactive simulation below first computes the derivative of the output $e$ with respect to $x$, then repeats the sweep for $y$:

❮ ❯

Step 1 / 14

Each forward sweep gives the derivative with respect to one chosen input. In the end, the variable $\delta e$ contains the full derivative. Forward mode never explicitly computes and stores the full Jacobian $\mathbf{J}_f$ of the program, but its cost grows with the number of inputs because the graph must be traversed again for each one.

Forward-mode differentiation can be formalized by using dual numbers. Similar to a complex number, a dual number $a + \epsilon b$ stores a real part $a$ and a dual part $b$. The symbol $\epsilon$ satisfies $\epsilon^2 = 0$ and hence the product of two dual numbers is:

The first part, $ac$, is the ordinary product. The coefficient of $\epsilon$, $ad+bc$, is exactly the product rule for its derivative. For a function $f$, we can use a Taylor expansion around $a$ to see that $f(a+\epsilon b)=f(a)+\epsilon b f'(a)$. All higher-order terms contain a factor $\epsilon^2$ and vanish. Therefore, we can compute both the used edge weights and their application to the gradient variables alongside the primal computation.

The main issue with forward-mode differentiation is that the entire derivative computation needs to be carried out separately for each input variable. Similar to finite differences, this does not scale to the large number of parameters encountered in inverse rendering.

Forward AD Code Example

Here is the code for the above example using operator overloading for dual numbers:

import math

class Dual:
    def __init__(self, real, dual):
        self.real = real
        self.dual = dual

    def __add__(self, other):
        # Handle addition with a normal scalar number
        other_real = other.real if isinstance(other, Dual) else other
        other_dual = other.dual if isinstance(other, Dual) else 0.0
        return Dual(self.real + other_real, self.dual + other_dual)

    def __mul__(self, other):
        # Handle multiplication with a normal scalar number
        other_real = other.real if isinstance(other, Dual) else other
        other_dual = other.dual if isinstance(other, Dual) else 0.0

        # (a + eb) * (c + ed) = (ac) + (ad + bc)e
        real = self.real * other_real
        dual = self.real * other_dual + self.dual * other_real
        return Dual(real, dual)

    def __rmul__(self, other):
        # Allows us to do 2 * Dual(...)
        return self.__mul__(other)

    def __repr__(self):
        return f"Dual(val={self.real:.4f}, grad={self.dual:.4f})"

# We define a custom sine function using the dual number Taylor expansion:
# sin(a + eb) = sin(a) + eb * cos(a)
def dual_sin(x):
    if isinstance(x, Dual):
        return Dual(math.sin(x.real), x.dual * math.cos(x.real))
    return math.sin(x)


x = Dual(2.0, dual=1.0) # Seed the derivative here!
y = Dual(3.0, dual=0.0)

# The forward pass
a = 2 * x
b = a * y
c = dual_sin(b)
d = x * x
e = d * c

print("Result:", e)
# Output: Result: Dual(val=-2.1463, grad=18.1062)

Reverse-mode Differentiation

The solution to the forward-mode scaling limitation is to traverse the computation graph in reverse order. Given a sequence of operations, reverse-mode AD will start by evaluating the chain rule for the last operation and proceed toward the input of the algorithm.

Mathematically, reverse-mode computation evaluates vector-Jacobian products (VJP) from the output end of the computation. For a function $f$, it evaluates:

The advantage of this evaluation order is that the gradient computation no longer needs to be duplicated for each input variable. A single reverse traversal computes the gradients for all inputs that affect the chosen output, which is why reverse mode (also known as backpropagation) is effective for optimization problems with millions of parameters.

For the graph $a=2x$, $b=ay$, $c=\sin(b)$, $d=x^2$, and $e=dc$, the complete reverse pass is:

The two terms in $\delta x$ are accumulated because $x$ reaches the output through both $a$ and $d$. Notice that the final line uses $\delta b$: $y$ is an input to $b=ay$, whereas $\delta a$ already encompasses the specific local factor $y$ and corresponds only to the $a$ branch.

Reverse mode is generally more difficult to implement than forward mode. Since it propagates gradients opposite to the primal program’s computation order, it requires storing some of the edge weights of the computation graph in memory to be able to run efficiently.

If we naïvely implemented reverse-mode AD without storing any edge weights, each gradient step would require re-running the primal computation up to the current node. This quadratic complexity is unusable in practice. Conversely, storing the entire graph can easily exceed system memory. The standard remedy is checkpointing, where the program state is only stored at a sparse set of points. As we will see, even checkpointing is often insufficient for physically-based differentiable rendering, requiring more specialized solutions.

❮ ❯

Step 1 / 6

Backward AD Code Example

import math

class Var:
    def __init__(self, val, _children=()):
        self.val = val
        self.grad = 0.0
        # Store the edges of the DAG
        self._prev = set(_children)
        self._backward = lambda: None

    def __mul__(self, other):
        other = other if isinstance(other, Var) else Var(other)
        # Pass (self, other) as children to the new output node
        out = Var(self.val * other.val, (self, other))

        def _backward():
            self.grad += other.val * out.grad
            other.grad += self.val * out.grad
        out._backward = _backward

        return out

    def __rmul__(self, other):
        return self * other

    def backward(self):
        # 1. Topological sort using DFS
        topo = []
        visited = set()

        def build_topo(v):
            if v not in visited:
                visited.add(v)
                for child in v._prev:
                    build_topo(child)
                topo.append(v)

        build_topo(self)

        # 2. Seed the output gradient
        self.grad = 1.0

        # 3. Apply chain rule to the sorted graph in reverse
        for node in reversed(topo):
            node._backward()

    def __repr__(self):
        return f"Var(val={self.val:.4f}, grad={self.grad:.4f})"

def var_sin(x):
    # Pass (x,) as the child
    out = Var(math.sin(x.val), (x,))

    def _backward():
        x.grad += math.cos(x.val) * out.grad
    out._backward = _backward

    return out


x = Var(2.0)
y = Var(3.0)

# Forward pass builds the DAG automatically in the background
a = 2 * x
b = a * y
c = var_sin(b)
d = x * x
e = d * c

# A single call handles the DFS and the entire backward pass
e.backward()

print("Result e:", e)
print("Grad x:", x)
print("Grad y:", y)

# Output:
# Result e: Var(val=-2.1463, grad=1.0000)
# Grad x: Var(val=2.0000, grad=18.1062)
# Grad y: Var(val=3.0000, grad=13.5016)

Why is Differentiable Rendering Difficult?

Symbolically differentiating a Monte Carlo path tracer does not generally work.

As the SIGGRAPH 2020 course notes on physically based differentiable rendering put it:

“Naïve combination of integral discretization and automatic differentiation does not compute the correct derivatives that converge in the limit.” [1]

There are two distinct reasons this happens, and we will look at both in turn: the integrand can be discontinuous in the parameter we are differentiating (the classic visibility problem), or the sampling process used to evaluate the integral can itself depend on that parameter.

The Illegal Swap

The core of the issue lies in the interchange of the derivative and the integral. In calculus, we often assume that:

However, this interchange is only valid under regularity conditions that justify differentiating under the integral sign, such as a suitable integrable bound on $\partial_\pi f$. Parameter-dependent jumps violate these conditions. In rendering, this frequently happens because of visibility: when an object moves, the color changes discontinuously across a moving boundary.

Example 1: Distributional Parameters

When a proposal distribution depends on the differentiated parameter, there are two valid viewpoints: detach both the generated sample and its density, or differentiate both through a parameter-independent primary sample. To compare them, consider estimating the derivative of an integral over an infinite domain.

Why are both correct? In the detached case, $x$ and the evaluated density are constants during differentiation, so

The attached case instead draws $\xi$ from a parameter-independent uniform distribution and differentiates the complete transformed weight. It estimates the same derivative but can have different variance. Bias is introduced by mixing the two viewpoints, for example by detaching $x$ while differentiating only $1/p(x,\lambda)$ and omitting the corresponding change in sampling probability.

Example 2: Discontinuities (The Visibility Problem)

For discontinuous integrands, the fundamental challenge is that the derivative and the integral cannot simply be swapped. Standard Monte Carlo sampling “misses” the boundary contribution entirely.

What goes wrong? The function $f(x, \pi) = (x < \pi\ ?\ 1 : 0.5)$ is a step function: constant everywhere except at the single point $x = \pi$, where it jumps. Any random sample $X$ almost surely lands away from that jump, where the derivative with respect to $\pi$ is exactly zero. The gradient information lives entirely at the moving boundary $x = \pi$, which has probability zero of being hit. So our estimator confidently returns zero every single time, while the true answer is $1/2$.

This is precisely the visibility problem in rendering: when a surface edge moves, the boundary between lit and shadowed regions shifts, but standard path tracing samples almost never land exactly on an edge. The gradient signal is invisible to naive AD.

Code: Generating translation gradient via Finite Difference

scene_path = 'scenes/teapot/teapot.xml'
if not os.path.exists(scene_path): scene_path = 'scenes/teapot/scene.xml'
scene, cam, _ = load_scene_from_xml(scene_path, device=device, override_res=512)

img_base = render_crn(scene, cam, integrator_fd).cpu().numpy()

mesh = scene.get_mesh("teapot")
h = 0.01
mesh.translate([h, 0.0, 0.0])
img_pert = render_crn(scene, cam, integrator_fd).cpu().numpy()
mesh.translate([-h, 0.0, 0.0])

fd_grad = np.mean((img_pert - img_base) / h, axis=-1)
plot_fd(img_base, img_pert, fd_grad, "X-shift", h, vmin=-25, vmax=25)

Motivating Example: Differentiating Visibility

Consider a simplified rendering problem with two constant-color 2D triangles that can occlude each other. The scene parameters are the six triangle vertices ($12$ numbers) and the two RGB colors ($6$ numbers). Given these 18 values as a vector $\boldsymbol{\pi}$, with vertex parameters $\boldsymbol{\pi}_v$ and color parameters $\boldsymbol{\pi}_c$, we want to generate an image $I(\boldsymbol{\pi})$ and compute $\nabla_{\boldsymbol{\pi}} \mathcal{L}(I(\boldsymbol{\pi}))$ for an image-space loss $\mathcal{L}$.

The triangles define an imaging function $m(x,y;\boldsymbol{\pi})$ that maps continuous image coordinates to a color according to the visible triangle. Point sampling this discontinuous function at pixel centers aliases its edges:

Instead, each pixel $I_i$ integrates the imaging function against a reconstruction filter $k$ around its center $(x_i,y_i)$:

The integral changes smoothly as a nondegenerate edge moves, even though its integrand jumps at that edge. We therefore need a differentiation rule that accounts for both changes inside the pixel support and motion of its discontinuity boundary. The next section develops exactly that rule before we return to this scene and implement its gradient.

Mathematical Preliminaries

The Leibniz Integral Rule (1D)

The two failure examples and the triangle scene share one root cause: differentiating only the sampled integrand omits motion of parameter-dependent boundaries. The Leibniz rule provides the formula for differentiating an integral whose limits, as well as its integrand, depend on a parameter $\pi$.

Regularity Conditions

One convenient set of sufficient regularity hypotheses for the Leibniz rule is the following (as detailed in standard real analysis and Delio Vicini’s PhD Thesis):

1. The integration limits $a(\pi)$ and $b(\pi)$ must be continuously differentiable functions of $\pi$.
2. The integrand $f(x, \pi)$ must be differentiable everywhere (specifically continuously differentiable, or $\mathcal{C}^1$) with respect to both $x$ and $\pi$ on the integration domain.
3. Under a measure-theoretic framework (using Lebesgue integration), the partial derivative $\partial f/\partial \pi$ must be Lebesgue-integrable and dominated by a Lebesgue-integrable function (enabling the use of the Lebesgue Dominated Convergence Theorem to swap differentiation and integration in the interior).

Without these hypotheses, for instance if $f$ has interior jump discontinuities that depend on $\pi$, the standard Leibniz rule cannot be applied directly.

For a 1D integral of the form $I(\pi) = \int_{a(\pi)}^{b(\pi)} f(x, \pi) dx$ satisfying these conditions, the derivative is:

Proof

We can derive the general Leibniz rule in two steps: first by assuming constant boundaries, and then generalizing to variable boundaries using the multivariable chain rule.

Part 1: Constant Boundaries

Consider an integral where the limits $a$ and $b$ are constant. We want to find the derivative:

$$\frac{d}{dt} \int_a^b f(t, x) dx$$

Step 1: Expand using the Taylor series We approximate $f(t+\Delta t, x)$ as $f(t, x) + {\color{#00d1b2}\frac{\partial f}{\partial t}\Delta t}$:

$$\frac{\int_a^b \left( f(t, x) + {\color{#00d1b2}\frac{\partial f}{\partial t}\Delta t} \right) dx - \int_a^b f(t, x) dx}{\Delta t}$$

Step 2: Cancel the original function terms The “old” area $\int f dx$ cancels out with the negative term:

$$\frac{\cancel{\int_a^b f(t, x) dx} + \int_a^b {\color{#00d1b2}\frac{\partial f}{\partial t}\Delta t} dx - \cancel{\int_a^b f(t, x) dx}}{\Delta t}$$

Step 3: Cancel the time step $\Delta t$ The $\Delta t$ in the numerator and denominator eliminate each other:

$$\frac{\int_a^b {\color{#00d1b2}\frac{\partial f}{\partial t}\cancel{\Delta t}} dx}{\cancel{\Delta t}} \implies \int_a^b {\color{#00d1b2}\frac{\partial f}{\partial t}} dx$$

Result (Part 1): For constant boundaries, the derivative of the integral is simply the integral of the derivative:

$$\frac{d}{dt} \int_a^b f(t, x) dx = \int_a^b {\color{#00d1b2}\frac{\partial f}{\partial t}} dx$$

Part 2: Variable Limits

Now consider the general case where boundaries depend on time: $I(t) = \int_{a(t)}^{b(t)} f(t, x) dx$.

Step 1: Decompose the integration domain We split the “new” integral $\int_{a+da}^{b+db}$ into the interior $[a,b]$ and the boundary changes:

$$\frac{{\color{#00d1b2}\int_a^b (f + \frac{\partial f}{\partial t}\Delta t) dx} - {\color{#ff6b6b}\int_a^{a + a'\Delta t} (f + \frac{\partial f}{\partial t}\Delta t) dx} + {\color{#4facfe}\int_b^{b + b'\Delta t} (f + \frac{\partial f}{\partial t}\Delta t) dx} - \int_a^b f dx}{\Delta t}$$

Step 2: Discard higher-order terms ($O(\Delta t^2)$) Terms like $\int \frac{\partial f}{\partial t} \Delta t dx$ in the boundary segments (which have width $\approx \Delta t$) become $\Delta t^2$ and vanish:

$$\frac{\int_a^b f dx + {\color{#00d1b2}\int_a^b \frac{\partial f}{\partial t}\Delta t dx} - {\color{#ff6b6b}\int_a^{a + a'\Delta t} f dx} + {\color{#4facfe}\int_b^{b + b'\Delta t} f dx} - \int_a^b f dx}{\Delta t}$$

Step 3: Cancel original function and evaluate boundaries Using the Fundamental Theorem of Calculus (or Mean Value Theorem), the boundary integrals become $f(t, a) a'\Delta t$ and $f(t, b) b'\Delta t$:

$$\frac{\cancel{\int_a^b f dx} + {\color{#00d1b2}\int_a^b \frac{\partial f}{\partial t}\Delta t dx} - {\color{#ff6b6b}f(t, a)a'\Delta t} + {\color{#4facfe}f(t, b)b'\Delta t} - \cancel{\int_a^b f dx}}{\Delta t}$$

Step 4: Final Cancellation of $\Delta t$

$$\frac{{\color{#00d1b2}\int_a^b \frac{\partial f}{\partial t}\cancel{\Delta t} dx} - {\color{#ff6b6b}f(t, a)a'\cancel{\Delta t}} + {\color{#4facfe}f(t, b)b'\cancel{\Delta t}}}{\cancel{\Delta t}}$$

Result (Part 2): For variable boundaries, the derivative of the integral includes the boundary movement terms:

$$\frac{dI}{dt} = {\color{#00d1b2}\int_a^b \frac{\partial f}{\partial t} dx} + {\color{#4facfe}f(t, b) \frac{db}{dt}} - {\color{#ff6b6b}f(t, a) \frac{da}{dt}}$$

Final Result: Combining both parts, we arrive at the general Leibniz Integral Rule (1D):

$$\boxed{\frac{d}{d\pi} \int_{a(\pi)}^{b(\pi)} f(x, \pi) dx = \underbrace{{\color{#00d1b2}\int_{a(\pi)}^{b(\pi)} \frac{\partial f}{\partial \pi}(x, \pi) dx}}_{\text{Interior Term}} + \underbrace{{\color{#4facfe}f(b(\pi), \pi) \frac{db}{d\pi}} - {\color{#ff6b6b}f(a(\pi), \pi) \frac{da}{d\pi}}}_{\text{Boundary Term}}}$$

Generalization: Reynolds Transport Theorem (3D)

In computer graphics, we deal with 2D images and 3D scenes. The 1D Leibniz rule generalizes to higher dimensions via the Reynolds Transport Theorem (RTT).

Regularity Conditions

As in the 1D case, a convenient sufficient set of regularity assumptions for RTT is (see Delio Vicini’s PhD Thesis [2]):

1. Differentiability everywhere in the subdomains: The integrand $f(\mathbf{x}, \pi)$ must be continuously differentiable ($\mathcal{C}^1$) with respect to both $\mathbf{x}$ and $\pi$ everywhere in the interior of the domains separated by the boundary/discontinuity surfaces $\Gamma(\pi)$.
2. Lipschitz Continuity: The boundary motion mapping (the trajectory of boundary points $\mathbf{x}(\pi)$) is Lipschitz continuous, so the boundary velocity field $\partial_\pi \mathbf{x}$ exists almost everywhere.
3. Lebesgue-Integrability: Both the integrand $f(\mathbf{x}, \pi)$ and the partial derivative $\partial_\pi f(\mathbf{x}, \pi)$ must be Lebesgue-integrable over the respective interior domains.

For an integral over a moving domain $X(\pi)$ satisfying these conditions:

Where:

• $X(\pi)$ is the integration domain, which moves as $\pi$ changes.
• $\Gamma(\pi)$ is the full boundary: the union of the external boundary $\partial X(\pi)$ and any internal surfaces where $f$ is discontinuous (e.g. silhouette edges of objects).
• On an internal interface, $\mathbf{n}$ is a consistently chosen unit normal pointing from the minus side to the plus side; on the external boundary it is outward-facing.
• $\partial_\pi \mathbf{x}$ is the velocity of the boundary: how fast each boundary point moves as $\pi$ changes.
• $\Delta f(\mathbf{x}, \pi) = f^-(\mathbf{x}) - f^+(\mathbf{x})$ is the jump in $f$ across $\Gamma$, with $f^-$ on the side from which $\mathbf n$ points and $f^+$ on the side toward which it points. On an external boundary, take the outside value $f^+$ to be zero.

Note that for points on $\Gamma$ where $f$ is actually continuous, $\Delta f = 0$ and they contribute nothing to the boundary integral, so it is safe to include more boundary points than strictly necessary. This matters in practice: when rendering, we do not always know in advance which edges are true silhouettes, so we can include all triangle edges and let the $\Delta f$ term naturally zero out the non-contributing ones.

This is the key formula for differentiable rendering. It tells us that the standard interior derivative must be supplemented with the boundary contribution. As we will see, explicit edge sampling evaluates this term directly; reparameterization and warped-area methods convert it into an equivalent interior estimator.

Continuing Example 2, let’s see how this resolves the failure of naïve AD. The function is:

The discontinuity is at $x = \pi$, so $\Gamma = \{\pi\}$, $\langle \partial_\pi x, \mathbf{n} \rangle = 1$, and the jump is $\Delta f = f^-(\pi) - f^+(\pi) = 1 - 0.5 = 0.5$. Applying the 1D Leibniz rule:

This matches the analytic derivative of $I(\pi) = 0.5\pi + 0.5$, confirming $\frac{dI}{d\pi} = 0.5$. Unlike naive AD, which returns zero by only seeing the interior term, the Leibniz rule correctly captures the contribution of the moving discontinuity by explicitly accounting for the jump $\Delta f$ at the boundary.

Applying the Theory: Two-Triangle Visibility

We now return to the filtered triangle image introduced above. We will not discuss the choice of reconstruction filter $k$ here; the PBRT book provides a detailed treatment of reconstruction filters.

Most renderers, whether real-time, offline, physics-based, differentiable or not, need to deal with the aliasing issue. Most of them solve the antialiasing integral numerically by evaluating the imaging function at sample locations. For a unit-area pixel and uniformly distributed or suitably equidistributed samples, the approximation is:

where $(x_j, y_j)$ are sample locations within the $i$-th pixel. For a nonuniform density $p$, each summand instead carries the importance weight $f(x_j,y_j)/p(x_j,y_j)$. The naive approach of evaluating at the pixel center can also be seen as a one-point quadrature rule with $N = 1$ and $x_1 = y_1 = 0.5$.

We say a discretization is consistent if it converges to the integral, i.e., $\lim_{N\rightarrow \infty} \frac{1}{N} \sum_{j=1}^N f(x_j, y_j; \boldsymbol{\pi}) = I_i$ under the unit-area uniform-sampling convention above. The samples need not be stochastic, but a deterministic sequence must still induce the correct integration measure. If the points are sampled uniformly at random, the estimator is unbiased when $\mathbb{E}[f(x_j, y_j)] = I_i$.

Integration is not limited to antialiasing. Motion blur integrates over the time for which the shutter is open, defocus blur integrates over the lens aperture, and area-light illumination integrates over the light source. The rendering equation similarly expresses global illumination through recursive integration over light-scattering directions.

Remember that our goal is to differentiate a scalar loss $\mathcal{L}$ with respect to the scene-parameter vector $\boldsymbol{\pi}$. The chain rule gives:

Here, $\partial \mathcal{L}/\partial I_i$ measures how the loss responds to pixel $i$, while $\nabla_{\boldsymbol{\pi}} I_i$ measures how that pixel responds to every scene parameter. For a fixed target image $\hat{I}$, a pixel-wise squared loss is

Its gradient is

We therefore need the derivative of each pixel color with respect to the scene parameters.

A common misconception is that a discontinuous visibility function makes the filtered pixel value non-differentiable everywhere. Recall that $I_i$ averages color over the filter support. Away from degenerate events such as topology changes or coincident edges, moving a triangle changes this average smoothly. The rendering integrand can be discontinuous even when its integral is differentiable. Rendering was not turned into an integral merely to obtain this property; image formation is already an integration problem, and rendering algorithms are numerical approximations of that integral.

How do we compute the derivatives of an integral? Recall that we wanted to compute the integral numerically (Equation $\eqref{eq:discretization}$). Unfortunately, we cannot just automatically differentiate the numerical integrator as we saw in Example 2. For vertex-position parameters and samples away from edges, naive AD returns a zero derivative almost surely.

However, the derivative of the integral with respect to a vertex position parameter $\mathbf{\pi}_v$ is not 0.

This is the same failure mode as Example 2: the discretization and the gradient operator do not commute for discontinuous integrands, since a uniformly placed sample has zero probability of landing exactly on the moving edge where the change actually happens. The fix is also the same: sample the boundary explicitly.

In general, we need to evaluate the Reynolds Transport Theorem (Equation $\eqref{eq:reynolds-transport-theorem}$) for this problem:

To intuitively understand the boundary derivative, we can visualize it as calculating the volume of an infinitesimal boundary wedge created by the movement of an edge.

For every point on a silhouette edge, as the parameter $\pi$ changes, the edge sweeps out a small parallelogram. The boundary integral accumulates these infinitesimal volumes along the entire discontinuity contour.

We can decompose the integrand into three intuitive geometric components:

1. Height ($f_- - f_+$): The difference in pixel color (or radiance) between the two sides of the edge (e.g., transitioning from the occluded blue background to the moving red foreground).
2. Width ($n \cdot v$): The distance the edge moves, projected along the normal direction $n$. Movement parallel to the edge simply slides along the boundary and doesn’t change the area; only perpendicular movement contributes to the derivative!
3. Length ($ds$): The differential line element along the boundary contour.

Thus the boundary integral becomes:

which can also be approximated with Monte Carlo sampling.

Complete 2D Triangle Code Implementation (Python)

Code for the above example

The following code is adapted from SIGGRAPH 2020 Course.

import numpy as np

class TriangleMesh:
    def __init__(self, vertices, indices, colors):
        self.vertices = np.array(vertices, dtype=np.float64)  # (N, 2) vertices
        self.indices = np.array(indices, dtype=np.int32)      # (M, 3) face indices
        self.colors = np.array(colors, dtype=np.float64)      # (M, 3) per-face RGB

def raytrace(mesh, pos):
    """
    Uses the half-plane test: a point is inside a triangle if it's
    on the same side of all three edges.
    """
    for i in range(len(mesh.indices)):

        # Extract the current triangle
        idx = mesh.indices[i]
        v0, v1, v2 = mesh.vertices[idx[0]], mesh.vertices[idx[1]], mesh.vertices[idx[2]]

        # Edge normals (2D perpendicular: normal of (dx,dy) = (-dy, dx))
        def normal_2d(v):
            return np.array([-v[1], v[0]])

        # Get edge normals for all edges of triangles
        n01 = normal_2d(v1 - v0)
        n12 = normal_2d(v2 - v1)
        n20 = normal_2d(v0 - v2)

        # Find in which side pos is for each edge
        side01 = np.dot(pos - v0, n01) > 0
        side12 = np.dot(pos - v1, n12) > 0
        side20 = np.dot(pos - v2, n20) > 0

        # if it is on same side for all edges, then it is inside (since this is 2D)
        if (side01 and side12 and side20) or (not side01 and not side12 and not side20):
            return mesh.colors[i], i

    return np.array([0.0, 0.0, 0.0]), -1  # background

def render(mesh, h, w, spp=4):
    """
    Forward pass: render the mesh into an image.
    """
    img = np.zeros((h, w, 3))    # setup the (H, W, 3) buffer for RGB image
    sqrt_spp = int(np.sqrt(spp)) # grid cells for stratified sampling

    # For each pixel
    for y in range(h):
        for x in range(w):
            # for each grid cell
            for dy in range(sqrt_spp):
                for dx in range(sqrt_spp):

                    # Offset the position within the pixel
                    xoff = (dx + np.random.rand()) / sqrt_spp
                    yoff = (dy + np.random.rand()) / sqrt_spp

                    # compute the color at that position
                    pos = np.array([x + xoff, y + yoff])
                    color, _ = raytrace(mesh, pos)
                    img[y, x] += color / spp
    return img

def compute_interior_derivatives(mesh, adjoint, spp=4):
    """
    Interior derivatives: ∂Loss/∂color.
    Standard AD works here because color changes are continuous.
    """
    img_h, img_w = adjoint.shape[:2]
    sqrt_spp = int(np.sqrt(spp))
    d_colors = np.zeros_like(mesh.colors)

    # For each pixel
    for y in range(img_h):
        for x in range(img_w):
            # For each grid cell
            for dy in range(sqrt_spp):
                for dx in range(sqrt_spp):

                    # Find the position within the cell within pixel
                    xoff = (dx + np.random.rand()) / sqrt_spp
                    yoff = (dy + np.random.rand()) / sqrt_spp

                    # compute the gradient at that position
                    pos = np.array([x + xoff, y + yoff])
                    _, hit_idx = raytrace(mesh, pos)
                    if hit_idx >= 0:
                        d_colors[hit_idx] += adjoint[y, x] / spp
    return d_colors


def collect_edges(mesh):
    """Collect unique edges."""
    edges = set() # Stores edges as tuples (u, v)

    for idx in mesh.indices:
        edges.add((min(idx[0], idx[1]), max(idx[0], idx[1])))
        edges.add((min(idx[1], idx[2]), max(idx[1], idx[2])))
        edges.add((min(idx[2], idx[0]), max(idx[2], idx[0])))

    # [(u, v) ...]
    return list(edges)

def build_edge_sampler(mesh, edges):
    """Build CDF for importance-sampling edges by length."""
    lengths = []

    # Store the lengths of the edges
    for v0_id, v1_id in edges:
        lengths.append(np.linalg.norm(mesh.vertices[v1_id] - mesh.vertices[v0_id]))

    lengths = np.array(lengths)

    # Use the edge lengths as weight for PDF and construct CDF
    pmf = lengths / lengths.sum()
    cdf = np.concatenate([[0], np.cumsum(pmf)])

    return pmf, cdf, lengths

def compute_edge_derivatives(mesh, adjoint, n_edge_samples=10000):
    """∂Loss/∂vertices via Reynolds Transport Theorem."""

    # Extract unique edges and build CDF for sampling
    img_h, img_w = adjoint.shape[:2]
    edges = collect_edges(mesh)
    pmf, cdf, lengths = build_edge_sampler(mesh, edges)

    d_vertices = np.zeros_like(mesh.vertices)
    screen_dx = np.zeros((img_h, img_w, 3))
    screen_dy = np.zeros((img_h, img_w, 3))

    for i in range(n_edge_samples):
        # 1. Pick an edge (importance sampling by length)
        u = np.random.rand()
        edge_id = np.searchsorted(cdf, u, side='right') - 1
        edge_id = np.clip(edge_id, 0, len(edges) - 1)
        u, v = edges[edge_id]

        # 2. Pick a point on the edge
        v0 = mesh.vertices[u]
        v1 = mesh.vertices[v]
        t = np.random.rand()   # t in [0, 1]
        p = v0 + t * (v1 - v0)

        xi, yi = int(p[0]), int(p[1])
        if xi < 0 or yi < 0 or xi >= img_w or yi >= img_h:
            continue

        # 3. Sample both sides of the edge (the "jump" / discontinuity)
        edge_dir = (v1 - v0) / np.linalg.norm(v1 - v0)
        n = np.array([-edge_dir[1], edge_dir[0]])  # outward normal
        eps = 1e-3

        color_in, _ = raytrace(mesh, p - eps * n)
        color_out, _ = raytrace(mesh, p + eps * n)

        # 4. Compute gradient contribution (Reynolds Transport Theorem)
        pdf = pmf[edge_id] / lengths[edge_id]
        weight = 1.0 / (pdf * n_edge_samples)
        color_diff = color_in - color_out  # the jump Δf
        adj = np.dot(color_diff, adjoint[yi, xi])

        # dp/dv0 = (1-t), dp/dv1 = t  (from p = v0 + t*(v1-v0))
        d_v0 = np.array([(1 - t) * n[0], (1 - t) * n[1]]) * adj * weight
        d_v1 = np.array([t * n[0], t * n[1]]) * adj * weight

        d_vertices[u] += d_v0
        d_vertices[v] += d_v1

        # Screen-space derivatives
        screen_dx[yi, xi] += -n[0] * color_diff * weight
        screen_dy[yi, xi] += -n[1] * color_diff * weight

    return d_vertices, screen_dx, screen_dy


# 1. Scene setup
c_blue =[15/255, 133/255, 165/255]
c_red  =[187/255, 37/255, 66/255]

scale = 2.0
mesh = TriangleMesh(
    vertices = np.array([
        # Tri 0 (Red)
        [10.0, 12.0],[26.0, 1.0], [31.0, 16.0],
        # Tri 1 (Blue)
        [2.0, 11.0],[16.0, 2.0], [20.0, 19.0],
    ]) * scale,
    indices = [[0, 1, 2], [3, 4, 5]],
    colors =[c_red, c_blue]
)

# Window setup
W, H, spp = 70, 45, 4
np.random.seed(48)

# 2. Forward Pass
print("Rendering...")
img = render(mesh, H, W, spp)

# 3. Backward Pass (Interior: ∂I/∂color)
adjoint = np.ones((H, W, 3)) # Uniform adjoint to pull gradients
d_colors = compute_interior_derivatives(mesh, adjoint, spp)

# 4. Backward Pass (Edges: ∂I/∂vertex via boundary sampling)
d_verts, screen_dx, screen_dy = compute_edge_derivatives(mesh, adjoint, n_edge_samples=W*H)

print("\nVertex Gradients (d_verts):")
print(np.round(d_verts, 4))

# Output:
# Vertex Gradients (d_verts):
# [[ -4.2248   2.533 ]
#  [  7.4785 -18.8305]
#  [ 13.7454  13.4763]
#  [-21.0542   4.3572]
#  [  0.4232 -20.9386]
#  [  2.0691  19.6481]]

While explicitly finding and sampling edges works well for 2D triangles, doing this for complex 3D meshes with secondary bounces such as shadows and reflections is much harder. We now turn to estimators designed for full Monte Carlo light transport.

Differentiable Monte Carlo Rendering

Our goal is to solve inverse rendering problems of the form:

where $g$ is an image-based objective function. To simplify the notation, we will consider only the intensity $I$ of a single pixel $j$ and one differentiable parameter $\pi$. The derivations generalize to differentiable rendering of RGB images and multiple parameters. (Note: from the estimator derivations onward we upgrade $\pi$ to the vector $\boldsymbol{\pi}$, to represent gradients with respect to the entire scene parameter space simultaneously.)

Objective Function Gradient

As in the two-triangle example above, the outermost step is just the chain rule. What’s new this time is that $I$ itself will be estimated by noisy Monte Carlo samples rather than computed exactly, and we need to handle that carefully. Using the simplified notation, our goal is to compute the derivative $\partial_\pi g(I(\pi))$. The chain rule allows writing this term as:

where $g'$ is the derivative of the objective function. We further declutter the notation by dropping the explicit dependency of $I$ on $\boldsymbol{\pi}$ from now on. We use Monte Carlo integration to estimate $I$. If we replace $I$ with a Monte Carlo estimator $\hat{I}$ in the equation above and take the expected value we get:

Generally, this is not an unbiased estimator of the true objective gradient. One source of bias is that $g'(\hat{I})$ and $\partial_{\boldsymbol{\pi}} \hat{I}$ use the same random samples, producing the covariance term. We can remove that covariance term, provided the random streams are independent, by using a primal estimator $\hat{I}^p$ for $g'$ and a separate derivative estimator $\partial_{\boldsymbol{\pi}} \hat{I}^a$:

In practice, this means rendering two images with independent random number streams. For nonlinear $g$, a second plug-in bias can remain because $\mathbb{E}[g'(\hat I^p)]$ need not equal $g'(I)$; increasing the primal sample count reduces this bias.

Detached Estimator

The remaining challenge is to estimate $\partial_{\boldsymbol{\pi}} I$ itself. The following derivation assumes that the integrand has no $\boldsymbol{\pi}$-dependent discontinuities. Mathematically, we need to differentiate a parameter-dependent, high-dimensional integral over light paths:

Here, a path $\mathbf{x}=(\mathbf{x}_0,\ldots,\mathbf{x}_k)$ is a sequence of sensor, surface, and emitter vertices, and $\mathcal{P}$ denotes the union of these path spaces over possible lengths $k$. A renderer also makes discrete choices, including path length, light or BSDF lobe selection, and Russian roulette; we return to those choices below. The function $f$ is the parameter-dependent contribution of a path. If $f$ does not contain parameter-dependent discontinuities, we can directly estimate its derivative using Monte Carlo integration. The derivative operator can be moved into the integral:

For this estimator, we need to differentiate the evaluation of $f$. We do not have to differentiate the sampling process that produces $\mathbf{x}_i$ or the corresponding PDF $p(\mathbf{x}_i)$. We call this estimator detached since both sampling and PDF evaluation are detached from the differentiation process. This is the most commonly used estimator in differentiable rendering. Zeltner et al. (2021) [12] provide the systematic study of this attached/detached distinction that the next two subsections summarize (see Fig. for the overall taxonomy).

If $f$ contains $\boldsymbol{\pi}$-dependent discontinuities, additional precautions are required (e.g., edge sampling or reparameterization). In this case, the detached estimator captures only the interior term of the Reynolds Transport Theorem; the boundary integral from moving discontinuities must be estimated separately and added to obtain the full derivative. Similarly, if the path space $\mathcal{P}$ is parameter-dependent, we need to account for changes in its geometry or switch to a parameterization of the integration domain that is independent of $\boldsymbol{\pi}$.

Attached Estimator

While conceptually simple, the detached estimator does not handle all potential use cases. In particular, it does not support perfectly specular BSDFs. Such BSDFs are delta functions, which do not yield valid derivatives. The solution to this problem is to also differentiate the BSDF sampling process. By doing so, we switch from differentiating the integrand by itself to differentiating the ratio of integrand to PDF. This avoids having to differentiate the delta function of the specular BSDF, as it cancels out with the sampling density.

Differentiating the sampling process can be interesting beyond perfectly specular surfaces. Many of the sampling steps in a Monte Carlo renderer are highly scene-dependent. For example, the roughness parameter of a microfacet BSDF will affect the sampling of the scattered direction. This and other sampling methods usually transform a set of uniformly distributed random numbers to the desired target distribution, e.g., using inverse transform sampling. We can interpret this transformation as a reparameterization of the original integral.

Because the sampling strategy may produce different distributions depending on the parameter $\boldsymbol{\pi}$, this choice also affects the variance properties of the resulting gradient estimator.

Formally, sampling strategies can be understood as a change of variables to new coordinates $\mathbf{u} \in \mathcal{U}$ parameterizing the integration domain $\mathcal{P}$ via a mapping $\mathcal{T} : \mathcal{U} \to \mathcal{P}$, where $\mathcal{U} =[0, 1]^n$ is a unit-sized hypercube of suitable dimension. The space $\mathcal{U}$ is called the primary sample space. The mapping $\mathbf{x} = \mathcal{T}(\mathbf{u})$ is constructed from a target density $p(\mathbf{x})$ so that its Jacobian determinant satisfies $|J_\mathcal{T}(\mathbf{u})| = p(\mathbf{x})^{-1}$. The reparameterized integral then takes the form:

This formulation is called attached, since samples geometrically follow the motion of $\mathcal{T}(\mathbf{u}, \boldsymbol{\pi})$ with respect to perturbations of $\boldsymbol{\pi}$. Similar to before, we can build an estimator of the derivative by applying Monte Carlo integration:

The attached estimator is primarily useful for perfectly specular surfaces, but it can also produce lower variance than the detached version for derivatives of BSDFs with low roughness. On the other hand, the additional motion of the samples might introduce more variance in the evaluation of other terms in the integrand.

Attached sampling handles a local delta interaction when the sampled specular direction changes smoothly with the scene parameters. It does not by itself resolve discontinuous visibility through a chain of specular events or changes in caustic-path topology. Those cases require specialized path-space or manifold techniques and are outside the surface-visibility methods developed here.

Finally, the attached estimator is more difficult to use as in practice it requires handling discontinuities in the sampling function $\mathcal{T}$. Examples of such discontinuities are discrete sampling decisions (such as BSDF component selection) or discontinuities due to sampled rays hitting different objects as $\boldsymbol{\pi}$ changes.

Discrete Sampling Decisions and MIS

The attached formula assumes a differentiable map $\mathcal{T}$, but practical path tracers also make discrete decisions: selecting a light or BSDF lobe, accepting a Russian-roulette continuation, or choosing among multiple importance sampling (MIS) techniques. A branch selected by a Bernoulli or categorical sample is locally constant, so ordinary pathwise AD cannot differentiate the change in its probability.

Russian roulette gives a useful example. If a path survives with probability $q(\pi)$, its surviving contribution is divided by $q(\pi)$. Differentiating the factor $1/q$ while treating the sampled survive/terminate decision as constant omits the derivative of the decision probability and is generally biased. Two consistent options are common:

1. Detach the proposal decision and its compensation. Sample survival using the current $q$, but stop gradients through both the discrete decision and $q$ in the Monte Carlo weight. The resulting detached estimator differentiates the underlying transport contribution rather than the proposal mechanism.
2. Differentiate the probability consistently. Add the corresponding score-function term, or use a valid continuous reparameterization when one exists. This is usually more expensive and can have high variance.

The same rule applies to light and lobe selection. MIS adds another layer because its weights depend on the PDFs of several techniques. Proposal PDFs and MIS weights should not be differentiated selectively: derive the complete estimator as either detached or attached, then apply that choice consistently to sampling, PDF factors, and weights. Selectively differentiating a PDF denominator or MIS weight while detaching the random choice that produced it is the mixed failure mode described in Example 1 (see Fig. ).

Explicit Edge Sampling (Li et al.)

While the interior term is straightforward to evaluate with the differentiable Monte Carlo estimators introduced above, the boundary integral poses a greater challenge: derivatives arising from visibility discontinuities must either be integrated explicitly over silhouette edges or reformulated as an equivalent smooth-domain integral.

Li et al. (2018) [10] model visibility with Heaviside step functions. Differentiating a step function yields a Dirac delta concentrated on the moving edge, so their estimator naturally decomposes the image derivative into two parts: the smooth interior term, handled by standard Monte Carlo sampling with AD, and the singular boundary term, estimated by a dedicated edge sampler. This decomposition is the distributional counterpart of the Reynolds Transport Theorem split derived above.

Primary Visibility

We begin with the $2D$ pixel filter integral, which for each pixel integrates the pixel filter $k$ against the incoming radiance $L$. The radiance itself may be a further integral over light sources or the hemisphere, but for convenience we absorb everything into a single scene function $f(x,y) = k(x,y)L(x,y)$, the $f$ used throughout the remainder of this section. The pixel color $I$ is then:

Under the paper’s assumptions of non-interpenetrating triangle meshes, finite-area emitters, non-delta BSDFs, and static scenes, the relevant visibility discontinuities occur at projected triangle edges. This makes it possible to integrate over them explicitly. Li et al. were the first to systematically study these discontinuities in the context of differentiable rendering, proposing Monte Carlo integration of the boundary term by directly sampling the edges responsible for visibility jumps. Open boundary edges, view-dependent silhouette edges, and sharp edges where neighbouring faces have differing normals can all define discontinuities; with smooth shading, only edges across which the scene function actually jumps contribute to the boundary estimator.

A $2D$ triangle edge partitions the domain into two half-spaces, $f_u$ and $f_l$ (illustrated below). The discontinuity across the edge can be modelled with the Heaviside step function $\theta$:

where $f_u$ represents the upper half-space, $f_l$ represents the lower half-space, and $\alpha$ defines the edge equation formed by the triangles. For each edge with endpoints $\mathbf{a} = (a_x, a_y)$ and $\mathbf{b} = (b_x, b_y)$, we construct the edge equation $\alpha(x, y) = Ax + By + C$. Since $\alpha(x, y) = 0$ along the line passing through both endpoints, substituting $\mathbf{a}$ and $\mathbf{b}$ yields:

A scene function $f$ can be rewritten as a sum of such Heaviside functions $\theta$, one per edge, and $f_i$ itself can contain further nested Heaviside terms (a single triangle is the product of three Heaviside step functions). This fact is also crucial for generalization to secondary visibility.

Analytically differentiating a single term in this sum, we apply the product and chain rules step by step:

The second (interior) term simply replaces $f_i$ with its gradient, which automatic differentiation handles directly. All of the new machinery developed below targets the first (boundary) term.

Because $\delta(\alpha_i(x,y))$ is nonzero only on the curve $\{\alpha_i(x,y) = 0\}$, the edge itself, the 2D area integral collapses to a 1D integral along that curve. This converts the Dirac delta into an ordinary arc-length integral over the edge:

2D Area Integral to 1D Arc-Length Line Integral

Assume that $\alpha \in C^1$ and that $\nabla\alpha \neq 0$ on the regular level set

$$ \mathcal{C}=\{(x,y)\mid \alpha(x,y)=0\}. $$

Near any point on $\mathcal{C}$, introduce local coordinates $(n,\sigma)$, where $n$ is the signed distance measured along the unit normal

$$ \mathbf n=\frac{\nabla\alpha}{\|\nabla\alpha\|}, $$

and $\sigma$ is the arc-length parameter along the curve.

Along the normal direction,

$$ \frac{\partial\alpha}{\partial n} =\nabla\alpha\cdot\mathbf n =\|\nabla\alpha\|, $$

so that

$$ d\alpha=\|\nabla\alpha\|\,dn, \qquad dn=\frac{d\alpha}{\|\nabla\alpha\|}. $$

Since the signed-distance coordinate has unit Jacobian on the curve, the area element becomes

$$ dx\,dy =dn\,d\sigma =\frac{d\alpha}{\|\nabla\alpha(x,y)\|}\,d\sigma. $$

Substituting this change of variables into the area integral gives

$$ \iint_{\mathbb R^2} \delta(\alpha)\,\mathbf g\;dx\,dy = \int_{\mathcal C} \left( \int_{-\infty}^{\infty} \delta(\alpha) \frac{\mathbf g(x,y)}{\|\nabla\alpha(x,y)\|} \,d\alpha \right) d\sigma. $$

Applying the one-dimensional sifting property of the Dirac delta, $\int_{-\infty}^{\infty}\delta(\alpha)h(\alpha)\,d\alpha =h(0)$, yields

$$ \iint_{\mathbb R^2} \delta(\alpha(x,y)) \,\mathbf g(x,y)\,dx\,dy = \int_{\alpha=0} \frac{\mathbf g(x,y)}{\|\nabla\alpha(x,y)\|} \,d\sigma. $$

Finally, setting $\mathbf g(x,y)=\nabla\alpha_i(x,y)\,f_i(x,y)$ immediately gives Equation $\eqref{eq:2d-delta-to-arclength}$.

The gradients of the edge equations $\alpha_i$ are:

Derivation of Edge Equation Partial Derivatives

The implicit edge equation $\alpha(x,y)$ for the line connecting endpoints $\mathbf{a} = (a_x, a_y)$ and $\mathbf{b} = (b_x, b_y)$ is given by the cross product determinant:

$$ \alpha(x,y) = (a_y - b_y)x + (b_x - a_x)y + (a_x b_y - b_x a_y). $$

Differentiating $\alpha(x,y)$ with respect to each variable step by step:

1. Spatial Derivatives ($x, y$)

   $$ \frac{\partial\alpha}{\partial x} = a_y - b_y, \qquad \frac{\partial\alpha}{\partial y} = b_x - a_x $$

   The norm of the spatial gradient equals the length of the edge segment:

   $$ \begin{aligned} \lVert \nabla_{x,y}\alpha \rVert &= \sqrt{\left(\frac{\partial\alpha}{\partial x}\right)^2 + \left(\frac{\partial\alpha}{\partial y}\right)^2} \\ &= \sqrt{(a_y - b_y)^2 + (b_x - a_x)^2} \\ &= \sqrt{(a_x - b_x)^2 + (a_y - b_y)^2}. \end{aligned} $$

2. Vertex $\mathbf{a}$ Derivatives ($a_x, a_y$)

   $$ \begin{aligned} \frac{\partial\alpha}{\partial a_x} &= \frac{\partial}{\partial a_x}\big[(a_y - b_y)x + (b_x - a_x)y + (a_x b_y - b_x a_y)\big] \\ &= -y + b_y = b_y - y \\[0.6em] \frac{\partial\alpha}{\partial a_y} &= \frac{\partial}{\partial a_y}\big[(a_y - b_y)x + (b_x - a_x)y + (a_x b_y - b_x a_y)\big] \\ &= x - b_x. \end{aligned} $$

3. Vertex $\mathbf{b}$ Derivatives ($b_x, b_y$)

   $$ \begin{aligned} \frac{\partial\alpha}{\partial b_x} &= \frac{\partial}{\partial b_x}\big[(a_y - b_y)x + (b_x - a_x)y + (a_x b_y - b_x a_y)\big] \\ &= y - a_y \\[0.6em] \frac{\partial\alpha}{\partial b_y} &= \frac{\partial}{\partial b_y}\big[(a_y - b_y)x + (b_x - a_x)y + (a_x b_y - b_x a_y)\big] \\ &= -x + a_x = a_x - x. \quad \end{aligned} $$

Gradients with respect to other scene parameters, such as camera pose, 3D vertex positions, or vertex normals, follow by applying the chain rule through the projection of the triangle vertices:

The image function across an edge $\alpha_i = 0$ is a step combination of the upper/interior color $f_u$ and lower/exterior color $f_l$:

Differentiating yields $\nabla \theta(\alpha_i) = \delta(\alpha_i) \nabla \alpha_i$, which isolates the color jump $\Delta f = f_u - f_l$ across the boundary.

To evaluate the 1D boundary integral over an edge $E$ with endpoints $\mathbf{a}$ and $\mathbf{b}$, we reparameterize the arc length via a line parameter $t \in [0,1]$ using $(x(t), y(t)) = (1-t)\mathbf{a} + t\mathbf{b}$. Since $\mathrm{d}\sigma = \lVert\mathbf{b}-\mathbf{a}\rVert\,\mathrm{d}t = \lVert E \rVert\,\mathrm{d}t$, the integral transforms as:

To estimate this integral via Monte Carlo, we first select a candidate edge $E$ from the scene with a discrete probability $P(E)$. We then draw a sample point $\mathbf{x}_j = (x_j, y_j)$ uniformly along the length of $E$.

The marginal probability density $p(\mathbf{x}_j)$ of selecting a specific boundary point $\mathbf{x}_j$ is evaluated using the law of total probability over the set of all edges $\mathcal{E}$:

(Note: we safely ignore the measure-zero set of vertices where edges intersect).

Substituting this probability density $p(\mathbf{x}_j)$ into the standard primary Monte Carlo estimator $\frac{1}{N}\sum_{j=1}^N \frac{f(\mathbf{x}_j)}{p(\mathbf{x}_j)}$ yields the unbiased boundary estimator:

where $\lVert E \rVert$ is the projected screen-space length of edge $E$, and $P(E)$ is the discrete probability of selecting that edge.

In practice, if we employ smooth shading, most triangle edges lie entirely in continuous regions where, by definition of continuity, $f_u(x,y) = f_l(x,y)$. The Dirac integral is therefore zero for these internal edges, and only silhouette edges have non-zero contributions to the gradients.

Given a camera viewpoint $\mathbf{v}$ and an edge associated with two adjacent faces with normals $\mathbf{n}_f$ and $\mathbf{n}_b$, the edge is a silhouette if for any point $\mathbf{p}$ on it, the view vector $\mathbf{p} - \mathbf{v}$ faces opposite directions with respect to the two normals:

To select the edges, all triangle meshes are projected into screen space and clipped against the camera frustum. We then select a silhouette edge $E$ with probability proportional to its screen-space length:

and uniformly pick a sample point $\mathbf{x}_j$ along the selected edge. Notice that choosing $P(E) \propto \lVert E \rVert$ neatly cancels the length term $\lVert E \rVert$ in the Monte Carlo weight $\frac{\lVert E \rVert}{P(E)}$, leaving simply the sum of all projected silhouette lengths.

This method handles occlusion naturally. If a sample is blocked by another surface, $(x, y)$ lands on the continuous part of the contribution function $f(x, y)$, producing $f_u(x, y) = f_l(x, y)$. Such samples evaluate to zero in the difference and contribute nothing to the gradient (Fig. b).

Secondary Visibility

This method can be generalized to handle shadows, reflections, and indirect illumination by integrating over the $3D$ scene.

Similar to the primary visibility case, an edge $(v_0, v_1)$ in 3D introduces a step function into the scene function $h$:

The 3D edge function $\alpha(m)$ is obtained by constructing a plane through the shading point $p$ and the two edge vertices. The sign of the dot product of $m - p$ with the plane normal assigns each point to one of the two half-spaces. Concretely, the edge equation is defined as

The gradient computation follows the same derivation as primary visibility, now applying the 3D counterparts of $\eqref{eq:2d-edge-derivation}$ and $\eqref{eq:2d-delta-to-arclength}$ with $x, y$ replaced by $p, m$. The resulting edge integral, the scene-surface analogue of the screen-space boundary integral, is:

where $n_m$ is the surface normal at $m$. Two key differences distinguish this 3D integral from its screen-space counterpart. First, the measure $\sigma'(m)$ is no longer the arc length along the 2D edge; instead it measures the projected length from the edge through the shading point $p$ onto the scene manifold (the semi-transparent triangle in Fig. (a) illustrates this projection). Second, an additional area-correction factor $\lVert n_m \times n_h \rVert$ appears because the scene surface element must be projected onto the infinitesimal width of the edge (Fig. (b)).

Here the Jacobian $J_m(t)$ is a 3D vector that captures how the edge $(v_0, v_1)$ projects onto the scene manifold as a function of the line parameter.

Derivation of the 3D edge Jacobian $J_m(t)$

We derive the Jacobian $J_m(t)$ in Equation $\eqref{eq:3d-edge-integral}$. The goal is to compute the derivatives of point $m(t)$ with respect to the line parameter $t$. The relation between $m(t)$ and $t$ is described by a ray-plane intersection. That is, we are intersecting a plane at point $m$ with normal $n_m$ with a ray of origin $p$ and unnormalized direction $\omega(t)$:

$$\omega(t) = v_0 + (v_1 - v_0)t - p$$$$\tau(t) = \frac{(m - p) \cdot n_m}{\omega(t) \cdot n_m}$$$$m(t) = p + \tau(t)\omega(t).$$

We can then derive the derivative $J_m(t) = \frac{\partial m(t)}{\partial t}$ as:

$$J_m(t) = \tau(t)\left((v_1 - v_0) - \omega(t)\frac{(v_1 - v_0) \cdot n_m}{\omega(t) \cdot n_m}\right)$$

The partial derivatives of $\alpha(p, m)$ required by the edge integral are:

These are the corrected forms from the paper’s published erratum. In particular, $p$ occurs in all three factors of the scalar triple product, so its derivative is not equal to $\nabla_m\alpha$.

Unlike primary visibility where the camera viewpoint is fixed, the shading point $\mathbf{p}$ varies across bounces throughout the scene. Consequently, an acceleration structure is required to prune non-contributing silhouette edges efficiently.

Selecting an Edge

Explicit edge sampling is attractive for primary visibility because the camera is fixed and projected silhouettes can be precomputed. Secondary visibility is harder: the shading point changes at every path vertex, and performance degrades with geometric and depth complexity. Containing that cost is what the acceleration structures below are for.

Three factors govern the importance of an edge at a given shading point: the geometric foreshortening (proportional to inverse squared distance to the edge), the material response between the shading point and the point on the edge, and the incoming radiance from the edge direction (e.g. whether it hits a light source or not).

To importance-sample edges efficiently, Li et al. build two acceleration hierarchies:

3D BVH for triangle edges associated with only one face (boundary edges) and meshes without smooth shading normals, built from the 3D positions of each edge’s two endpoints.

6D BVH for the remaining edges, built from the two endpoint positions and the two normals of the adjacent faces.

Each hierarchy node stores a cone direction and opening angle covering all possible normal directions within it, enabling quick rejection of non-silhouette edges. The directional components are scaled by $\frac{1}{8}$ the diagonal of the scene’s bounding box, and during construction the node is split along the dimension with the longest extent.

The hierarchy is traversed twice. The first traversal focuses on edges that overlap with the cone subtended by the light source at the shading point, using a box-cone intersection to quickly discard edges that do not intersect the light sources. The second traversal samples all edges. The two sets of samples are combined using multiple importance sampling.

During traversal, for each node an importance value is computed by upper-bounding the contribution: total edge length $\times$ inverse squared distance $\times$ a Blinn-Phong BRDF bound. Nodes that do not contain any silhouette receive zero importance. Both children are traversed if the shading point lies inside both bounding boxes, the BRDF bound exceeds a threshold (set to $1$), or the angle subtended by the light cone is smaller than $\cos^{-1}(0.95)$.

Importance Sampling a Single Edge

Once an edge is selected, a point along it must be chosen. With a highly specular BRDF, only a small portion of the edge carries significant contribution. The Linearly Transformed Cosine (LTC) distribution provides a closed-form solution for the integral between a point and a linear light source, weighted by BRDF and geometric foreshortening. The integrated CDF is numerically inverted via Newton’s method for importance sampling, using a precomputed table of fitted LTC lobes for the target BRDFs.

Compared to the baseline of uniformly sampling edges by length, this importance sampling strategy is far more effective at capturing rare events (shadows cast by a small light source or very specular reflections of edges) and produces images with much lower variance. The problem is structurally similar to next-event estimation with many light sources, where the set of important sources depends on the current shading point.

Limitations

Performance. Explicit edge sampling is expensive, especially for secondary visibility and does not scale well to complex scenes with many edges.

Other light transport phenomena. As noted above, the method assumes static scenes with no participating media. Differentiating motion blur requires sampling on 4D edges with an extra time dimension.

Interpenetrating geometries and parallel edges. Dealing with the derivatives of interpenetration of triangles requires a mesh splitting process and its derivatives. Interpenetration can happen if the mesh is generated by some simulation process. This method also does not handle the case where two edges are perfectly aligned as seen from the center of projection (camera or shadow ray origin). However, these are zero-measure sets in path space, and as long as the two edges are not perfectly aligned to the viewport, we will be able to converge to the correct solution.

Shader discontinuities. The method assumes the BSDF models and shaders are differentiable and does not handle discontinuities in the shaders. They also don’t handle the discontinuities at total internal reflection and some other BRDFs relying on discrete operations.

Edge sampling sets the template that the rest of this section varies: find the boundary, sample it, add its contribution. One direction is to keep that template and make the search better. Zhang et al. (2020) [9] lift it into path space, sampling points and directions on silhouette edges and expanding them outwards into complete paths connecting sensor and emitter, which is harder to implement but produces high-quality edge gradients under difficult lighting. The other direction, taken next, is to stop searching for the boundary altogether.

Reparameterizing Visibility Discontinuities (Loubet et al.)

Explicit edge sampling does not always scale efficiently and is difficult to generalize to implicit surface representations, where discontinuities are not simply a discrete set of mesh edges.

Loubet et al. [11] instead apply a change of variables that removes or reduces the parameter-dependence of discontinuity locations. If the transformation fixes every moving discontinuity exactly, the derivative operator can be moved inside the transformed integral and accounting for its Jacobian yields an unbiased gradient estimator. Their practical rotations approximate this ideal transformation, however, so the paper describes the resulting gradients as low-bias rather than unbiased.

Given a transformation $\mathcal{T}:\mathcal{Y}\rightarrow\mathcal{X}$, the reparameterized integral

can be differentiated under the integral sign in the usual way.

1D Example

Consider differentiating the integral derivative:

The integrand has a step function at position $\pi$ multiplied by a smooth function $g$. The step discontinuity prevents moving $\partial_\pi$ inside the integral. Substituting $y = x - \pi$ (with unit Jacobian $|\operatorname{det} J_\mathcal{T}| = 1$) yields:

The indicator $\mathbb{1}_{[0, \infty]}(y)$ no longer depends on $\pi$, making the integrand differentiable with respect to $\pi$ for almost every fixed sample $y$:

Two equivalent interpretations:

• Instead of integrating a function with a moving discontinuity, we integrate in a reparameterized domain where the discontinuity location is fixed.
• This is equivalent to importance sampling $\int f(x) \, \mathrm{d}x$ using samples $x_i(\pi) = y_i + \pi$ that follow the movement of the discontinuity.

To preserve the primal computation of $I$, the transformation $\mathcal{T}$ should be the identity map at the current parameter value $\pi_0$, i.e., $\mathcal{T}(y, \pi) = y + \pi - \pi_0$. The step location is fixed at $y = \pi_0$, allowing automatic differentiation to evaluate the smooth motion of $g$ without differentiating through a moving visibility test. Note that the sampling density $p(y_i)$ must not depend on $\pi$, otherwise parameter dependencies are reintroduced into the integrand.

Removing Discontinuities Using Rotations

A typical shading integral can contain complex parameter-dependent discontinuities. However, when the integrand has small angular support (e.g., narrow pixel reconstruction filters, glossy BSDF lobes, or small light sources), the discontinuity within the support reduces to the silhouette of a single object, as shown in Fig. .

The displacement of a silhouette on $\mathbb{S}^2$ under infinitesimal perturbations of $\pi$ is well approximated by a spherical rotation (the spherical counterpart of a planar domain translation). As the support shrinks, this approximation improves, becoming exact in the limit. Assuming a suitable rotation $R(\boldsymbol{\omega}, \pi)$ exists, the change of variables

makes $f(R(\boldsymbol{\omega}, \pi), \pi)$ continuous with respect to $\pi$ for each direction $\boldsymbol{\omega}$. The rotation determinant is $|\operatorname{det} J_R| = 1$, and $R$ depends explicitly on $\pi$.

Rotations are simple to compute and accurately track local boundary movements. Using $R$ to reparameterize the integral yields the Monte Carlo estimator:

where $\boldsymbol{\omega}_i \sim p(\cdot, \pi_0)$ are drawn from the default sampling distribution (e.g., BSDF sampling) evaluated at $\pi_0$ rather than $\pi$, removing sample dependency on $\pi$.

Generalizing to Functions with Large Support

When integrands have large support on $\mathbb{S}^2$, they contain multiple interacting silhouettes that violate the single-object assumption, causing bias in local rotation estimates.

To resolve this, we leverage the property that the integral of a function $f$ equals the integral of its spherical convolution:

where $k$ is a normalized spherical convolution kernel. With the argument order used above, preservation of the integral requires:

By choosing $k$ to be a smooth, concentrated distribution (such as a von Mises-Fisher distribution) with small angular support, the inner integral is restricted to a small domain, restoring compatibility with local rotations.

von Mises-Fisher (vMF) Distribution & Sampling

The von-Mises-Fisher (vMF) distribution is commonly used to describe directional data and can be regarded as akin to the isotropic Gaussian distribution on the sphere in $D$ dimensions, $\mathbb{S}^{D-1}$. It is parametrized by a mean direction $\mu \in \mathbb{S}^{D-1}$ and a concentration $\tau > 0$. Its density is defined as:

$$ \text{vMF}(n; \mu, \tau) = Z(\tau) \exp(\tau \mu^T n) , \quad Z(\tau) = (2\pi)^{-D/2} \frac{\tau^{D/2-1}}{\mathrm{I}_{D/2-1}(\tau)} , $$

where $\mathrm{I}_\nu$ is the modified Bessel function of the first kind of order $\nu$. In $D = 3$ dimensions, the normalizer simplifies to:

$$ Z(\tau) = \frac{\tau}{4\pi \sinh(\tau)} = \frac{\tau}{2\pi(\exp \tau - \exp(-\tau))} . $$

A numerically more stable way of writing the vMF distribution in 3D is:

$$ \text{vMF}(n; \mu, \tau) = \frac{\tau \exp \left( \tau (\mu^T n - 1) \right)}{2\pi (1 - \exp(-2\tau))} . $$

Depiction of the von-Mises-Fisher (vMF) distributions on the unit sphere in 3D, $\mathbb{S}^2$, and 2D. The color encodes the probability density function value of the vMF over the whole sphere. As $\tau \to \infty$ the von-Mises-Fisher distribution approaches a delta function on the sphere at its mode $\mu$. From the coloring it can be observed that the von-Mises-Fisher distribution is isotropic. (Drag on the widget to move the mode direction $\mu$ and adjust the concentration $\tau$ using the slider).

Sampling in 3D To sample a random vector from a vMF distribution with mode $m = (0, 0, 1)$, first sample variables $v$ and $u$:

$$ v \sim \text{Unif}(\mathbb{S}^{1}), \qquad u \sim p(u; \tau) = \frac{\tau}{2 \sinh \tau} \exp(\tau u) $$

In practice, $v$ is obtained by sampling a zero-mean isotropic Gaussian, and $u$ is sampled via inverse transform sampling using a uniform random variable $\xi \sim \text{Unif}(0, 1)$:

$$ u = F^{-1}(\xi) = 1 + \tau^{-1} \log \left( \xi + (1 - \xi) \exp(-2\tau) \right) $$

The vMF-distributed direction $n$ in local coordinates is then:

$$ n = \begin{pmatrix} \sqrt{1 - u^2} v & u \end{pmatrix} $$

Finally, rotate $n$ from mode $m$ to mean direction $\mu$ using Rodrigues’ rotation formula ${}^\mu R_m = \text{Exp}\big([\theta\,\hat{w}]_\times\big)$ with unit axis $\hat{w} = \frac{m \times \mu}{\lVert m \times \mu \rVert}$ and angle $\theta = \arccos(\mu^T m)$.

Based on these observations, Loubet et al. introduced an additional convolution to handle integrals with large support as illustrated in Fig. (right).

To evaluate Equation $\eqref{eq:reparam_conv}$ numerically, we sample outer directions $\omega_i$ and offset directions $\mu_i \sim k(\cdot, \omega_i)$, giving the combined estimator:

The size of the spherical kernel $k$ provides a trade-off between variance and bias: narrower kernels model local edge displacements accurately but decrease the chance of finding discontinuities (increasing variance), while wider kernels reduce variance but increase rotation-approximation bias. For a fixed kernel, the practical estimator is not consistent: tracing more paths reduces variance but does not remove bias from an imperfect rotation or a missed discontinuity.

Determining Suitable Rotations Using Ray Tracing

Determining rotation matrices that track boundary motion without explicitly searching for silhouette edges is key to making this reparameterization practical for high scene complexity.

A central insight of Loubet et al. is that finding a suitable change of variables does not require identifying silhouette edges or even knowing whether an integrand contains a discontinuity. The only required information is how surface points move under infinitesimal perturbations of scene parameters $\pi$.

Because the integrand has small support, the displacement of points on silhouette edges closely approximates the displacement of nearby surface positions on the same object. We exploit this by tracing a small batch of auxiliary rays within the integrand’s support (where the ray count controls the trade-off between variance and the probability of missing a discontinuity). Using distance and surface normal information, a heuristic selects a candidate occluder point whose motion under parameter changes tracks that of the silhouette.

Projecting the selected point onto $S^2$ gives direction $\omega_P(\pi)$, with $\omega_{P_0} = \omega_P(\pi_0)$. A differentiable rotation matrix $R(\pi)$ is then constructed to satisfy:

Thus $R(\pi_0) = I$ (leaving primal ray tracing unchanged) while its derivative tracks the occluder’s motion to first order. Appendix B of Loubet et al. provides a closed-form formula for $R(\pi)$.

Differentiable Rotation Matrices

Evaluating local reparameterizations requires differentiable rotation matrices $R$ that map direction $\mathbf{\omega}_a$ to $\mathbf{\omega}_b$ such that $R\mathbf{\omega}_a = \mathbf{\omega}_b$. In particular, these matrices must be evaluated in the limit $\mathbf{\omega}_a \to \mathbf{\omega}_b$—where $R$ equals the identity matrix, but its gradients with respect to $\mathbf{\omega}_a$ and $\mathbf{\omega}_b$ remain non-trivial. Standard Rodrigues’ rotation [Belongie 2019] [16] suffers from a singularity at $\mathbf{\omega}_a = \mathbf{\omega}_b$ because it divides by the norm of the rotation axis. Loubet et al. factor out and simplify this norm, yielding the smooth expression:

$$R = \alpha I + [\mathbf{\beta}]_\times + \frac{1}{1 + \alpha} \mathbf{\beta} \mathbf{\beta}^T$$

where $\alpha = \mathbf{\omega}_a \cdot \mathbf{\omega}_b$, $\mathbf{\beta} = \mathbf{\omega}_a \times \mathbf{\omega}_b$ and

$$[\mathbf{u}]_\times = \begin{bmatrix} 0 & -u_z & u_y \\ u_z & 0 & -u_x \\ -u_y & u_x & 0 \end{bmatrix}$$

This formulation is well-defined and differentiable even when $\mathbf{\omega}_a = \mathbf{\omega}_b$.

Crucially, this construction requires only standard ray intersection queries (well suited for hardware acceleration) and the auxiliary rays can often be reused for Monte Carlo integration.

Variance Reduction Using Control Variates

A naive implementation of the change-of-variables estimator introduced above exhibits significant gradient variance. Loubet et al. resolve this by leveraging control variates constructed from correlated path pairs.

The Origin of Gradient Variance

Consider a general reparameterized Monte Carlo estimator evaluating an integral over domain $\mathcal{Y}$:

At $\pi = \pi_0$, the transformation is the identity $\mathcal{T}(y, \pi_0) = y$, giving $w_i(\pi_0) = 1$. The primal estimate $E = \frac{1}{N}\sum f(y_i, \pi_0)$ is unaffected.

Differentiating $E$ with respect to $\pi$ via the product rule yields:

Now consider a constant integrand $f(x) = c$. The true derivative of the integral $\int c \, \mathrm{d}x$ is zero. Consequently:

While the expected derivative of the weights $w_i(\pi)$ is zero for any distribution $k$, individual sample weight gradients $\frac{\partial w_i(\pi)}{\partial \pi}$ are non-zero and fluctuate randomly. These fluctuations introduce severe variance into gradient estimates.

Control Variates Formulation

The classical control variates method reduces the variance of an estimator $E$ using a correlated estimator $F$ with known expectation $\mathbb{E}[F]$:

where optimal variance reduction is achieved when $\alpha = -\frac{\operatorname{Cov}(E, F)}{\operatorname{Var}(F)}$.

Since $\mathbb{E}\left[ \frac{\partial w_i(\pi)}{\partial \pi} \right] = 0$, we construct the control variate $F(\pi) = \frac{1}{N}\sum_{i=1}^N w_i(\pi)$, whose gradient has zero expectation and whose value at $\pi_0$ is known exactly, $F(\pi_0) = 1$. This modifies the estimator to:

Differentiating $E'$ at $\pi = \pi_0$ gives:

If $f(x) = c$ is constant, setting $\alpha = -c$ makes $\frac{\partial E'}{\partial \pi} = 0$ for every sample point, completely eliminating gradient variance. For a general smooth function, $\alpha$ should therefore approximate the negative average value of $f$. It may reduce variance substantially without introducing bias, provided it is independent of the weight gradient to which it is applied.

Cross-Reduction with Correlated Path Pairs

To determine $\alpha$ without introducing bias (as $\alpha$ must remain independent of sample weights $w_i$), Loubet et al. employ a cross-weighting scheme using pairs of correlated paths ($r_0$ and $r_1$).

Path $i \in \{0, 1\}$ accumulates radiance across scattering bounces $l$:

where $W_{i,l}(\pi)$ is the product of reparameterization weights along path $i$ up to bounce $l$, and $f_{i,l}(\pi)$ represents throughput and emitter radiance.

By using $\alpha=-f_{1,l}$ for path 0 and $\alpha=-f_{0,l}$ for path 1, the cross-reduced path contribution $r'$ is:

Correlated path pairs reuse the random numbers for path-construction steps except those that affect the local reparameterization weights. Those samples remain independent so that $f_{0,l}$ is uncorrelated with $\partial_\pi W_{1,l}$ and vice versa. Under this independence condition, cross-reduction lowers gradient variance for direct and multi-bounce illumination without adding bias, at the cost of tracing paired paths.

Limitations

Variance and Bias of gradient estimates. The gradients have significant variance in several cases including high-order diffuse interreflections. Other cases exhibiting variance include rapidly changing radiance values, e.g., in very sharp shadows in locations where an occluder directly touches another surface that is being shadowed.

For a fixed convolution kernel size, the estimators are not consistent: they do not converge to the correct gradients when more rays are traced. While it is generally assumed that stochastic gradient descent requires unbiased gradients, accurate and low-variance gradients are preferable to unbiased but uninformative gradients for approaching a local minimum within few steps. However, bias could potentially deteriorate the convergence when approaching a solution during the last steps of an optimization.

Other light transport phenomena. The method relies on several simplifying assumptions that may not hold in particular scenarios. For example, the assumption that there is a single type of discontinuity within the support of an integrand may not hold when two discontinuities are very close to each other.

The method does not support perfectly specular materials and degenerate light sources containing Dirac delta functions.

Unbiased Warped-Area Sampling (Bangaru et al., 2020)

Bangaru et al. [8] ask whether the boundary term can be estimated using the same area samples as an ordinary path tracer. Their answer is yes: apply the divergence theorem to replace flux through visibility boundaries by divergence throughout the smooth interior. The resulting method does not enumerate or sample silhouette edges. This is different from merely smoothing visibility; the construction specifies conditions under which the area estimator represents the exact boundary derivative.

Boundary Integral in Differentiable Rendering

Let $D$ be an angular integration domain and partition it, only for the derivation, into disjoint regions $D_i(\boldsymbol{\pi})$ such that $f(\boldsymbol{\omega};\boldsymbol{\pi})$ is continuous within the boundaries of each piece, and all the discontinuities are at the boundaries of the domain. Therefore the domains $D_i$ are dependent on the scene parameters $\boldsymbol{\pi}$:

By applying Reynolds transport theorem $\eqref{eq:reynolds-transport-theorem}$, we get:

where $D_i'=D_i\setminus\partial D_i$. The first term is the usual interior derivative. The second measures the flux produced by moving discontinuities. Note that the jump $\Delta f$ of Equation $\eqref{eq:reynolds-transport-theorem}$ does not appear explicitly here: once the domain is partitioned, each piece contributes its own one-sided value of $f$ against its own outward normal, and the jump reappears when adjacent pieces are summed.

This partition is only a device used in the proof; evaluating the estimator does not require clipping the scene into the regions $D_i$ or enumerating their boundaries.

Area Form of the Boundary Derivative

The boundary term is still an integral over curves we would have to find. Bangaru et al.’s move is to introduce a vector field $\mathcal{V}_{\boldsymbol{\pi}}(\boldsymbol{\omega})$ that interpolates the boundary velocity into the interior, then apply the divergence theorem to $f\mathcal{V}_{\boldsymbol{\pi}}$ to convert that curve integral back into an integral over the interior.

Divergence theorem (Gauss-Ostrogradsky). Several vector calculus results (Green’s theorem, divergence theorem, Stoke’s theorem) relate the boundary integral to the interior integral:

Bangaru et al. apply the divergence theorem to rewrite the boundary derivative integral in the same domain as the interior derivative integral:

where the warp field $\mathcal{V}_{\boldsymbol{\pi}}(\boldsymbol{\omega})$ is a smooth interpolation of the boundary velocity $\partial_{\boldsymbol{\pi}}\boldsymbol{\omega}$. The last step expands the divergence product rule to separate the advection and domain-change terms for the Monte Carlo estimation set up later.

Consequently, an area sample contributes three conceptually different derivatives:

The first is the ordinary interior derivative. The second moves the sample with the warp. The third accounts for local expansion or contraction of the warped domain. Dropping the divergence term is only valid for volume-preserving warps.

Validity Conditions for the Warp

There are infinitely many warp fields that satisfy the divergence equivalence. To ensure the applicability of the divergence theorem, a warp field $\mathcal{V}_{\boldsymbol{\pi}}(\boldsymbol{\omega})$ is defined as valid for a boundary velocity $\partial_{\boldsymbol{\pi}}\boldsymbol{\omega}^{(b)}$ if and only if it satisfies two strict criteria:

1. Interior Continuity ($C^0$ Continuity): The field must be continuous everywhere inside the open interior domain $D' = D \setminus \partial D$:

   $$\mathcal{V}_{\boldsymbol{\pi}}(\boldsymbol{\omega}) \in C^0(D'), \quad \forall \boldsymbol{\omega} \in D'$$

2. Boundary Consistency ($\epsilon$-$\delta$ Limit): As an interior direction $\boldsymbol{\omega}$ approaches any boundary point $\boldsymbol{\omega}_b \in \partial D$, the warp field must converge to the true boundary velocity $\partial_{\boldsymbol{\pi}}\boldsymbol{\omega}_b$:

   $$ \begin{aligned} &\forall \boldsymbol{\omega}_b \in \partial D \quad \text{and} \quad \forall \delta \in \mathbb{R}^+, \quad \exists \epsilon > 0 \quad \text{such that} \\ &\forall \boldsymbol{\omega} \in D' \text{ with } \|\boldsymbol{\omega} - \boldsymbol{\omega}_b\| < \epsilon \implies \|\mathcal{V}_{\boldsymbol{\pi}}(\boldsymbol{\omega}) - \partial_{\boldsymbol{\pi}}\boldsymbol{\omega}_b\| < \delta \end{aligned} $$

   Or equivalently in limit notation:

   $$\lim_{\boldsymbol{\omega} \to \boldsymbol{\omega}_b} \mathcal{V}_{\boldsymbol{\pi}}(\boldsymbol{\omega}) = \partial_{\boldsymbol{\pi}}\boldsymbol{\omega}_b, \quad \forall \boldsymbol{\omega}_b \in \partial D$$

The continuity condition guarantees that no spurious interior jump discontinuities are introduced during differentiation. The boundary consistency condition enforces that the field smoothly interpolates the exact boundary motion near silhouettes. Satisfying both conditions guarantees that inserting $\mathcal{V}_{\boldsymbol{\pi}}(\boldsymbol{\omega})$ into the area integral yields an exact, unbiased representation of the boundary derivative.

The equality is required as a limit from the smooth regions; the field need not be defined on the measure-zero boundary itself. These conditions are the central correctness criterion of the paper. A smooth field with the wrong boundary value remains biased, and a field that is correct at a silhouette but discontinuous in the interior cannot be inserted into the area formula above.

Valid warp fields for the rendering integral

The smooth boundary interpolation problem of constructing the warp field $\mathcal{V}_{\boldsymbol{\pi}}(\mathbf{x})$ presents a unique challenge in rendering. In many other fields such as geometry processing or numerical computation, it is often possible to identify all the boundaries a priori, and then discretize the interior in order to construct the smooth field. However, in the area-based rendering formulation, the goal is to avoid explicit discontinuity enumeration (like explicit boundary sampling) in the first place.

More concretely, this is a blind interpolation problem: the method must construct a valid warp field that can be computed without explicit samples on the boundary.

Exploiting structure in the scene. Bangaru et al.’s approach exploits the manifold structure of the 3D scene, by observing that the derivative of a 3D scene point $\partial_{\boldsymbol{\pi}} \mathbf{x}$ is continuous for all surface points $\mathbf{x} \in \mathcal{X}$. The discontinuities in the visibility term only arise when the geometry is projected to the solid angle space $\mathcal{X} \to \Omega$ of the point where the radiance is being evaluated.

Selecting the warp field. The goal is to define a field $\mathcal{V}_{\boldsymbol{\pi}}(\boldsymbol{\omega})$ for all $\boldsymbol{\omega} \in \Omega$ such that the continuity and boundary consistency conditions are satisfied. Bangaru et al.’s strategy is to construct a field that satisfies the boundary consistency, but is not necessarily continuous, by differentiating the ray-geometry intersection procedure.

The rendering integral maps a solid angle $\boldsymbol{\omega}$ at position $\mathbf{x}$ to a scene point $\mathbf{y}$ through the ray-scene intersection operator, denoted $\mathbf{y} = \text{INTERSECT}(\mathbf{x}, \boldsymbol{\omega}; \boldsymbol{\pi})$. Specifically, the derivatives of the intersection function with respect to the scene parameters $\boldsymbol{\pi}$ serve as the initial (invalid) warp field, obtained by automatically differentiating the intersection function to get $\mathbf{y}, \partial_{\boldsymbol{\pi}} \mathbf{y}, \partial_{\boldsymbol{\omega}} \mathbf{y} = \text{DIFF-INTERSECT}(\mathbf{x}, \boldsymbol{\omega}; \boldsymbol{\pi})$. Concretely, the direct warp field is:

The division by the Jacobian $\vert{}\partial_{\boldsymbol{\omega}} \mathbf{y}\vert{}$ converts the measure of $\partial_{\boldsymbol{\pi}} \mathbf{y}$ from area measure to solid angle measure. This projection term is the same as the geometry term for converting between solid angle and area formulations in path-space rendering methods. It projects the derivative instead of the radiance.

The warp field $\mathcal{V}_{\boldsymbol{\pi}}^{(\text{direct})}(\boldsymbol{\omega})$ satisfies the boundary consistency criterion, since at the points close to the boundary, the derivative $\frac{\partial_{\boldsymbol{\pi}} \mathbf{y}}{\vert{}\partial_{\boldsymbol{\omega}} \mathbf{y}\vert{}}$ approaches the boundary derivative $\partial_{\boldsymbol{\pi}} \boldsymbol{\omega}_b$.

Intuitively, this states that the rate at which a given point $\boldsymbol{\omega}$ moves is equal to the motion of the corresponding 3D point adjusted by the Jacobian of the projection between the spaces (Fig. ). Unfortunately, this warp field is not valid since it breaks the continuity criterion. For example, consider two angles close together on either side of a boundary, but which intersect different surfaces, and therefore have very different warps.

The discontinuity of the warp field $\mathcal{V}_{\boldsymbol{\pi}}^{(\text{direct})}(\boldsymbol{\omega})$ is also why it is difficult to differentiate the visibility function. In general, for a discontinuity in the visibility function, only one side is representative of the edge. Evaluating the warp at a point on the other side of the discontinuity yields the motion of a completely independent 3D point. Still, in order to create a continuous warp, the field must somehow be made continuous.

Bangaru et al.’s solution is to convolve the direct warp field with weights that ensure boundary consistency. This rectifies the underlying discontinuous field $\mathcal{V}_{\boldsymbol{\pi}}^{(\text{direct})}(\boldsymbol{\omega}')$ to produce a smooth and continuous field $\mathcal{V}_{\boldsymbol{\pi}}^{(\text{filtered})}(\boldsymbol{\omega}; w)$:

Next, the question is how to choose the weights $w$ such that boundary consistency is preserved.

Boundary-aware convolution. It is not immediately obvious what the weights should be. As a counter example, a warp field obtained using weights from a normal distribution deviates heavily from the true warp at the boundary, and this field is still not valid since it violates boundary consistency. The warp at a point very close to the boundary would be the average of the warp on both sides, which is not, in general, equal to (or even close to) the warp at the boundary.

The weights need to converge to the derivative at points close to the boundary to produce a valid interpolation. That is, $w(\boldsymbol{\omega}, \boldsymbol{\omega}')$ should grow to infinity when $\boldsymbol{\omega}$ is on the boundary while $\boldsymbol{\omega}'$ approaches $\boldsymbol{\omega}$. The weight should also be small when $\boldsymbol{\omega}$ is far from the boundary. Drawing inspiration from harmonic interpolation, Bangaru et al. select weights using the inverse distance to the boundaries.

Unfortunately, the true distance to the nearest boundary point is unknown because finding the boundary is difficult. Instead, the paper introduces the concept of a boundary-test $\mathcal{B}(\boldsymbol{\omega}')$, which serves as a weaker definition of a boundary. $\mathcal{B}(\boldsymbol{\omega}')$ is essentially a soft indicator function which takes zero value at the boundary and non-negative everywhere else, i.e., $\boldsymbol{\omega}' \in \partial\Omega \implies \mathcal{B}(\boldsymbol{\omega}') = 0$. This is the only key requirement that $\mathcal{B}(\boldsymbol{\omega}')$ needs to satisfy, which gives great flexibility over its form.

The harmonic weights then take the form:

where $\mathcal{D}(\boldsymbol{\omega}, \boldsymbol{\omega}')$ is the distance between $\boldsymbol{\omega}$ and $\boldsymbol{\omega}'$. The paper uses the von-Mises Fisher distance $\mathcal{D}(\boldsymbol{\omega}, \boldsymbol{\omega}') = e^{\kappa (1 - \langle\boldsymbol{\omega}, \boldsymbol{\omega}'\rangle)} - 1$, but any smooth function that satisfies $\mathcal{D}(\boldsymbol{\omega}, \boldsymbol{\omega}) = 0$ can be used.

These weights satisfy the handy asymptotic property that in the limit where a point $\boldsymbol{\omega}^{(b)}$ approaches a boundary (denoted by $\partial\Omega$),

where $\delta(\cdot)$ is the Dirac delta function. This implies that for a point $\boldsymbol{\omega}^{(b)}$ very close to the boundary, the weights are such that the resulting warp at $\boldsymbol{\omega}^{(b)}$ is equal to the direct warp. As shown above, the direct warp is consistent at the boundaries, which implies that the convolution with harmonic weights is also consistent.

The analysis uses a limit near the boundary and not the boundary itself, since the warp field only needs to be defined in the smooth interior region $\Omega - \partial\Omega$. This allows the warp field to be undefined at points on the boundary, and the proposed harmonic weights are infinite at such points.

Pixel prefiltering. While, in principle, this method could handle discontinuous pixel prefiltering such as a box filter, Bangaru et al. always opt for a continuous filter. For a discontinuous filter, the pixel filter integral $\int_A f(\boldsymbol{\omega})\mathrm{d}\boldsymbol{\omega}$ over the pixel support $A$ needs to take the discontinuity at the support’s boundary into account for the divergence theorem to hold. This means the boundary test $\mathcal{B}(\boldsymbol{\omega}')$ needs to return $0$ at such boundaries. To simplify the boundary test implementation, the paper always uses a Gaussian pixel filter truncated at a radius of 4 times the standard deviation, where the kernel contribution is indistinguishable from the floating point error.

Relation to the Reparameterization Method

Before moving to the Monte Carlo sampling algorithm for solving the harmonic convolution integral, Bangaru et al. discuss the relation of their method to Loubet et al.’s reparameterization.

Although the area form of the boundary term is derived in a different way, Loubet et al.’s method of transforming the samples using a rotation is actually a special case of this formulation. Bangaru et al. show that it can be interpreted as applying a particular warp field that does not satisfy the boundary consistency requirement, thus introducing bias to the result.

The reparameterization method applies a transformation $\mathbf{x} = \mathcal{T}(\mathbf{y}; \boldsymbol{\pi})$ in an attempt to remove the discontinuities of the integral:

Bangaru et al. prove that this transformation can be converted to an equivalent warp field $\mathcal{V}_{\boldsymbol{\pi}}(\mathbf{x})$ using the relationship,

where $\boldsymbol{\pi}_0$ is the point at which the derivative is to be computed.

Given a warp field, it can also be converted back into a transformation by finding a solution for the right-hand side of this differential equation. There exist many possibilities for $\mathcal{T}(\mathbf{x}; \boldsymbol{\pi})$’s basic form, depending on the coordinate system (e.g., linear solutions for 2D Euclidean or rotational solutions for 3D unit-sphere coordinate systems). The rotational solution represents the basic transformation that Loubet et al. uses.

Denoting the rotational solutions by $\mathcal{R}(\boldsymbol{\omega}; \boldsymbol{\pi})$, this limits the set of possible transformations to those with a unit Jacobian $\vert{}\partial_{\boldsymbol{\omega}} \mathcal{R}(\boldsymbol{\omega}; \boldsymbol{\pi})\vert{} = 1$, independent of $\boldsymbol{\pi}$. It follows that the corresponding warp satisfies $\nabla_{\boldsymbol{\omega}} \cdot \mathcal{V}_{\boldsymbol{\pi}}(\boldsymbol{\omega}) = 0$. Unfortunately, this means that simple rotation-based reparameterization cannot always be a valid warp, as some scenarios explicitly require a non-zero divergence.

Loubet et al. recognize this drawback and propose a more complex transformation that samples nearby rays using a von-Mises Fisher distribution (spherical analog of the Gaussian) and constructs a transformation that is the average of these rotations. This can generally be expressed using $\mathcal{R}'(\boldsymbol{\omega}; \boldsymbol{\pi}) = \int k(\boldsymbol{\omega}, \boldsymbol{\omega}') \mathcal{R}(\boldsymbol{\omega}'; \boldsymbol{\pi})$. Since the transformation and the warp field are linearly related, the warp field of a convolved transformation $\mathcal{V}_{\boldsymbol{\pi}}'(\mathbf{x}) = [\partial_{\boldsymbol{\pi}} \mathcal{R}'(\boldsymbol{\omega}; \boldsymbol{\pi})]_{\boldsymbol{\pi}=\boldsymbol{\pi}_0}$ is equivalent to the convolution over the warp field of the original transformation $\mathcal{V}_{\boldsymbol{\pi}}'(\mathbf{x}) = \int k(\boldsymbol{\omega}, \boldsymbol{\omega}') \mathcal{V}_{\boldsymbol{\pi}}(\mathbf{x})$. This leads to a scenario where the resulting transformation is smooth but fails to meet the boundary consistency criterion.

To compensate for this bias, Loubet et al. introduce a heuristic on top of this convolution. Since the heuristic involves discrete operations such as sorting and comparing object IDs, it is difficult to analytically express the resulting warp field and study its properties.

Monte Carlo Estimation of the Derivative

With the theory of area sampling and the formulation for the convolutional warp field established, the next step is to develop a Monte Carlo estimator for the divergence area integral.

Because the convolutional warp field is itself an integral, a nested Monte Carlo estimator is required. The process works as follows: a primary sample $\boldsymbol{\omega}$ is first generated for the outer divergence integral, and then a set of auxiliary samples $\{\boldsymbol{\omega}_1' \cdots \boldsymbol{\omega}_N'\}$ is generated to estimate the inner convolution integral precisely at $\boldsymbol{\omega}$.

A major mathematical hurdle arises here: the convolution integral of the warp field contains a division for normalization. A naïve Monte Carlo estimator of the reciprocal of an integral is inherently biased, because the expected value of a reciprocal is not equal to the reciprocal of the expected value ($\mathbb{E}[1/f] \neq 1/\mathbb{E}[f]$).

Fortunately, an unbiased Monte Carlo estimator can be constructed from a consistent one using Russian Roulette de-biasing.

Estimating the Warp Field $\mathcal{V}_{\boldsymbol{\pi}}(\cdot)$

The goal is to estimate the divergence area integral, whose integrand is:

For each direction sample $\boldsymbol{\omega}$ for the outer divergence integral, $N'$ auxiliary directions $\boldsymbol{\omega}_i'$ are sampled to estimate the warp field $\mathcal{V}_{\boldsymbol{\pi}}(\boldsymbol{\omega})$. The warp field relies on a kernel of harmonic weights. While the distribution of these weights depends on the configuration of silhouette edges, it is heavily correlated with a normal distribution centered exactly at the outer directional sample $\boldsymbol{\omega}$.

Therefore, Bangaru et al. importance-sample from a normal distribution around the outer directional sample. For spherical or hemispherical sampling, the von Mises-Fisher (vMF) distribution is used.

If we use a finite, fixed number of auxiliary rays $N'$, the estimator is consistent (meaning it converges to the true derivative as $N' \to \infty$) but strictly biased due to the division by the integral in the harmonic convolution. In practice, a fixed auxiliary ray count between 4 and 64 provides a highly robust, albeit biased, estimate.

De-biasing the Estimator

To completely eliminate bias, we utilize a Russian Roulette technique introduced by McLeish, commonly used in Bayesian inference and physically-based rendering.

The intuition is to take a converging sequence $T_0, T_1, T_2, \dots$ produced by a consistent estimator and construct an infinite series:

Since $\hat{T}$ algebraically collapses to $T_\infty$ (the right answer), building an unbiased estimator for this infinite series turns our consistent estimator into a perfectly unbiased one.

Step-by-Step Derivation: McLeish Debiasing

1. Let $T_i$ be our biased estimate of $1/\mathbb{E}[f]$ using $i$ samples.
2. Let $\Delta_i = T_i - T_{i-1}$ be the difference between successive estimates.
3. The true value is $T_\infty = T_0 + \sum_{i=1}^{\infty} \Delta_i$.
4. To estimate this infinite sum without infinite compute time, we introduce a random integer $N'$ drawn from a distribution with probability mass function $\mathbb{P}(N'=n) = p_n$, and cumulative survival probability $\mathbb{P}(N' \ge i) = P_i$.
5. We construct the randomized estimator: $\hat{T} = T_0 + \sum_{i=1}^{N'} \frac{\Delta_i}{P_i}$.
6. Taking the expected value of $\hat{T}$, the probabilities cancel out: $\mathbb{E}[\hat{T}] = T_0 + \sum_{i=1}^{\infty} P_i \left( \frac{\Delta_i}{P_i} \right) = T_0 + \sum_{i=1}^{\infty} \Delta_i = T_\infty$.
7. Thus, $\hat{T}$ is entirely unbiased.

To implement this, $N'$ is treated as a random variable following a geometric distribution. For every path bounce, during the auxiliary sampling step, a new $N'$ is sampled for the estimation of the warp field, which makes the derivative estimate unbiased.

Below are the algorithms detailing both the standard consistent estimator and the fully unbiased Russian Roulette formulation.

(Note: $\text{Geom}(i; p)$ in Algorithm 3 denotes the survival probability $\mathbb{P}(N' \ge i)$ used in the derivation above, not the geometric probability mass function. Also note that in practice, using the consistent version with a high $N'$ to strictly reduce bias is often faster than Russian Roulette, since managing a dynamic number of rays per pixel can cause warp divergence and performance bottlenecks on memory-constrained GPUs).

Variance Reduction

Even when unbiased, some parts of this estimator exhibit exceptionally high variance if used without explicit variance reduction. Surprisingly, smooth regions suffer significantly due to the structural variation of both the harmonic weights and the warp field divergence.

To understand why, consider an infinite, flat emitter facing the camera directly, where the scene parameter $\pi$ is its spatial translation. Since it faces the camera uniformly, the expected gradient of the image is exactly 0. However, the spatial weight derivatives $\nabla_{\boldsymbol{\omega}} w(\boldsymbol{\omega}, \boldsymbol{\omega}')$ evaluated at individual samples can be extremely large depending on their location relative to the distribution center. Thus, the estimator bounces around wildly with high variance, even in a scene with zero geometric complexity and a true gradient of 0.

Antithetic Variates. To counter this, Bangaru et al. apply antithetic variates, pairing each auxiliary sample with its 180°-rotated counterpart about the distribution center. The resulting pair is negatively correlated, so when averaged, the symmetric noise components of the harmonic weight derivatives cancel each other, dramatically collapsing variance in smooth regions.

Control Variates. When the underlying geometry is at an angle (so the warp field varies approximately linearly rather than uniformly), antithetic variates alone are insufficient. In this regime, the intersection derivative $\partial_\pi \boldsymbol{\omega}$ varies linearly, and the divergence of a locally planar surface can be computed analytically. Bangaru et al. use this closed-form linear approximation as a control variate, subtracting its known-mean residual from the Monte Carlo estimator to absorb the dominant source of noise without introducing bias. Together, antithetic and control variates are indispensable for the low sample counts typical of iterative inverse rendering.

How Little the Boundary Test Has to Do

Before turning to limitations, one structural observation is worth extracting, because it explains why such a loose boundary test suffices. Writing $Z(\boldsymbol{\omega})=\int w(\boldsymbol{\omega},\boldsymbol{\omega}')\,\mathrm d\boldsymbol{\omega}'$ for the normalisation of the convolution, the quotient rule gives

Only the derivative of $w$ with respect to the primary direction $\boldsymbol{\omega}$ appears. Therefore $\mathcal{B}(\boldsymbol{\omega}')$ itself need not be differentiable, or even continuous, across visibility boundaries: a differentiable distance function $\mathcal{D}(\boldsymbol{\omega},\boldsymbol{\omega}')$ alone is enough to keep the divergence finite at interior points. This is what makes the boundary test cheap to design for a new geometric representation.

Generalization and Limitations

With that in place, the paper is explicit about where its unbiasedness guarantee does and does not apply:

• Implicit edges. The triangle-mesh boundary test $\mathcal{B}$ relies on face normals and open/silhouette edges. It does not automatically detect implicit edges created by triangle self-intersections, where geometric boundaries do not coincide with mesh edges. Near such intersections, $\mathcal{B}$ fails to vanish, so the convolution averages neighbouring warps rather than applying singular weights. The result remains visually accurate and behaves much like Loubet et al., but the field loses its strict boundary-consistency guarantee and the method is no longer provably unbiased.
• Unbounded work. The Russian-roulette estimator has theoretically unbounded work and storage, so in practice it must be truncated (the paper caps $N' \leq 512$). The resulting truncation bias is normally indistinguishable from floating-point error, but can become non-negligible with very small kernel sizes.
• Incoherent workload. Because each path vertex may draw a different number of auxiliary rays, the debiased variant produces thread divergence and dynamic allocation that complicate GPU parallelism. This is a large part of why the fixed-$N'$ consistent estimator is preferred in practice.
• Domain extensions. The divergence argument extends in principle to motion blur and depth of field by enlarging the integration domain and redefining $\mathcal{B}$. Extension to path space is less immediate, since a single path-space point can cross several occluders and the correct boundary definition in that setting remains open.
• Universal boundary test. The method augments a unidirectional path tracer and naturally supports secondary transport, but the paper does not claim its particular triangle boundary test is universal. Other representations (SDFs, Bézier curves) require their own test satisfying the same limiting condition.

Projective Sampling for Differentiable Rendering of Geometry (Zhang et al., 2023)

Boundary-sampling methods like Li et al.’s edge sampling place additional Monte Carlo samples directly on silhouettes, while the reparameterization methods above trade that explicit search for a warp field that absorbs the boundary term into the interior integral. Zhang, Roussel, and Jakob (2023) [5] propose a hybrid of the two: instead of building a dedicated sampler for the boundary, or reparameterizing the whole domain, they repurpose the ordinary primal samples that a path tracer already generates (BSDF, emitter, and MIS samples) by projecting each resulting ray segment onto a nearby silhouette. Whatever primal sampling strategy is good for the image is thereby automatically put to work on its derivative too.

A Local Boundary Integral

To place Projective Sampling in context, consider the path-space differentiable rendering approach of Zhang et al. (2020) [9], which explicitly samples points on silhouette edges and connects them to full light paths to handle secondary visibility. Finding these relevant edges in a complex 3D scene is computationally demanding, often requiring the construction of heavy auxiliary data structures (such as spatio-angular photon maps) prior to rendering to guide samples toward important boundaries. (For the complete derivation of the path-space framework, please see Zhang et al. (2020) [9]).

The theoretical starting point for Zhang et al. (2023) is this same path-space formulation, which decouples the effect of boundaries from their interior. Expressed with respect to a path segment $(\mathbf{x}_a, \mathbf{x}_c)$, the pixel derivative with respect to a scene parameter $\boldsymbol{\pi}$ states:

The integral is over segments $(\mathbf{x}_a, \mathbf{x}_c)$ that make contact with surface boundary along the way. The domain $\mathcal{B}(\mathbf{x}_a)$ describes the “shadow” of this boundary and is further composed of $\mathcal{B}(\mathbf{x}_a) = \cup_i \mathcal{B}_i(\mathbf{x}_a)$ (one set $\mathcal{B}_i$ per edge) when the scene consists of discrete geometry. $L_i$ and $W_i$ refer to incident radiance and importance, and $G$ is the standard geometric term. The inner product measures the perpendicular speed of the “shadow” at $\mathbf{x}_c$.

In earlier path-space formulations, the final boundary integral expressions retained implicit derivatives through ray-tracing operations evaluated via automatic differentiation, which obscured potential analytical simplifications.

Zhang et al. (2023) re-derive the boundary integral into an explicit local formulation for both perimeter and interior components, yielding the following reduced expression:

The term $L_d$ stands for the radiance difference between foreground and background $L_d(\mathbf{x}_b, \boldsymbol{\omega}) = L_o(\mathbf{x}_b, \boldsymbol{\omega}) - L_i(\mathbf{x}_b, -\boldsymbol{\omega})$. In the perimeter term (top integral), $\phi$ is the angle between $\boldsymbol{\omega}$ and the boundary tangent $\mathbf{t}_b$. In the interior term, the angle $\phi \in \mathcal{S}^1$ parameterizes all relevant quantities over tangential directions at the surface position $\mathbf{x}_b$, and $\kappa(\phi)$ denotes the normal curvature.

The following aspects are noteworthy:

1. Locality: Neither term involves the complex non-local domain $\mathcal{B}(\mathbf{x}_a)$.
2. Absence of the Geometric Term: Contrasting with \eqref{eq:zhang2020_boundary}, the geometric term $G$ is absent in both integrals, which means that variance arising from this factor can be avoided in Monte Carlo methods.
3. Shape Specialization: When the scene geometry is smooth and closed, $\partial\mathcal{A} = \emptyset$ removes the first term. Polygonal meshes do not have curved interior ($\kappa=0$) and hence do not require the second term.
4. Resemblance to Rendering Equation: The sine (perimeter) and curvature (interior) terms resemble the cosine factor in the rendering equation. The influence of an edge with tangent $\mathbf{t}_b$ exerted in direction $\boldsymbol{\omega}$ tends to zero as $\mathbf{t}_b \to \boldsymbol{\omega}$, since the edge becomes invisible. Likewise, the curvature term accounts for foreshortening in the mapping between silhouette positions and scene positions observing this silhouette.

Sampling Strategies and Guiding

Because the method consumes ordinary primal samples, the first question is which primal distributions to feed it. Zhang et al. find that the established intuition transfers directly: projected BSDF sampling fails to produce enough silhouette samples in rough reflections, projected emitter sampling has the mirror-image problem in smooth reflections, and combining both through multiple importance sampling behaves exactly as it does in the primal setting.

The projection also does not have to be exact. Since the projected samples only guide a subsequent integration phase that accounts for all surface positions simultaneously, it suffices to find an approximate boundary segment $\mathbf{x}'_a \to \mathbf{x}'_b$ near the original $\mathbf{x}_a \to \mathbf{x}_b$, allowing the origin to shift slightly as well. This buys two practical freedoms: when a projection fails outright, the tentative endpoint can be snapped to the nearest local boundary and paired with the tangential direction closest to the original; and when a projection needs root-finding, the iteration can stop after a few steps rather than converging to machine precision. The first reduces variance, the second reduces the cost of building the guiding distribution.

Given a set of projected samples, the next step is condensing them into a structure that supports efficient sampling and density evaluation. For polygonal meshes the boundary sample space is three-dimensional: a single parameter $t \in [0,1]$ enumerating the mesh edges $\partial\mathcal{A}$, plus a direction $(\theta, \phi)$.

Grid-based guiding. The simplest option is a dense 3D density grid over $(t, \theta, \phi)$, which is also what Zhang et al.’s earlier path-space method relies on. Grids have one attractive property: the projection can accumulate density straight into voxels without ever storing the samples, which allows high-quality statistics in memory-constrained settings. The cost is that a high-resolution grid is memory-hungry and becomes the bottleneck in complex scenes that need fine discretization to resolve sparse features.

Hierarchical guiding. The integrand on boundary sample space is extremely sparse, and its sharp features grow more pronounced with scene complexity: narrow BSDF peaks, strongly peaked emitters, second-order visibility where silhouettes are themselves occluded, and plain discontinuities in the edge parameterization. Zhang et al. therefore replace the grid with a set of octrees, one per equal-sized interval of the $t$ axis. The critical difference from prior adaptive schemes is the construction order. Yan et al. (2022) build kd-trees top-down, subdividing until a local smoothness criterion is met, which inherits the failure mode of adaptive quadrature: a criterion based on local evaluations will sometimes miss a sharp peak and stop early. Zhang et al. instead build bottom-up, starting from the projected samples, which already concentrate at exactly those sparse features, and subdividing each node until it reaches a maximum depth or holds at most one sample. The projection is thus used not to accumulate density but to cheaply construct a high-fidelity adaptive partition; the actual integrand value is then estimated by drawing a fixed number of uniform samples inside each leaf, which requires no further projections.

Projection

Supporting a new geometric representation requires two operations: a projection mapping a ray segment $\mathbf{x}_a \to \mathbf{x}_b$ onto a nearby segment $\mathbf{x}'_a \to \mathbf{x}'_b$ that is tangential at $\mathbf{x}'_b$, and a parameterization of the perimeter (where one exists) and of the interior (for curved geometry). Spheres are the trivial case and make a good implementation test, since being smooth and closed they need only a surface mapping. Meshes are more involved.

Triangle meshes

Two projection strategies are used, both relying only on local differential information.

JUMP is a Newton-style iteration on a local linear model. The neighbourhood of a surface position $\mathbf{p}$ is approximated by $\tilde{\mathbf{p}}(u, v) = \mathbf{p} + u\,\partial_u \mathbf{p} + v\,\partial_v \mathbf{p}$ with an interpolated normal $\tilde{\mathbf{N}}(u, v) = \mathbf{N} + u\,\partial_u \mathbf{N} + v\,\partial_v \mathbf{N}$. Holding the viewing direction fixed, the silhouette of this approximation is a line in $(u, v)$ parameter space, so the solution $(u', v')$ follows analytically. A perpendicular ray is then traced from the inferred silhouette point to find the next intersection, from which the procedure could be repeated.

WALK is a greedy search across the mesh. It visits the three neighbouring triangles, computes the angle between each of their normals and the viewing direction, and steps toward increasing angle until reaching a perpendicular ($90^\circ$) or back-facing triangle. Consistently walking to the neighbour with the largest angle, however, attracts too many projections onto a handful of mesh edges, so the two largest-angle neighbours are instead used as discrete probabilities and one is picked at random.

The two are complementary. JUMP is aggressive and can bypass plateaus and heavily tessellated regions in one step, but requires a ray trace. WALK takes smaller steps and robustly detects nearby silhouettes even on bumpy geometry where the extrapolated local model is deceptive, and needs no ray tracing at all, making individual steps much cheaper. The hybrid used in practice therefore front-loads WALK and reaches for JUMP only when it stalls:

def find_silhouette(ray_origin, initial_triangle):
    # Walk for 30 steps
    current_triangle = initial_triangle
    for _ in range(30):
        current_triangle = WALK(ray_origin, current_triangle)
        if is_silhouette(current_triangle): return current_triangle

    # Jump once if no silhouette was found
    current_triangle = JUMP(ray_origin, current_triangle)
    if is_silhouette(current_triangle): return current_triangle

    # Followed by an additional 30 Walk steps
    for _ in range(30):
        current_triangle = WALK(ray_origin, current_triangle)
        if is_silhouette(current_triangle): return current_triangle

    return None

The whole procedure costs a single ray-tracing step. If it still fails to find a silhouette $\mathbf{x}'_b$, the tentative position is kept and paired with the tangential segment closest to the original direction, per the relaxation described above.

The local formulation of the boundary derivative is what makes the remaining two representations tractable at all, since both derive gradients from curved interiors rather than discrete edges.

Fiber curves

Fibers are modelled with a parametric base curve $\mathbf{C}(v)$ and radius $r(v)$, both cubic B-spline interpolants. Fixing $v$ gives a circular cross-section with centre $\mathbf{C}(v)$, radius $r(v)$, and normal $\mathbf{C}'(v)$; assigning an azimuth parameter $u$ to that circle yields a $C^1$-continuous surface $\mathbf{M}(u, v)$. Setting the curve endpoints aside (they can be capped with spheres), only the curved-interior half of the local boundary integral is needed.

Given a viewpoint $\mathbf{O}$ and a surface position $\mathbf{P} = \mathbf{M}(u_0, v_0)$, the projection holds $v_0$ fixed and solves for the $u$ that lands on a silhouette, $\langle \mathbf{M}(u, v_0) - \mathbf{O}, \mathbf{n}(u, v_0) \rangle = 0$. This expands into

with $A, B, C, D, E$ analytically computable and independent of $u$. A constant fiber radius admits a closed-form solution; otherwise 20 bisection iterations are fast and robust enough.

Implicit functions

The final representation is a signed distance function stored as a grid-based trilinear interpolant, chosen for two reasons. Intersections are fast, since they can be found using a mesh proxy plus analytic solutions inside each voxel rather than sphere tracing; and the low-order representation keeps the mapping from $\mathbb{R}^2$ to the surface simple. (Note that the signed-distance property is never actually used, so much of this generalises to other implicit surfaces.)

The trilinear scheme has one clear cost: geometric normals jump across voxel boundaries, so both the interior and perimeter terms of the local boundary integral are live. Parameterizing them relies on two observations. For the interior, trilinear interpolation guarantees that a given $x, y$ coordinate inside a voxel has at most one root along $z$. For the perimeter, the surface’s intersection with a voxel face is uniquely determined by the face index and a perpendicular coordinate, which flattens into a 1D mapping over all perimeter curves much like the mesh case. Together these give a per-voxel mapping that is globally discontinuous but locally well-behaved. Since the mapping depends on which dimension parameterizes the curve or surface, and becomes unstable when that dimension is close to perpendicular, the most numerically stable axis is assigned per voxel from the SDF gradient.

Limitations

The method’s main drawbacks all stem from the central projection step: the conversion of interior to boundary samples only works well when there is sufficient surface area to collect samples, and when the projection itself is well-behaved.

A thin shape such as a blade of grass, for instance, has a large ratio of boundary arclength to surface area and may never receive enough samples; this remains a fundamental limitation of the approach. Separately, a function that maps every interior point onto a single fixed silhouette location is technically a projection, but not a useful one. The authors answer this second concern by exhibiting high-quality projections for the representations above, with one gap worth noting: no projection operator was implemented for SDFs, so those results fall back on a uniformly initialised grid.

Comparing the Four Boundary Estimators

All four methods estimate the same object, the boundary term of the Reynolds transport theorem. What separates them is how they locate it, and what they pay for doing so.

The trend across the four is a steady retreat from explicit geometric search: from enumerating silhouette edges, to never looking for them, to letting the primal sampler find them as a side effect.

Physics-Based Differentiable Rendering Theory

Everything so far has concerned a single integral and its moving boundary. Rendering, however, is a recursive integral, and gradients must propagate through every bounce, not just the one being differentiated. We now differentiate the rendering equation itself to obtain that recursion in closed form. The result, the differential rendering equation, is what the efficient reverse-mode algorithms of the next section actually solve; the boundary estimators above supply one of its source terms.

Differentiable Rendering of Surfaces

Physically-based rendering of surfaces has been a central topic in computer graphics for decades and is governed by the well-known rendering equation (RE). The RE is an integral equation stating that the (steady-state) outgoing radiance $L_o$ at any surface point $\mathbf{x}$ with direction $\boldsymbol{\omega}_o$ is given by:

where $L_i$ is the incident radiance, $\mathrm{d}\sigma$ is the solid-angle measure, and $f_s$ denotes the BSDF multiplied by the cosine factor $|\mathbf{n}\cdot\boldsymbol{\omega}_i|$. We use a BSDF rather than a reflection-only BRDF because general scenes can include both reflection and transmission, for example at refractive interfaces. Light directions are written as bold unit vectors (e.g., $\boldsymbol{\omega}_i, \boldsymbol{\omega}_o$).

The RE has no analytical solution in general, and numerous numerical methods have been developed. Some of the widely adopted examples include unbiased methods like unidirectional and bidirectional path tracing, as well as biased ones such as photon mapping and lightcuts.

Direct Illumination

Before differentiating the full RE $\eqref{eq:rendering-equation}$, we will first consider the case of direct illumination as a warm-up. Specifically, the radiance $L_r$ resulting from exactly one reflection at a surface point $\mathbf{x}$ with direction $\boldsymbol{\omega}_o$ equals

where $\mathbf{y}$ represents the closest intersection of a light ray originating at $\mathbf{x}$ with direction $\boldsymbol{\omega}_i$, i.e., $\mathbf{y} = \operatorname{rayTrace}(\mathbf{x}, \boldsymbol{\omega}_i)$. Unlike RE $\eqref{eq:rendering-equation}$, which takes the form of an integral equation, Eq. $\eqref{eq:direct-illumination}$ is a simple spherical integral as its right-hand side involves only known quantities.

We now consider the problem of calculating the derivative of $L_r(\mathbf{x}, \boldsymbol{\omega}_o)$ with respect to the scene parameter vector $\boldsymbol{\pi}$. Given $\mathbf{x}$ and $\boldsymbol{\omega}_o$, let $f_{direct}(\boldsymbol{\omega}_i; \mathbf{x}, \boldsymbol{\omega}_o) := L_e(\mathbf{y}, -\boldsymbol{\omega}_i) \; f_s(\mathbf{x}, \boldsymbol{\omega}_i, \boldsymbol{\omega}_o)$. It holds that

By applying Reynolds transport theorem $\eqref{eq:reynolds-transport-theorem}$, we obtain:

where $\mathrm{d}\ell$ is the curve-length measure. This is exactly the RTT split from before, specialized to a spherical integral: the interior term integrates over the ($\boldsymbol{\pi}$-independent) sphere $\mathbb{S}^2$, and the boundary term picks up the jump of $f_{\text{direct}}$ across the 1D discontinuity curves $\Delta \mathbb{S}^2$, the silhouette-induced jumps in $L_e$, as they move with $\boldsymbol{\pi}$. For any $\boldsymbol{\omega}_i \in \mathbb{S}^2$, $\mathbf{n}^{\perp}(\boldsymbol{\omega}_i)$ is, as before, the tangent-space vector at $\boldsymbol{\omega}_i$ perpendicular to the discontinuity curve.

Assuming the (cosine-weighted) BSDF $f_s(\mathbf{x}, \boldsymbol{\omega}_i, \boldsymbol{\omega}_o)$ to be continuous with respect to $\boldsymbol{\omega}_i$, which is usually the case except for perfectly specular BSDFs, the discontinuities of the integrand $f_{direct}$ fully emerge from those of incident emission $L_e(\mathbf{y}, \boldsymbol{\omega}_i)$, which is generally discontinuous due to occlusions. Therefore,

Differential Rendering Equation

Based on the analysis above, we now differentiate the full rendering equation (RE) $\eqref{eq:rendering-equation}$ using the Reynolds transport theorem ($\eqref{eq:reynolds-transport-theorem}$). This yields another integral equation, which we call the differential rendering equation.

We begin by applying the derivative operator $\partial_{\boldsymbol{\pi}}$ to both sides of the standard rendering equation:

Since the scene geometry may move as the parameter $\boldsymbol{\pi}$ changes, the integration domain inherently contains moving boundaries (i.e., silhouettes). Applying the Reynolds Transport Theorem ($\eqref{eq:reynolds-transport-theorem}$) splits the derivative of this integral into an interior and a boundary component:

Note on the boundary term: Notice that we wrote the jump of the integrand as $f_s \Delta L_i$ rather than $\Delta(L_i f_s)$. This assumes that the cosine-weighted BSDF $f_s$ evaluates smoothly and continuously with respect to the incoming direction $\boldsymbol{\omega}_i$. For typical materials such as diffuse and rough microfacet models, this holds: the discontinuity is caused by the incoming radiance $L_i$ abruptly jumping when an integration ray sweeps past a silhouette edge or shadow boundary. The notable exception is perfectly specular materials, whose BSDFs are Dirac delta functions; handling those requires a different mathematical approach, such as attached sampling.

Inside the interior integral, we expand the derivative of the product $\partial_{\boldsymbol{\pi}} (L_i f_s)$ using the standard product rule:

Substituting this expansion back into the equation allows us to regroup the terms into two distinct transport components. We collect the terms acting as “sources” of differential radiance into one group, and the term representing scattered differential radiance into another:

This final equation shares the exact same structure as the original rendering equation. Instead of standard light emission and scattering, it describes the emission and scattering of differential radiance (gradients).

The differential emission term $Q(\mathbf{x}, \boldsymbol{\omega}_o)$ acts as the source of gradients. It evaluates to a non-zero value at any point where the primary emission changes ($\partial_{\boldsymbol{\pi}} L_e$), the material scattering properties change ($\partial_{\boldsymbol{\pi}} f_s$), or a silhouette edge moves to uncover a different object ($f_s \Delta L_i$). During Monte Carlo integration, we compute this local change $Q$ at every path vertex and add it to the running gradient estimate.

The differential scattering term $\int (\partial_{\boldsymbol{\pi}} L_i) f_s \mathrm{d}\sigma$ handles the propagation of these gradients. Here, $\partial_{\boldsymbol{\pi}} L_i$ represents the derivative of the incident radiance arriving from the previous bounce. Just like standard radiance, this incoming differential radiance is multiplied by the material’s BSDF ($f_s$) and scattered towards the camera.

Consequently, when a Monte Carlo path tracer simulates this process, it unrolls the recursion identical to forward rendering. For a path traced outward from the camera through vertices $\text{Camera} \to x_1 \to x_2 \to x_3$, the total gradient expands mathematically via the chain rule as:

In practice, this means we trace a standard light path and, at each bounce, compute the local differential emission $Q$, add it to the accumulated gradient, and multiply the running total by the surface BSDF as the path continues.

The differential rendering equation tells us what to compute; it says nothing about doing so efficiently. Evaluating this expansion naively ie., recording every bounce of every traced path onto an autodiff tape and replaying it backward, reintroduces exactly the memory and runtime blowup that made naive AD unsuitable for rendering in the first place (see Why is Differentiable Rendering Difficult?). For a light path of length $D$, that tape costs $\mathcal{O}(D)$ memory, and with millions of paths per frame it becomes the bottleneck long before the renderer does.

Efficient Reverse-Mode Differentiable Rendering

This section evaluates the differential rendering equation derived above, just without paying for the tape. Radiative Backpropagation and Path Replay Backpropagation both reformulate the backward pass as a second, physically-grounded transport simulation, so reverse-mode gradients can be computed with the same $\mathcal{O}(1)$-memory, single-pass character as forward rendering.

Radiative Backpropagation (Nimier-David et al., 2020)

Nimier-David et al. [7] introduce Radiative Backpropagation (RB). It begins from the same differential rendering equation derived above, but reorganizes reverse-mode differentiation as a second physical transport simulation. The objective is not a Jacobian image for one parameter. It is the vector-Jacobian product required by optimization: the derivative of one scalar objective with respect to all active scene parameters.

The Three Stages

Write the ordinary renderer as $\mathbf y=f(\boldsymbol{\pi})$ and the scalar objective as $g(\mathbf y)$. Radiative backpropagation separates one optimization iteration into:

1. an ordinary, non-differentiable render $\mathbf y=f(\boldsymbol{\pi})$;
2. differentiation of the comparatively small image-space objective, producing the adjoint rendering

3. a radiative-backpropagation simulation estimating

$\delta\mathbf y$ says how each rendered pixel should change to reduce the objective. The sensor emits this quantity into the scene as adjoint radiance. Scattering transports it to emitters, materials, and media whose local derivatives contribute to $\delta\boldsymbol{\pi}$.

In pseudocode:

def grad(x):
    # 1. Ordinary rendering (no AD)
    y = f(x)
    # 2. Differentiate objective at y (manually or w/ AD)
    δ_y = J_gᵀ(y)
    # 3. Estimate δ_x = J_fᵀ δ_y using radiative backpropagation
    return radiative_backprop(x, δ_y)

Adjoint Radiance and Operator Formulation

We begin directly from the differential rendering equation $\eqref{eq:diff-rendering-equation}$ derived above:

Expanding the integrand terms explicitly yields:

Note on Static Visibility Boundaries: Nimier-David et al. [7] never actually write down the boundary term above, the general, boundary-aware equation is machinery imported from the Zhang et al. [14] framework developed concurrently in the literature, not something the RB paper derives and then discards. The RB paper’s own derivation assumes static geometry from the outset ($\partial_{\boldsymbol{\pi}} \boldsymbol{\omega}_i = \mathbf{0}$), so for them the Boundary Integral is simply absent, leaving Direct Emission, Diff. Scattering, and Material Emission. The paper is explicit that this is a limitation of its prototype rather than something it resolves: visibility-related gradients are left to future work, pointing at Li et al. [10] and Loubet et al. [11] as compatible options (Section 3.6 of the paper).

Grouping the non-scattering gradient source terms into the Differential Emission term $Q(\mathbf{x}, \boldsymbol{\omega}_o)$:

the differential transport equation simplifies to:

This equation expresses an energy balance for differential radiance: $Q$ acts as a source of differential emission wherever a scene parameter directly alters emission or reflectance, while the integral describes how incident differential radiance undergoes standard physical scattering. Radiative backpropagation leverages this structure by evaluating gradients through an adjoint transport simulation governed by the same linear transport operators as primal rendering, with $L_e$ replaced by $Q$.

To make that reuse precise and to set up reverse-mode propagation, the paper packages the two remaining physical processes (scattering at a surface, and propagating along a ray to the next one) into two linear operators. Using Nimier-David et al.’s own notation [7] (Section 3.4):

1. Scattering operator $\mathcal{K}$. Takes an incident directional field $h$ and scatters it through the BSDF, exactly the way the ordinary scattering equation treats $L_i$:

   $$ (\mathcal{K} h)(\mathbf{x}, \boldsymbol{\omega}_o) = \int_{\mathbb{S}^2} h(\mathbf{x}, \boldsymbol{\omega}_i) \, f_s(\mathbf{x}, \boldsymbol{\omega}_i, \boldsymbol{\omega}_o) \,\mathrm{d}\sigma(\boldsymbol{\omega}_i). $$

2. Propagation operator $\mathcal{G}$. This is the paper’s own name for it; you’ll also see it called a ray transport operator, since all it does is walk backward along a ray to the next surface. It turns outgoing radiance at the point you hit into incident radiance at the point you came from:

   $$ (\mathcal{G} h)(\mathbf{x}, \boldsymbol{\omega}_i) = h(r(\mathbf{x}, \boldsymbol{\omega}_i), -\boldsymbol{\omega}_i), \qquad\text{so that}\qquad \partial_{\boldsymbol{\pi}} L_i = \mathcal{G}\, \partial_{\boldsymbol{\pi}} L_o. $$

Because differential radiance scatters and propagates exactly like ordinary radiance, $\mathcal{K}$ and $\mathcal{G}$ are the very same operators Veach used to analyze primal light transport (nothing new had to be invented here, which is precisely the point). Substituting $\partial_{\boldsymbol{\pi}} L_i = \mathcal{G}\partial_{\boldsymbol{\pi}} L_o$ into the scattering term folds the whole differential rendering equation into one compact line:

where $\mathcal{S} = \sum_{k=0}^\infty (\mathcal{K}\mathcal{G})^k$ sums over paths of every length, the operator equivalent of “trace a one-bounce path, then a two-bounce path, then a three-bounce path, and so on.”

If $A_e$ is the emitted adjoint radiance obtained by back-projecting the loss gradient $\delta\mathbf{y} = \mathbf{J}_g^T(\mathbf{y})$ from the sensor into the scene, the vector-Jacobian product we actually want is the ray-space inner product

For reciprocal, energy-conserving BSDFs, Veach [1997] established that the propagation operator $\mathcal{G}$, the scattering operator $\mathcal{K}$, and the composite transport operator $\mathcal{G}\mathcal{S}$ are self-adjoint under the ray-space measure: $\mathcal{G}^\ast = \mathcal{G}$, $\mathcal{K}^\ast = \mathcal{K}$, and $(\mathcal{G}\mathcal{S})^\ast = \mathcal{G}\mathcal{S}$. Note that $\mathcal{S}$ alone is not generally self-adjoint because $\mathcal{K}$ and $\mathcal{G}$ do not commute ($\mathcal{S}^\ast = (\mathcal{I}-\mathcal{K}\mathcal{G})^{-\ast} = (\mathcal{I}-\mathcal{G}\mathcal{K})^{-1} \ne \mathcal{S}$). However, the composite operator expands as $\mathcal{G}\mathcal{S} = \sum_k (\mathcal{G}\mathcal{K})^k \mathcal{G}$, which inherits self-adjointness directly from the individual symmetry of $\mathcal{G}$ and $\mathcal{K}$.

Using the self-adjointness of $\mathcal{G}\mathcal{S}$, the inner product shifts from forward differential emission to backward adjoint transport:

Here $A$ is the incident adjoint radiance field. This identity establishes the central computational advantage of the adjoint formulation: rather than propagating high-dimensional differential emission $Q$ forward through the scene for every parameter, a single scalar adjoint field $A_e$ is propagated backward from the sensor, and the parameter gradient $\delta\boldsymbol{\pi}$ is accumulated locally as an inner product $\langle A, Q \rangle$ at active surface interactions.

The corresponding incident and outgoing adjoint radiance satisfy the same recursive balance as ordinary light:

The paper obtains $A_e$ directly from the adjoint image: if $W_k(\mathbf{x},\boldsymbol{\omega}_o)$ is pixel $k$’s sensor importance and $\delta y_k$ is that pixel’s objective derivative,

turning the discrete sum over pixel derivatives into the ray-space inner product $\langle A_e,\partial_{\boldsymbol{\pi}}L_i\rangle$. For a pinhole camera, $A_e$ can be pictured as a textured “spotlight” that projects the adjoint image back into the scene from the camera.

Summary of Operator Roles: $\mathcal{K}$ represents local BSDF scattering, $\mathcal{G}$ performs ray propagation to the next surface, $\mathcal{S} = (\mathcal{I}-\mathcal{K}\mathcal{G})^{-1}$ computes the infinite path-length Neumann series, and $Q$ evaluates local parameter derivatives of emission and reflectance. Radiative backpropagation executes adjoint transport backward from the sensor ($A_i = \mathcal{G}A_o$, $A_o = A_e + \mathcal{K}A_i$), accumulating $\langle A, Q\rangle$ into $\delta\boldsymbol{\pi}$ at differentiable surface hits.

Sampling the Adjoint Transport Problem

These equations can be sampled by an ordinary path-tracing random walk launched from the sensor. At a surface hit, the algorithm:

• backpropagates the current adjoint weight through $L_e$ into active emitter parameters;
• samples a BSDF direction $\boldsymbol{\omega}_i$;
• backpropagates the weight $A_i L_i/p(\boldsymbol{\omega}_i)$ through the sampled BSDF value into active material parameters; and
• continues the adjoint path with throughput $f_s/p$.

The local reverse-mode operation is sparse. A texture lookup, for example, contributes only to nearby texels rather than constructing a dense derivative with respect to every scene parameter.

def radiative_backprop(π, δ_y):
    # Initialize parameter gradient(s) to zero
    δ_π = 0
    for _ in range(num_samples):
        # Importance sample a ray from the sensor
        x, ω_o, weight = sensor.sample_ray()
        # Evaluate the adjoint emitted radiance
        weight *= A_e(δ_y, x, ω_o) / num_samples
        # Propagate adjoint radiance into the scene
        δ_π += radiative_backprop_sample(π, x, ω_o, weight)
    # Finished, return gradients
    return δ_π


def radiative_backprop_sample(π, x, ω_o, weight):
    # Find an intersection with the scene geometry
    y = r(x, ω_o)
    # Backpropagate to parameters of emitter, if any
    δ_π = adjoint([[ L_e(y, -ω_o) ]], weight)
    # Sample a ray from the BSDF
    ω_i, bsdf_value, bsdf_pdf = sample_f_s(y, -ω_o, ·)
    # Backpropagate to parameters of BSDF, if any
    δ_π += adjoint([[ f_s(y, -ω_o, ω_i) ]], weight * L_i(y, ω_i) / bsdf_pdf)
    # Recurse
    return δ_π + radiative_backprop_sample(π, y, ω_i, weight * bsdf_value / bsdf_pdf)

Code: Radiative Backpropagation (RB) PyTorch Implementation

import torch
from scene import Scene
from camera import Camera
from ray import Ray


def _flip_normal(n, d):
    """Flip shading normal to face against the ray direction."""
    return torch.where((n * d).sum(-1, keepdim=True) > 0, -n, n)


def L_i(scene: Scene, x, wi, max_depth):
    """Estimate incoming radiance without building an autograd graph."""
    with torch.no_grad():
        L = torch.zeros_like(x)
        throughput = torch.ones_like(x)
        ray = Ray(x, wi)
        for _ in range(max_depth):
            si = scene.intersect(ray)
            valid = si.is_valid()
            n = _flip_normal(si.n, ray.dirs)

            L += torch.where(valid, throughput * si.emission, 0.0)
            wi, bsdf_value, bsdf_pdf = si.bsdf.sample(-ray.dirs, n)
            throughput = torch.where(valid, throughput * bsdf_value / bsdf_pdf, 0.0)
            ray = Ray(si.p + n * 1e-3, wi)
        return L


class RBPathTracer:
    """Radiative Backpropagation (Nimier-David et al. 2020)."""

    def __init__(self, max_depth=5, num_samples=128):
        self.max_depth = max_depth
        self.num_samples = num_samples

    def sample_path(self, scene: Scene, camera: Camera, seed: int = 42):
        """Primal render, independent of the adjoint pass."""
        torch.manual_seed(seed)
        accum = torch.zeros_like(camera.origins)

        for _ in range(self.num_samples):
            rays = camera.sample()
            accum += L_i(scene, rays.origins, rays.dirs, self.max_depth)

        return accum / self.num_samples

    def radiative_backprop_sample(self, scene: Scene, x, wo, weight):
        """radiative_backprop_sample(π, x, ω_o, weight), unrolled over max_depth bounces."""
        weight = weight.detach()
        ray = Ray(x, wo)

        for depth in range(self.max_depth):
            si = scene.intersect(ray)
            valid = si.is_valid()
            n = _flip_normal(si.n, ray.dirs)

            Le = torch.where(valid, si.emission, torch.zeros_like(si.emission))
            if Le.requires_grad:
                (Le * weight).sum().backward()

            if depth + 1 == self.max_depth:
                break

            wo = -ray.dirs.detach()
            wi, bsdf_value, bsdf_pdf = si.bsdf.sample(wo, n.detach())
            f_s = si.bsdf.eval(wo, n, wi.detach())
            y = si.p.detach() + n.detach() * 1e-3
            Li = L_i(scene, y, wi.detach(), self.max_depth - depth - 1)

            adjoint = torch.where(valid, weight * Li / bsdf_pdf.detach(), 0.0)
            if f_s.requires_grad:
                (f_s * adjoint.detach()).sum().backward()

            with torch.no_grad():
                weight = torch.where(valid, weight * bsdf_value / bsdf_pdf, 0.0)
                ray = Ray(si.p + n * 1e-3, wi)

    def radiative_backprop(self, scene: Scene, camera: Camera, dL):
        """radiative_backprop(π, δ_y): seed each sensor ray with weight = δ_y / num_samples."""
        for _ in range(self.num_samples):
            rays = camera.sample()
            weight = dL / self.num_samples
            self.radiative_backprop_sample(scene, rays.origins, rays.dirs, weight)

The paper does not claim that camera-path sampling is the only estimator. It points out that $Q$ may have large components at specific differentiable objects, suggesting connection strategies analogous to next-event estimation. If many of those connections are occluded, scattering them before connection leads to a family of bidirectional strategies. Casting reverse differentiation as transport is useful precisely because ordinary rendering tools such as importance sampling, next-event estimation, and bidirectional connection strategies become available.

Assumptions and Loose Ends

The derivation and prototype make several explicit assumptions:

• Static Sensor Importance: Sensor importance $W_k$ is treated as static. A differentiable camera would contribute an additional local derivative term.
• Self-Adjoint Operators: Surface operators are self-adjoint under reciprocal, energy-conserving BSDFs and the ray-space measure. Nonreciprocal transport (e.g. non-reciprocal BSDFs or camera lenses) would require true adjoint operators rather than reusing primal ones.
• Volumetric Extension: The same adjoint construction extends directly to participating media by replacing surface transport operators with their volumetric counterparts (RB paper, Appendix A.1).
• Visibility Derivatives (Edge Sampling vs. Reparameterization):
  • Edge Sampling (Li et al., 2018) [10]: Explicitly samples silhouette edges, which requires modifying the theoretical formulation to add extra 1D boundary integral terms to the differential emission source $Q$.
  • Reparameterization / Change of Variables (Loubet et al., 2019) [11]: Performs a parameter-dependent coordinate warp so that discontinuities remain static under scene perturbations ($\partial_{\boldsymbol{\pi}}\boldsymbol{\omega}_i = \mathbf{0}$). This allows RB to compute visibility-aware gradients without changing any of the adjoint derivations.

At an intersection $\mathbf{y} = r(\mathbf{x}, \boldsymbol{\omega}_o)$, an adjoint path contributes two local VJPs corresponding to the terms of $Q$:

and, for $\boldsymbol{\omega}_i \sim p(\boldsymbol{\omega}_i)$,

The path then continues with adjoint throughput multiplied by $f_s/p$. This is the Monte Carlo realization of $\langle\mathcal G\mathcal S A_e,Q\rangle$, not a reversal of stored primal vertices.

Primal Radiance Dependence and Acceleration

A key practical bottleneck is that the differential emission term $Q$ depends on the unknown primal incident radiance $L_i$. Evaluating $Q$ at every differentiable interaction requires launching a recursive primal path-tracing query. Along a path of depth $D$ with differentiable surfaces at each bounce, these suffix queries have lengths $D, D-1, \ldots, 1$, leading to quadratic time complexity $\mathcal{O}(D^2)$ (a bottleneck also reported by Zhang et al. (2019) [14] for forward AD). Non-differentiated interactions do not trigger this extra work.

To mitigate this quadratic overhead in long light paths, one can precompute an approximate spatio-directional data structure during the primal phase (such as a Path Guiding tree [Müller et al., 2017] [15]) to perform fast $\mathcal{O}(1)$ interpolant queries of $L_i$ during adjoint backpropagation.

Thus, RB solves the reverse-mode storage and transport problem; it is not by itself a solution to moving visibility discontinuities unless paired with reparameterization or boundary sampling.

Biased Variants

The paper studies two faster approximations:

• Biased I: replace $L_i$ in $Q$ by $1$. This removes recursive radiance queries and reduces time to $O(D)$, but changes the gradient.
• Biased II: pipeline optimization by using the previous iteration’s adjoint rendering with the current iteration’s rendering Jacobian. This overlaps primal and adjoint work but introduces an intentional one-iteration mismatch.

The RB paper originally claimed that Biased I preserves gradient signs. The authors’ published errata withdraw this claim and defer to Section 3.2 of Vicini et al. [13] for a detailed explanation: each local product preserves its sign when multiplied by positive radiance, but globally summing many differently weighted signed contributions need not preserve the sign. Biased I can still work well and often reduces variance, but it must be described as a heuristic rather than a directionally correct gradient estimator.

For Biased II, if superscripts denote optimization iterations, the propagated quantity is

instead of using $\delta\mathbf y^{(i)}$. It is reasonable only when the rendering Jacobian and objective gradient vary slowly between iterations. The paper labels the combined approximation Biased I + II.

Limitations

While the unbiased algorithm requires constant memory with respect to path length, repeated evaluation of $L_i$ leads to quadratic time complexity. Specifically, along an adjoint path of depth $D$ where every bounce interacts with a differentiable surface, evaluating $Q$ requires launching a recursive primal path-tracing query of depth $D - d$ at each bounce $d$. Summing over all bounces yields $(D - 1) + (D - 2) + \dots + 1 = \mathcal{O}(D^2)$ ray intersections. Although one could attempt to mitigate this overhead by randomly sampling a single evaluation branch per bounce, the resulting estimator suffers from exponential variance growth over deep paths. This is prohibitive for highly scattering media with thousands of events. Biased I reduces that cost to linear time but is not a correct derivative. The formulation above omits moving-visibility derivatives and does not handle derivatives through ideal specular BSDF sampling. Faster gradient evaluation also does not remove nonconvexity or poor conditioning from the inverse problem itself.

Path Replay Backpropagation (Vicini et al., 2021)

Radiative backpropagation achieves a constant memory footprint by computing a fresh primal suffix at each differentiable interaction. However, this nested recursion causes computation time to grow quadratically ($\mathcal{O}(D^2)$) with the number of scattering events.

Vicini et al. [13] propose an elegant alternative called Path Replay Backpropagation (PRB). By leveraging the mathematical invertibility of local light transport Jacobians, PRB computes exact gradients in linear time ($\mathcal{O}(D)$) and constant memory ($\mathcal{O}(1)$).

PRB splits gradient evaluation into two separate passes:

1. Primal Pass: Light paths are sampled as usual, but instead of building a massive automatic differentiation (AD) graph, the renderer only records the total path radiance and the random seed.
2. Adjoint Replay Pass: The random sequence is replayed to trace the exact same path. As the path unfolds, local derivatives are backpropagated to the scene parameters on the fly by dynamically reconstructing the incident illumination.

Dynamic Suffix Reconstruction

To evaluate exact gradients without incurring a prohibitive memory overhead, the adjoint replay pass must accurately reconstruct the incident illumination suffix arriving at each vertex. Explicitly storing this information per vertex would require $\mathcal{O}(D)$ memory. Instead, PRB algebraically recovers this suffix on the fly in $\mathcal{O}(1)$ time by sequentially peeling off emitted contributions.

Consider the total accumulated radiance $L_N$ over a path of $N$ vertices:

During the adjoint replay pass, the exact same sequence of vertices is visited in forward order. At the first vertex ($k=1$), the current suffix $L_{\text{current}}$ is obtained by subtracting the local emission $L_{e,1}$ from the total radiance:

Progressing to the second vertex ($k=2$), the subsequent emission $\beta_1 L_{e,2}$ is subtracted to yield the next suffix:

In general, the suffix dynamically tracks the remaining path radiance by sequentially removing the local emission at each step $k$:

This tracking variable maps directly to the L_reconstructed = L_total - throughput * L_e(...) operation within the adjoint pseudocode.

To formally connect this algebraic tracking variable to incident illumination, we first define the physical incident radiance $L_{i,k}$ actually arriving at vertex $k$. It is the sum of all future emissions, weighted by the relative scattering throughput from that point onward:

If we multiply this physical incident radiance by the accumulated path throughput up to and including the scattering at $k$ (where $\beta_k = \beta_{k-1} \frac{f_k}{p_k}$), we project it into sensor space. This yields the remaining path radiance $L_k$:

Notice that this expression for $L_k$ is mathematically identical to our algebraically reconstructed suffix $L_{\text{current}}$. Therefore, we establish the direct relationship:

Proof: Exactness of Constant-Memory Path Replay Backpropagation

We seek the gradient of the scalar image loss $\mathcal{J}$ with respect to a specific BSDF evaluation $f_k$. By the chain rule:

$$ \frac{\partial \mathcal{J}}{\partial f_k} = \delta L \frac{\partial L_N}{\partial f_k} $$

where $\delta L$ is the incoming adjoint radiance from the sensor.

Because path throughputs prior to bounce $k$ ($j \leq k$) are independent of $f_k$, their derivatives vanish identically. Applying the derivative strictly to the subsequent terms of the total radiance sum yields:

$$ \frac{\partial L_N}{\partial f_k} = \frac{\partial}{\partial f_k} \sum_{j=k+1}^{N} \beta_{j-1} L_{e,j} $$

Since each downstream throughput $\beta_{j-1}$ depends linearly on $f_k$, its partial derivative is simply $\frac{\beta_{j-1}}{f_k}$. Factoring out $\frac{1}{f_k}$ isolates the reconstructed suffix $L_{\text{current}}$:

$$ \frac{\partial L_N}{\partial f_k} = \frac{1}{f_k} \sum_{j=k+1}^{N} \beta_{j-1} L_{e,j} = \frac{1}{f_k} L_{\text{current}} $$

Substituting the identity $L_{\text{current}} = \beta_{k-1} \frac{f_k}{p_k} L_{i,k}$ into the chain rule expression:

$$ \frac{\partial \mathcal{J}}{\partial f_k} = \delta L \cdot \frac{1}{f_k} \left( \beta_{k-1} \frac{f_k}{p_k} L_{i,k} \right) = \delta L \cdot \beta_{k-1} \frac{L_{i,k}}{p_k} $$

The factor $f_k$ cancels out. This algebraic result proves that dynamically tracking $L_{\text{current}}$ and dividing by $f_k$ evaluates the exact gradient multiplier on the fly, eliminating the need to construct or retain an automatic differentiation graph.

Continuous Path Unrolling and Intermediate Gradient Flow

This discrete index cancellation has a direct physical analogue in the continuous formulation of light transport. Differentiating outgoing radiance $L_o$ at a surface point $\mathbf{x}_0$ splits the integral into two components via the product rule:

$$ \partial_{\boldsymbol{\pi}} L_o(\mathbf{x}_0, \boldsymbol{\omega}_0) = \int_{\mathbb{S}^2} \Bigg[ {\color{#3b82f6}\underbrace{(\partial_{\boldsymbol{\pi}} f_s) \cdot L_i}_{\text{Term A}}} + {\color{#ff6b6b}\underbrace{f_s \cdot (\partial_{\boldsymbol{\pi}} L_i)}_{\text{Term B}}} \Bigg] \mathrm{d}\boldsymbol{\omega}_1 $$

Here, ${\color{#3b82f6}\text{Term A}}$ captures the local material derivative under unperturbed incident illumination, while ${\color{#ff6b6b}\text{Term B}}$ accounts for the recursive change in incoming radiance.

To see why detaching the incident illumination state during local backpropagation does not lose multi-bounce gradient flow, consider unrolling a two-bounce path ($\mathbf{x}_0 \to \mathbf{x}_1 \to \mathbf{x}_2$):

$$ L_o(\mathbf{x}_0, \boldsymbol{\omega}_0) = \int_{\mathbb{S}^2} f_s(\mathbf{x}_0, \boldsymbol{\omega}_1, \boldsymbol{\omega}_0) \left[ \int_{\mathbb{S}^2} f_s(\mathbf{x}_1, \boldsymbol{\omega}_2, -\boldsymbol{\omega}_1) L_e(\mathbf{x}_2, -\boldsymbol{\omega}_2) \, \mathrm{d}\boldsymbol{\omega}_2 \right] \mathrm{d}\boldsymbol{\omega}_1 $$

Differentiating the unrolled integral under the assumption of a static emitter ($\partial_{\boldsymbol{\pi}} L_e = 0$) evaluates ${\color{#ff6b6b}\text{Term B}}$ at $\mathbf{x}_0$, for a single fixed direction $\boldsymbol{\omega}_1$, as:

$$ {\color{#ff6b6b}\text{Term B}_{\mathbf{x}_0}} = f_s(\mathbf{x}_0) \int_{\mathbb{S}^2} {\color{#3b82f6}\underbrace{\partial_{\boldsymbol{\pi}} f_s(\mathbf{x}_1) \cdot L_e(\mathbf{x}_2)}_{\text{Term A at vertex } \mathbf{x}_1}} \mathrm{d}\boldsymbol{\omega}_2 $$

Notice the equivalence: ${\color{#ff6b6b}\text{Term B}}$ at vertex $\mathbf{x}_0$ is identically equal to ${\color{#3b82f6}\text{Term A}}$ at the subsequent vertex $\mathbf{x}_1$, scaled by the local throughput $f_s(\mathbf{x}_0)$.

Evaluating ${\color{#ff6b6b}\text{Term B}}$ recursively at $\mathbf{x}_0$ is therefore redundant: when the adjoint replay pass advances to $\mathbf{x}_1$, evaluating ${\color{#3b82f6}\text{Term A}_{\mathbf{x}_1}}$ with accumulated throughput $\beta_0 = f_s(\mathbf{x}_0)$ automatically computes the exact contribution required by ${\color{#ff6b6b}\text{Term B}_{\mathbf{x}_0}}$. Consequently, PRB detaches the illumination state at each bounce without discarding any downstream gradient signal.

---

Implementation in Adjoint Replay

In practical code, these mathematical properties translate directly into two localized operations per bounce:

1. Dynamic Suffix Peeling: The total radiance accumulator is decremented by local emission, L = L - Le.detach(), maintaining $L_{\text{current}} = \beta_k \frac{f_k}{p_k} L_{i,k}$.
2. Local Backward Step: Dividing $L_{\text{current}}$ by $f_k$ via relative_grad(f_s) isolates $\beta_{k-1} \frac{L_{i,k}}{p_k}$, accumulating ${\color{#3b82f6}\text{Term A}_k}$ directly into parameter gradients with zero graph retention.

Algorithmic Implementation (Detached PRB)

This algebraic collapse maps elegantly into two constant-memory passes.

Primal Phase: The first phase determines the total radiance L accumulated by the path without tracking gradients.

# Pass 1: Primal path tracing pass
def sample_path(ray):
    L = 0
    β = 1
    for i in range(N):
        L += β * L_e(...)
        ω_i, bsdf_value, bsdf_pdf = sample_bsdf(...)
        β *= bsdf_value / bsdf_pdf
    return L

Adjoint Phase: The second phase replays the exact same random walk. By dynamically reconstructing the required radiance suffix on the fly, PRB accumulates the local parameter gradients directly in constant memory, completely bypassing the need to store a global AD computation graph.

# Pass 2: Adjoint replay pass
def sample_path_adjoint(ray, L, δL):
    β = 1
    δ_π = 0
    for i in range(N):
        L -= β * L_e(...)
        ω_i, bsdf_value, bsdf_pdf = sample_bsdf(...)
        δ_π += backward_grad(bsdf_value, δL * L / bsdf_value)
        β *= bsdf_value / bsdf_pdf
    return δ_π

Code: Path Replay Backpropagation (PRB) PyTorch Implementation

import torch
from scene import Scene
from camera import Camera
from ray import Ray


def relative_grad(x, eps=1e-10):
    x_d = x.detach()
    safe = x_d.abs() > eps
    denom = torch.where(safe, x_d, torch.ones_like(x_d))
    return torch.where(safe, x / denom, torch.zeros_like(x))


def _flip_normal(n, d):
    """Flip shading normal to face against the ray direction."""
    return torch.where((n * d).sum(-1, keepdim=True) > 0, -n, n)


class PRBPathTracer:
    """Path Replay Backpropagation (Vicini et al. 2021)."""

    def __init__(self, max_depth=5, num_samples=128):
        self.max_depth = max_depth
        self.num_samples = num_samples

        self.seed = 42
        self._primal_samples = []

    def sample_path(self, scene: Scene, camera: Camera, seed: int = 42):
        self.seed = seed
        torch.manual_seed(seed)

        self._primal_samples.clear()
        accum = torch.zeros_like(camera.origins)

        with torch.no_grad():
            for _ in range(self.num_samples):
                ray = camera.sample()
                L = torch.zeros_like(ray.origins)
                throughput = torch.ones_like(ray.origins)

                for _ in range(self.max_depth):
                    si = scene.intersect(ray)
                    valid = si.is_valid()
                    n = _flip_normal(si.n, ray.dirs)

                    L += torch.where(valid, throughput * si.emission, 0.0)
                    wi, bsdf_value, bsdf_pdf = si.bsdf.sample(-ray.dirs, n)
                    throughput = torch.where(valid & (bsdf_pdf > 1e-8), throughput * bsdf_value / torch.clamp(bsdf_pdf, min=1e-8), 0.0)
                    ray = Ray(si.p + n * 1e-3, wi)

                self._primal_samples.append(L)
                accum += L

        return accum / self.num_samples

    def sample_adjoint(self, scene: Scene, camera: Camera, _primal_img, dL):
        torch.manual_seed(self.seed)

        scale = dL / self.num_samples

        for sample_idx in range(self.num_samples):
            ray = camera.sample()
            ray = Ray(ray.origins.detach(), ray.dirs.detach())
            L = self._primal_samples[sample_idx]
            throughput = torch.ones_like(ray.origins)

            for _ in range(self.max_depth):
                si = scene.intersect(ray)
                valid = si.is_valid()
                n = _flip_normal(si.n, ray.dirs)

                # L -= β · L_e  ->  L is now suffix radiance R_k
                Le = throughput.detach() * si.emission
                Le = torch.where(valid, Le, torch.zeros_like(Le))
                L = L - Le.detach()

                # same random stream -> identical (wi, w)
                wi, bsdf_value, bsdf_pdf = si.bsdf.sample(-ray.dirs, n)

                # differentiable f_s re-evaluation
                f_s = si.bsdf.eval((-ray.dirs).detach(), n, wi.detach())

                # dπ += J_{Le}^T(dL)  +  J_{f_s}^T(dL * R_k / f_s)
                Lo = Le + torch.where(valid, L * relative_grad(f_s), 0.0)
                (scale.detach() * Lo).sum().backward()

                # advance (fully detached)
                with torch.no_grad():
                    throughput = torch.where(valid & (bsdf_pdf > 1e-8), throughput * bsdf_value / torch.clamp(bsdf_pdf, min=1e-8),
                                             torch.zeros_like(throughput))
                    ray = Ray((si.p + n * 1e-3).detach(), wi.detach())

Iterative Jacobian Inversion

The replay principle has an elegant algebraic interpretation in terms of loop state transitions and local Jacobian inversion.

Consider path tracing as the repeated composition of a loop transition function $h$. In a basic path tracer, the loop state is $\mathbf{z} = (L, \beta)$ initialized to $\mathbf{z}_0 = (0, 1)$. At step $k$, $h$ updates the state:

(omitting the division by sampling density for clarity). After $N$ bounces, the accumulated state is the $N$-fold composition:

Differentiating this composition with respect to scene parameters $\boldsymbol{\pi}$ via the multivariate chain rule gives:

There are three architectural ways to evaluate this chain-rule expression:

1. Conventional Reverse-Mode AD: Store all intermediate activations $(L_k, \beta_k)$ in memory and backpropagate in reverse order (requiring $\mathcal{O}(D)$ memory).
2. Primal Program Inversion (Reversible Computing): Invert the entire forward computation to re-evaluate states during the backward pass (as in Reversible Residual Networks [Gomez et al. 2017]).
3. Local Jacobian Inversion (PRB): Rather than inverting the primal program, invert the low-dimensional local Jacobian matrix $J_h$ relating adjacent loop states.

For the transition function $h$, the single-step Jacobian is:

Multiplying these matrices over the suffix from vertex $k+1$ to $N$ yields the suffix Jacobian product:

where $L_{k,N}$ and $\beta_{k,N}$ denote the accumulated incident radiance and throughput from vertex $k$ to $N$. This reveals that the physical indirect illumination suffix is literally an entry of the accumulated Jacobian product.

Expanding the parametric derivative $\partial_{\boldsymbol{\pi}} h(\boldsymbol{\pi}, L_{k-1}, \beta_{k-1}) = \begin{pmatrix} \beta_{k-1}\partial_{\boldsymbol{\pi}} L_e \\ \beta_{k-1}\partial_{\boldsymbol{\pi}} f_s \end{pmatrix}$ into the chain rule yields:

In rendering, only the first component ($\partial_{\boldsymbol{\pi}} L$) is needed. Extracting the top row gives:

During the adjoint replay pass, as we sequentially subtract the current emitted radiance and divide by $f_s$, we are iteratively applying the inverse Jacobian matrix:

Because $J_h$ is only $2\times 2$ (or $4\times 4$ in attached rendering with ray differential Jacobians), inverting it is exact, numerically stable, and computationally trivial.

Attached Sampling and Specular Paths

Detached differentiation cannot optimize parameters that move ideal specular samples, such as the index of refraction, geometry, or normals of smooth dielectrics and conductors. To include these dependencies, write inverse-transform sampling as a parameter-dependent map $\mathbf{x} = T(\mathbf{u}, \boldsymbol{\pi})$ from $\mathbf{u} \in [0, 1]^n$ to path space. The reparameterized pixel integral is

Differentiating while keeping $\mathbf{u}$ fixed gives the attached estimator

The numerator contains both the explicit material derivative and the motion of the sampled path:

A perturbation at one interaction therefore changes all later vertices and their BSDF and emission factors. PRB reconstructs this non-local dependence using the differential relationship between adjacent path segments. A ray can be represented by two local surface coordinates at its origin and two at its endpoint, reducing the Jacobian between adjacent segments to a $4 \times 4$ matrix. Explicitly propagating these small ray, throughput, and radiance Jacobians carries all derivative information between interactions without storing a path-length AD graph. A position-position parameterization keeps the matrix entries dimensionally compatible and is better conditioned than mixing positions and angles.

# Pass 1: Primal path tracing pass (Attached)
def sample_path(ray):
    L = 0
    β = 1
    J_L = 0_3,4
    J_β = 0_3,4
    J_ray = I_4
    for i in range(N):
        L += β * L_e(...)
        ω_i, bsdf_value, bsdf_pdf = sample_bsdf(...)
        bsdf_weight = bsdf_value / bsdf_pdf
        ray_prime = spawn_ray(ω_i, ...)

        # Compute the directional radiance derivative
        J_ray_prime, J_bsdf, J_L_e = forward_grad(ray, {ray_prime, bsdf_weight, L_e})
        J_bsdf = J_bsdf @ J_ray
        J_L_e  = J_L_e  @ J_ray
        J_ray  = J_ray_prime @ J_ray
        J_L   += β * J_L_e + L_e * J_β
        J_β    = bsdf_weight * J_β + β * J_bsdf
        β     *= bsdf_weight
    return L, J_L


# Pass 2: Adjoint replay pass (Attached)
def sample_path_adjoint(ray, L, J_L, δL):
    β = 1
    J_ray = I_4
    δ_π = 0
    for i in range(N):
        L -= β * L_e(...)
        ω_i, bsdf_value, bsdf_pdf = sample_bsdf(...)
        bsdf_weight = bsdf_value / bsdf_pdf
        ray_prime = spawn_ray(ω_i, ...)

        J_ray_prime, J_bsdf, J_L_e = forward_grad(ray, {ray_prime, bsdf_weight, L_e})
        J_bsdf = J_bsdf @ J_ray
        J_L_e  = J_L_e  @ J_ray
        J_ray  = J_ray_prime @ J_ray

        # Update the directional radiance derivative
        J_L -= L / bsdf_weight * J_bsdf + β * J_L_e
        J_L_prime = J_L @ (J_ray)^-1

        # Backpropagate gradients of the current BSDF value
        δ_π += backward_grad(bsdf_weight, δL * L / bsdf_weight)
        # Backpropagate through shading frame and BSDF sampling calculation
        δ_π += backward_grad(ray_prime, δL @ J_L_prime)

        β *= bsdf_weight
    return δ_π

Stochastic Regularization and Moving Discontinuities

The adjacent-ray Jacobian is singular when a sampling map does not depend on the incident segment, as with diffuse scattering. PRB regularizes it with

The perturbation has zero mean, and the same draw is reused by the correlated evaluations. This regularizes the inversion while preserving the estimator expectation.

Attached sampling may also turn a static discontinuity into one that moves with $\boldsymbol{\pi}$. Replay does not account for the resulting boundary term. An auxiliary reparameterization can slow or stop the sampling-map motion near discontinuities; alternatively, a boundary estimator must supply the missing term if unbiased geometry derivatives are required.

Differentiable Delta Tracking (Volumetric Rendering)

Path replay backpropagation is especially crucial for volumetric transport in participating media (e.g., clouds, smoke, tissue), where paths easily reach thousands of scattering interactions.

The $\mathcal{O}(D^2)$ complexity of standard Radiative Backpropagation is computationally prohibitive for these deep volumes. Furthermore, unbiased null-collision methods (like delta tracking) introduce discrete random decisions regarding real vs. fictitious collisions. By applying PRB to the volumetric radiative transfer equation, PRB isolates gradients with respect to the continuous absorption and scattering coefficients, ignoring the discontinuous sampling decisions. This unlocks unbiased volumetric derivatives in strictly linear time.

Complexity and Scope

The computational complexity and capabilities of the Path Replay Backpropagation framework are summarized below. For a light path of length $D$:

None of the rows above account for moving visibility discontinuities on their own; that term has to come from one of the boundary estimators of the previous section.

PRB strictly removes the path-length memory and time bottlenecks, successfully bringing the computational cost of unbiased differentiable rendering down to match that of standard forward path tracing.

Conclusion

Differentiable rendering bridges the gap between physics-based light transport and gradient-based optimization. As we have seen, taking the derivative of a rendering algorithm is far from trivial, and the field’s progress splits cleanly along two axes.

The first is correctness in the presence of moving discontinuities. The Reynolds transport theorem says what the missing gradient is; the disagreement is over how to estimate it. Edge sampling finds the silhouette and integrates over it directly, reparameterization and warped-area sampling refuse to look for it and absorb it into the interior instead (the latter fixing the former’s bias via the divergence theorem), and projective sampling recovers it from samples the renderer was already generating.

The second is cost. Once the differential rendering equation makes the recursion explicit, the problem becomes propagating gradients through arbitrarily long paths without recording them. Radiative Backpropagation reframes that propagation as a second transport simulation and buys $\mathcal{O}(1)$ memory at the price of quadratic time; Path Replay Backpropagation recovers linear time by inverting the local loop Jacobian instead of storing it.

While the foundational theory operates on standard geometric representations and path-tracing operators, translating these concepts into robust, scalable, and low-variance algorithms remains an exciting area of active research.

Further Reading

For those interested in exploring state-of-the-art developments and modern inverse rendering frameworks, consider the following resources:

• Mitsuba 3 Documentation: A highly flexible, retargetable forward and differentiable renderer. Read the docs
• Differentiable Signed Distance Function Rendering (Vicini et al., 2022) [3]: Optimizing geometry using implicit SDFs to allow topological changes. DOI: 10.1145/3528223.3530139
• A Simple Approach to Differentiable Rendering of SDFs (Wang et al., 2024) [4]: Exchanging unbiasedness for low variance and structural simplicity in SDF rendering. DOI: 10.1145/3680528.3687573
• Many-Worlds Inverse Rendering (Zhang et al., 2025) [6]: Avoiding local minima by evaluating a superposition of independent surface hypotheses. DOI: 10.1145/3767318