Differentiable Rendering with Exact Replay

1 Introduction

Users judge a renderer first by its images. They judge a differentiable renderer also by its gradients, which say how each pixel would change if a material, a pose or a dimension changed. Those gradients feed optimisers, so an error in them does not show as a wrong pixel. It shows as an optimiser that fails to converge, and the user cannot tell whether the estimator is biased, the device rounded differently, or the samples were noisy.

Path replay backpropagation ¹ made physically based differentiable rendering practical by replacing the stored path tape with a second, replayed pass. That only works if the replay regenerates exactly the paths the forward pass drew; otherwise the gradient belongs to a different estimator from the one rendered. Dr.Jit ² meets this contract by replaying with the same seeded generator, and accumulates parameter gradients with atomic adds whose order can vary between runs. We wanted a GPU tier and an f64 CPU reference that run the same material source and agree to rounding, with each tier bit-identical across runs. A gradient test with fixed random numbers then points at the estimator, not at the replay.

This paper describes irys, the renderer behind vcad and our robotics tools, and how it keeps that contract. Its contributions are:

• One BSDF source, two lowerings, two tiers. Materials are tang tensor graphs. tang lowers each one to an evaluation program and an adjoint program, emitting WGSL kernels for the GPU and running the same graph in its f64 interpreter for the CPU reference. The path logic (ray generation, traversal, bounce, accumulation) is hand-written (Section 3).
• PRB on a deterministic wavefront with hero wavelengths and NEE. The estimator is path replay backpropagation ¹ with next-event estimation, as in Mitsuba 3; what we add is a wavefront whose replay and gradient are bit-identical across runs, with no float atomics, with four wavelengths per path and the light-sampled lanes folded into the same graph batch (Section 4).
• Bit-exact replay by construction. The replay is exact because the integer RNG is exact on both tiers, reductions run in a fixed order with no float atomics, and the BVH walk depends only on the ray (Section 5).
• Pose gradients through instanced CAD meshes. irys samples primary edges with crease ownership and visibility-aware importance, differentiates the interior term by forward-mode duals with attached frames, and samples secondary edges for moving shadows and reflections (Section 6).
• Validation at every layer. We check gradients against central differences with common random numbers, and the GPU against the CPU’s f64, on boxes, on a shaded robot arm read from its CAD source, and on the shadows it casts (Section 7).
• Gradients for CAD parameters. irys computes a per-vertex image gradient that the vcad BRep kernel pulls back to the dimensions of the part. vcad’s paper ³ shares it as a benchmark (Section 7.5).

1.1 Scope and limits

• Determinism is per device and driver. The GPU tier runs in f32, so its results agree with the f64 reference to rounding, not bitwise; it is bit-identical against itself on one device and driver.
• Materials. The differentiable materials are a GGX conductor with a three-band spectral Fresnel, a diffuse surface’s three-band spectral albedo, and each Gaussian splat’s ⁴ DC colour and base opacity (Section 7.3). irys does not differentiate the environment or a splat’s position or shape. On the GPU a replay records at most 256 splat crossings per path.
• Lights. The emitters are the environment (a sky and a von Mises–Fisher sun lobe) and area emitters, sampled with next-event estimation and multiple importance sampling on both tiers.
• Pose parameters are rigid. Joint angles and the camera move instances rigidly. Non-rigid shape change enters only through the BRep pullback of Section 7.5.
• No comparison with Mitsuba 3 here. Every check in this paper is internal, against irys’s own f64 tier and central differences. The repository has an earlier comparison with Mitsuba 3 (bench/mitsuba), measured on different hardware and not repeated in this paper’s run.
• Hardware. The GPU tier runs on Metal and on Vulkan. Every number in this paper comes from one run on one machine: an NVIDIA RTX 4090 through Vulkan, with 12 vCPUs of an AMD EPYC 7K62 for the CPU tier (Section 9). We have not measured timings on other GPUs.

2 Background and related work

Differentiable path tracing. Li et al. introduced edge sampling to account for visibility discontinuities in Monte Carlo rendering ⁵. Radiative backpropagation ⁶ and path replay backpropagation ¹ reduced the memory of the adjoint pass to a constant per path. Zeltner et al. classify estimators for differential transport as attached or detached ⁷, the vocabulary we use. Dr.Jit provides a tracing compiler for such methods ², on which Mitsuba 2 and 3 are built ^8,9. Path-space differentiable rendering derives the boundary term as an integral over path space ¹⁰. Rasterization-based differentiable renderers trade global illumination for speed ¹¹. Other methods handle visibility by reparameterisation ¹², warped-area sampling ¹³ or projective sampling ¹⁴. We compared warped estimators for the secondary term and chose edge sampling (Section 6).

Spectral rendering. Hero wavelength sampling ¹⁵ traces several wavelengths along one path, the hero choosing directions. We use it unchanged and make each wavelength its own PRB lane.

Compilers for rendering. SLANG.D adds automatic differentiation to a shading language ¹⁶ and is the closest analogue to writing BSDFs as tang graphs. Laine et al. argue for splitting a path tracer into wavefront kernels ¹⁷, the structure we use.

Differentiable rendering systems. Mitsuba 3 ⁹ implements path replay with next-event estimation and multiple importance sampling, projective sampling and reparameterisation on Dr.Jit, and is the closest system to irys’s differentiable core. irys differs in what it guarantees rather than in the estimator: a replay that is bit-identical to the forward pass by construction, an f64 reference that runs the same material source, spectral transport with hero wavelengths, and visibility terms adapted to CAD tessellations (coplanar slivers, crease ownership, attached frames).

Inverse CAD and measurement from images. Recovering an editable CAD description from observations is usually posed as program synthesis over meshes, point clouds or drawings ^18–20, and close-range photogrammetry measures manufactured parts by triangulating targets and features ²¹. For photographs, irys adds what both lack, the gradient of a rendered image with respect to a part’s pose, finish and dimensions. Our companion paper’s BRep kernel ³ differentiates the CAD parameters upstream of the mesh, and Section 7.5 joins it to irys.

3 Materials as compiled graphs

The differentiable material has parameters : GGX roughness with height-correlated Smith masking, and a Schlick Fresnel whose is interpolated between three spectral hats. We write it once as a tang TensorGraph whose inputs are per-lane columns (the cosines at a vertex and the hero wavelength’s hat weights) and whose is broadcast. tang lowers it two ways. The eval program returns and its pdf. The adjoint program is the backward pass of , so the batch sum over lanes is the reduction to . tang emits each program as WGSL, one kernel per materialised comprehension, or as compiled C in f32 and an interpreter in f64 for the CPU. The interpreter matches a hand-written f64 GGX to , and the compiled C matches the interpreter to in the worst relative error.

The path logic is hand-written, in WGSL for the GPU and in Rust for the CPU reference: ray generation, BVH traversal, bounce, the throughput update, the final splat and accumulation, and the device-side loss scoring. tang’s comprehensions are straight-line bodies over fixed shapes, so they cannot express data-dependent path loops or a traversal stack today. Sampling the GGX visible normals is detached, since it needs no derivative under PRB, and the divisor pdf is the sampler’s, which we test against the graph’s pdf to .

4 Path replay on a wavefront

Estimator. For a path with vertices and radiance , BSDF-only PRB accumulates at each vertex, where is the radiance collected before vertex . With next-event estimation and multiple importance sampling the gradient becomes

where is the BSDF-sampled direction, the light-sampled one, the NEE contribution at and the radiance up to and including it. We detach the pdfs in the MIS weights, so the weights are constants under ; this keeps the estimator unbiased because the weights still sum to one for every path. In the wavefront the light-sampled lanes are the second half of one graph batch, so a single adjoint dispatch serves both strategies. A nonlinear loss needs the primal image before the adjoint, so irys scores the loss between the forward and replay passes.

Hero wavelengths. Each path carries four stratified wavelengths, the hero and its rotations by of the visible range ¹⁵. Only the hero draws directions. The wavefront’s lanes are (path, wavelength) pairs, four times as many as paths, and each lane is its own PRB path. At equal paths, four wavelengths cut per-pixel variance 7.4× in total and 28.7× in chroma. That costs four times the lanes, so per lane the gain is 1.9× in total and 7.2× in chroma.

Memory. PRB stores nothing per vertex. With four wavelengths a path costs about 71 f64 words at depth 2 and at depth 8 alike. Memory scales with paths, not pixels. irys traces frames larger than one dispatch in chunks, and a 1024² image at 4 spp, 16 777 216 lanes, replays bit-identically in 9 chunks.

5 Determinism

PRB’s replay must redraw the forward pass’s paths. irys makes this hold in four ways.

• Same pipelines. Forward and replay dispatch the same compute pipelines over the same bind groups; they differ only in a uniform that enables the store of .
• An integer RNG on both tiers. SplitMix64 is ported onto pairs of 32-bit words built from 16-bit limbs, so the GPU, which lacks 64-bit integers, produces exactly the CPU’s stream. A uniform keeps the top 24 of 53 bits, so both tiers walk the same paths up to f32 rounding of the geometry.
• No float atomics. The reduction over lanes is tang’s monoid reduction, split 256 ways with one workgroup per parameter and a fixed combining tree. Splat accumulation runs in bounce order, loss scoring uses a fixed stride and tree, chunks are summed in order, and the device prefix sum used by edge sampling has a fixed order.
• A walk that depends only on the ray. BVH traversal uses a fixed-size stack, visits the near child first by slab entry with a deterministic tie-break, and culls by the closest hit so far. The walk, and therefore the floating-point sequence, is a function of the ray alone, not of thread scheduling.

The result is that the replay reproduces the forward radiance bit for bit on all 65 536 lanes at depths 2 and 8, with NEE, and through the two-level instanced scenes and a garden scene of meshes and Gaussian splats; two runs give bitwise-identical loss and . The guarantee holds per device and driver. A different GPU or compiler may round differently, and then agrees with the CPU reference only to f32 precision.

6 Pose gradients through instanced CAD

Scenes. A scene is a two-level BVH: one bottom-level tree per CAD part or group, and a top-level tree over instances, each with an object-to-world transform. irys maps rays into an instance’s space without renormalising them, so hit distances compare across instances. A pose , such as joint angles, is a map from parameters to the transforms of the moving instances. The full robot in this paper has 2.4 million triangles, tessellated by vcad from its CAD source (a text format of parametric solids).

Primary edges. Moving an instance moves the curves across which the radiance arriving at the camera plane jumps: silhouettes, boundaries and creases between faces that are not coplanar. By Reynolds transport, a pixel’s derivative includes

over the projected edges , where is the edge’s outward normal and the normal speed of the material point under ⁵. We sample edges with density proportional to screen length times an importance (1 for silhouettes, a dihedral weight for creases, a small weight for back-facing edges), trace the two sides pixels apart with common random numbers, and require the inner ray to hit the edge’s own triangle or a coplanar neighbour. CAD tessellations contain many slivers thinner than , and without the coplanar rule the estimator loses most silhouette samples. The face that meets the eye more squarely owns a crease between two front faces, so the estimator counts its jump once. A low-resolution buffer of which instance each pixel sees lets the sampler demote edges it shows hidden, which cut variance 2× at equal samples. Sampling runs on the device, with a fixed-order prefix sum; it took the edge term on the arm from 273 ms to 15.3 ms per iteration. The device sampler is biased low on the full robot by a few percent, so the robot checks below use host samples (Section 8).

Interior term and attached frames. Away from discontinuities the estimator is smooth in at fixed random numbers, and we differentiate it by forward-mode dual numbers through a re-trace of the frame: ray–triangle intersection with the corners as duals, the shading frame, GGX sampling and evaluation, NEE and the MIS weights. The renderer builds a vertex’s tangent frame from its normal with a construction that flips sign when the normal crosses a plane, a discontinuity in at fixed random numbers that no edge term sees. We instead differentiate the estimator whose frame on a moving instance is rigidly attached, . At it is the rendered estimator exactly, and for every it is an equally valid estimator of the same image, so, wherever the integrand is smooth in and differentiation passes under the integral, its derivative is an unbiased estimate of from the paths the forward pass already drew; the edge terms handle the discontinuities where it does not. The duals reproduce the forward pass’s values to .

Secondary edges. When an instance moves, its shadow and its reflection move too. As in Li et al. ⁵, we sample the edges seen from each path vertex with a hierarchy over edges, here a per-tree edge BVH, under one-sample multiple importance over three densities: geometric, toward the sun’s cone, and toward the mirror lobe. The boundary term includes the velocity of the BSDF-sampled direction, obtained by inverting the sampler at the path’s random numbers and redrawing it in duals. Against the best setting of a warped-area estimator on the same configuration, at about the same cost per estimate, edge sampling had at least 21× lower variance on every parameter.

Camera. The same machinery differentiates a seven-parameter camera (translation, rotation, field of view), with the eye, axes and field of view as duals and every instance’s edges as candidates. This paper reports no check of the camera gradient.

7 Results

Unless stated, we check gradients against central differences (CD) of the same estimator with common random numbers; “frozen” means the random numbers are fixed so the CD is noise-free. Monte Carlo comparisons report a z-score, the difference over its standard error. Timings are release builds on the RTX 4090 machine above. The number in each table is read from data/measured.tsv, whose rows give the source test and commit.

7.1 Correctness of the material gradient

7.2 Variance and throughput

At equal samples per pixel, NEE with MIS lowers per-pixel variance 7.5× against BSDF sampling alone on a scene of two instanced boxes, with image means agreeing within noise. A gradient of the 96×96, 8 spp frame at depth 4 (73 728 paths, four wavelengths each) costs a median of 1.1 ms on the GPU without NEE (1.06–4.67 ms over four runs; the slowest run is the first, which includes pipeline compilation) and 3.88 ms with NEE (2.61–5.51 ms). The forward pass alone takes a median of 0.82 ms, about 89 million paths per second, so the replay and adjoint add 0.27 ms; NEE doubles the batch. The same gradient through the CPU tier’s compiled C costs a median of 555 ms (551–608 ms), measured in a separate test against the GPU gradient with the loss scored on the device, 1.04 ms (median of four runs), so the GPU tier is about 532× faster end to end on this machine. This compares irys with itself; it is not a comparison with another renderer.

7.3 Diffuse albedo, splats and emitters

The same replay differentiates two more kinds of parameter. A Lambertian BSDF is linear in its albedo, so a diffuse vertex adds to band of the albedo’s gradient, with the radiance before the vertex. A splat’s colour is linear too. Its opacity also scales everything behind it, so , where is the path throughput at the segment, the transmittance in front of splat , its colour, the emission composited behind it on the segment and the light beyond the segment. On the GPU the replay writes one record per splat crossed and the host sums the records in a fixed order, so there are still no float atomics. Area emitters with next-event estimation and multiple importance sampling run on both tiers, and their 16 gates pass on Vulkan.

7.4 Pose gradients

On a box against the sky, where the edge term is the whole derivative, edge sampling matches CD on two of its three pose parameters (z = 0.29 and 0.04). On the third it is 5.8% high (z = 3.41), in every run we have made; the test’s gate, four standard errors plus 3% of the reference, passes, but we have not found the cause. The device’s f32 edge term matches the CPU’s f64 to . On the shaded arm over a floor, the edge and interior terms alone miss CD by z = 3.56, −4.71, 3.71. The secondary term, mostly the arm’s shadow on the floor, is 13%, 34% and 23% of the three derivatives. With secondary edges the same check gives z = 0.23, −0.71, 0.39 at 64 seeds; at 256 seeds the check fails, because the secondary estimate is heavy-tailed (Section 8).

7.5 From pixels to CAD parameters

The edge sampler’s velocity is a projection of the material point’s motion, so replacing the pose’s velocity with a unit translation per axis yields for every mesh vertex. vcad’s differentiable seam ³ takes that cotangent on the nodes of a frozen tessellation and returns the derivative with respect to a CAD parameter in one reverse pass. The shared benchmark uses rose’s left shoulder-pitch output plate, a 3 mm rounded plate with a 49 mm bearing bore and twelve bolt holes, with the bore diameter. irys renders vcad’s frozen tessellation unchanged (4 942 nodes, 9 932 triangles). This is the mesh whose connectivity the kernel holds fixed while it differentiates node positions. The render is against the sky at 160×128, under a fixed smooth image loss.