Procedural UV Derivatives Evaluation in SORT Renderer

72-minute read

In this post, I want to talk about a topic I have been meaning to write about for quite a few years: texture UV derivatives evaluation in my offline renderer, SORT (Simple Open Source Ray Tracer), which I’ll use throughout this post. It normally is not a major challenge. However, my own custom shading language[1] specifically designed for my renderer adds quite a lot of complexity in it. It was previously named as Tiny Shading Language and renamed as SORT Shading Language (SSL for short) as I fused it into my rendererer project. This article may not be useful to everyone, but it should be relevant if you are interested in implementing an offline renderer and shading language yourself. If that describes you, I hope this post is useful, because this particular combination of topics is rarely discussed in one place.

Why Do We Care?

Derivatives are a powerful mathematical tool used widely across many scientific fields, such as deep learning, physics, biology, and many others. At first glance, they do not seem like a must-have for computer graphics. In practice, you can build a toy offline renderer that produces great images without them. SORT renderer used to lack derivative support entirely, yet it could still render respectable results.

As a matter of fact, derivatives do show up throughout computer graphics in many places. Signed distance fields (SDFs) use them to approximate surface normals. Geometry processing often relies on them for operations like smoothing and deformation. Smoothed particle hydrodynamics (SPH) methods simulate fluids with them[2]. Gradient-domain rendering is another interesting area built on the same foundation[3]. The Jacobians in ReSTIR’s shift mapping are essentially derivatives as well[4]. The list goes on. Still, the application every graphics programmer knows best is mipmapping and that is the main reason I spent so much effort figuring out how to evaluate derivatives in my offline renderer.

Mipmap Level Selected
Intel Sponza
Left, a render from the SORT renderer. Right, the mip level chosen by the method described in this post. Note that the right image shows, per primary view pixel, the average mip level across all textures sampled at that pixel. Asset courtesy of Intel.

Coming from a real time rendering background, I want to start with mipmapping in that context. Mipmapping is essential in real time rendering because it helps avoid texture aliasing. To hit real time frame rates, renderers can usually afford only a very limited number of samples per pixel, often just one. In fact, most modern game engines also use some form of upscaling[5, 6], which further lowers the effective sampling rate per pixel. With such a low sampling rate, any frequency above the Nyquist limit must be prefiltered to avoid artifacts. That is what mipmaps provide. They also help performance. Lower mip levels use less texture memory, are more likely to stay in the GPU cache, and can speed up texture sampling.

As useful as mipmaps are in real time rendering, their benefits are less decisive for offline renderers. Most Monte Carlo path tracers fight noise by increasing sample count, so with a much higher effective texture sampling rate, prefiltered textures matter less. Lower mip levels can still improve cache behavior, but that alone is not why production offline renderers invest in mipmapping.

There is another wrinkle. Mipmapping is not strictly unbiased. Even when the correct mip level is chosen from texture coordinate derivatives, filtering still introduces a small amount of error in theory. One source of bias is geometric, the sampling footprint induced by a path is rarely an axis aligned square in texture space, so approximating that footprint with filtering on a mipmap pyramid is not exact in general. Even when the footprint is axis aligned, if its boundary does not align exactly with a mip level’s texel grid, we usually blend the two nearest mip levels. That interpolation approximates detail between mip resolutions with a linear blend, which is typically biased. Setting that aside, consider an idealized setup, an orthographic camera, a viewport that matches a single quad (two triangles), and a texture with twice as many texels as screen pixels along each dimension. In that simple case, we should select the second most detailed mip level. To make the remaining bias concrete, suppose albedo is the product of two same-resolution textures, $T_A$ and $T_B$, with UVs $(u(x,y), v(x,y))$ at image location $(x,y)$. In theory, the integrated shading result is:

$$ \tag{1} F = \int f\!\left(x, y,\, T_A\!\left(u(x,y), v(x,y)\right)\, T_B\!\left(u(x,y), v(x,y)\right)\right) \mathrm{d}x\,\mathrm{d}y $$

Note that $f$ could simply be a function of $(x,y)$ in the equation above. I write it with explicit texture lookups because that form makes the comparison to the mipmapped case clearer.

With prefiltered mipmapped lookups at the second most detailed mip level, each texture is averaged over the pixel footprint $\mathcal{P}(x,y)$ in texture space before they are multiplied, giving

$$ \tag{2} F_{\mathrm{mip}} = \int f\!\left(x, y,\, \left(\int_{\mathcal{P}(x,y)} T_A(u,v)\,\mathrm{d}u\,\mathrm{d}v\right) \left(\int_{\mathcal{P}(x,y)} T_B(u,v)\,\mathrm{d}u\,\mathrm{d}v\right) \right) \mathrm{d}x\,\mathrm{d}y $$

It is apparent that $F_{\mathrm{mip}}$, the value obtained with mipmap filtering, is biased relative to $F$. Mipmapping is not unbiased, and unbiased estimation is critical in offline rendering, so production renderers need a compelling reason to use it anyway.

The main reason most commercial offline renderers adopt texture mipmaps is memory consumption. Instead of loading every texture at the start of a render, which is what most toy offline renderers do, a production renderer often touches only a tiny fraction of the data needed for the full image. During path tracing, when a mip level is requested, the renderer first checks whether it is already in the texture cache. If it is, it fetches the data and continues, much like a simple ray tracer. If not, the requesting thread is paused while an I/O thread loads the data from disk, and another worker can use the core in the meantime. When the load finishes, the original thread resumes. That may sound like overhead, but if physical cores stay busy and thread switching is cheap, which can be achieved through fibers[7], this cache can greatly reduce the texture memory needed to render a frame. In effect, the theoretical upper bound on memory usage for a shot is not only determined by asset size, but also output resolution[8, 9].

A texture cache in an offline renderer frees artists from being constrained by physical memory, it is possible only a fraction of the textures on disk is actually needed for a given shot. When the cache budget is exhausted, older data is evicted to make room for new requests. Mip levels that are never sampled are never loaded at all, which sharply cuts the renderer’s memory footprint. The cache budget still matters, you want enough headroom to avoid thrashing, but the system makes it practical to render scenes whose total texture data far exceeds the physical memory of the machine doing the work.

Texture coordinate derivatives determine which mip level to use for a given sample and therefore govern filtering quality. Because a texture cache loads mip levels on demand, the cache relies on accurate derivative estimates, incorrect derivatives cause over-filtering or under-filtering and can force unnecessary loads or evictions. Reliable derivative evaluation is thus a prerequisite for combining mipmapped textures with a paging cache, a requirement that motivates the work described in the remainder of this post.

Challenges in SORT Renderer

UV coordinate derivative evaluation is not rocket science, but it does require some calculus. There were several articles on evaluating UV derivatives in Ray Tracing Gems II[10, 11, 12]. Many open source renderers, including PBRT[13], also use partial derivatives for mip selection. On the surface, the problem does not look complicated.

To fully understand why this becomes a major challenge in SORT, it is important to understand the rationale behind how SORT processes materials in the first place.

Material System in SORT Renderer

Ever since I built a Blender plugin for SORT, I realized the renderer needed a shader graph based material system. It was an interesting challenge at the time, and such workflows are common in both film and games.

In game industry, engines typically gather all nodes and emit a single shader kernel (possibly split across several files) for the shader compiler. From the compiler’s perspective, there is no shader graph, only shaders. Though not every game engine uses shader graphs. For example, Naughty Dog’s in-house engine uses shader packages[14]. Offline renderers often take a different path. Open Shading Language accepts shader segments, and the compiler wires them together according to the host program. In effect, the compiler does the gathering, not the rendering engine. Other compilers may work differently. I do not have full visibility into all of them.

Initially, I used OSL in my renderer, then introduced my own shading language. A primary motivation was Apple Silicon support. At the time, OSL had no official Apple Silicon build, and it was unclear when that would arrive. By piggybacking on LLVM, I implemented a compiler that targets multiple CPU architectures (x64 and Arm64) and operating systems (Windows, Ubuntu, and macOS). Below is a brief overview of how SSL fits into SORT.

  • The Blender plugin exports shader segments as source code into a binary asset file based on material shader graph information.
  • It then spawns the renderer with an argument pointing to that asset. This step can run asynchronously.
  • During startup, SORT walks materials like other renderers, but instead of loading baked parameters from the asset, it loads shader code, similar to how real time engines compile HLSL. SSL compilation runs in a multithreaded environment. Each material stores the resulting JIT compiled function pointer.
  • During path tracing, whenever a material is needed, the renderer invokes that JIT compiled function, transferring control from C++ to compiled SSL code.
  • For texture sampling, SSL calls a C++ interface defined by SORT so the renderer handles sampling details.
    • This is intentional. Decoupling texture sampling from SSL leaves room for a texture cache system later, which would require pausing threads, outside SSL’s scope.
  • After any texture sampling, SSL resumes, completes its instructions, and returns a closure tree to the renderer. That tree becomes a BSDF, a stack of blended BxDF layers. Parameters such as albedo and roughness are evaluated entirely in SSL, which is its main purpose. The renderer then continues like any other path tracer.

As this workflow shows, SSL works in a similar fashion to OSL. It takes shader segments and wires them together. Its role is to move BxDF parameter evaluation from a hard coded fixed pipeline into artist programmable logic, not to own lighting or shading, despite the name “shading language”. Those parameters can depend on texture lookups, which in turn depend on UV coordinates, which ultimately depend on SSL’s global inputs. SSL global structure is an analogue to root signature concept in D3D12, it is a structure that passes data from the renderer (C++) into the shader (SSL). The exact formulas for texture coordinates are authored during content creation, not fixed at renderer compile time. That separation is the root cause of the UV derivative challenges discussed in this post. Because texture coordinates are not evaluated during renderer compilation, an analytical derivative solution cannot be implemented as a fixed function pipeline either.

Please check out my previous post if you are interested in learning more about the custom shading language that I used in SORT.

Why Existing Work don’t Apply

All of the aforementioned work, however, assumes a fixed function pipeline for texture coordinate generation, usually simply passing through the UVs stored in the mesh to texture sampling interface. That assumption does not hold in SORT renderer. With SSL, texture coordinates can be computed from essentially any shader input, such as position, normal, and so on. The exact mathematical formulation can be anything authored in the shader graph, which is not known at C++ compile time. There is no straightforward way to hard-code derivative evaluation for UV generation in the renderer.

I am not the first to run into this. The problem is largely solved in production. RenderMan’s RSL[15] can evaluate texture derivatives for mip selection, and OSL supports partial derivatives as well. Without RSL source code, I cannot simply integrate it, I would give up control over cross platform support. If I need a platform RSL does not support, I am back to looking for alternatives. I explained in an earlier post why I did not adopt OSL. I will not repeat that here. NVIDIA’s Slang[16] is another option. It has grown popular, supports derivatives natively, and can run on the CPU with the right setup[17]. But Slang has a blocker for SORT, texture sampling. The point of mipmapped textures is eventually a texture cache system in the renderer[9]. I need to pause a thread inside a texture sampling call from C++, not from the shading language. With Slang, sampling lives in the language, which is a deal breaker in my case. None of these existing paths fit cleanly, so the remaining option is to implement derivatives in SSL myself.

Along the way I picked up backpropagation from deep learning[18], a way to compute gradient given a few outputs, many inputs, and a multi layer network in between. Backpropagation is essentially reverse mode automatic differentiation. That led me to automatic differentiation more broadly, and I studied Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation, Second Edition[19]. The book is quite academic but insightful. I did not finish reading the whole book, but what I read was enough to shape how I might implement partial derivatives in my SSL.

Different ways to Evaluate Derivatives

Given my goal is to support derivative evaluation in SSL, the next step is to clarify how this can be achieved. Before looking for a concrete engineering solution, it helps to review the common approaches to derivative evaluation[20]. To keep the comparison fair, we will use the same problem throughout this section and compute $\partial f / \partial x_1$ with each method.

$$ \tag{3} f(x_1, x_2) = \Bigl( \sin\!\Bigl(\frac{x_1 x_2}{x_1 + 1}\Bigr) + \ln\!\Bigl(1 + \frac{x_1 x_2}{x_1 + 1}\Bigr) - e^{x_2} \Bigr) \Bigl( \frac{x_1 x_2}{x_1 + 1} - \tanh(x_2) \Bigr) $$

Manual Differentiation

Let’s start with the one we already know, manual differentiation. In essence, you examine the underlying math and derive a closed-form expression for the derivatives of the original evaluation.

Applying the chain rule to equation 3 gives

$$ \tag{4} \frac{\partial f}{\partial x_1} = \frac{x_2}{(x_1 + 1)^2} \left[ \bigl(\sin u + \ln(1+u) - e^{x_2}\bigr) + \bigl(u - \tanh x_2\bigr)\Bigl(\cos u + \frac{1}{1+u}\Bigr) \right] $$

To keep the result compact, I introduce an intermediate variable $u$.

$$ \tag{5} u = \frac{x_1 x_2}{x_1 + 1} $$

Much existing work that assumes a fixed texture coordinate pipeline can adopt this approach. But is it feasible in a shader graph context? Two workflows are conceivable in theory.

  • One option is to ask artists to do the math themselves, implementing derivatives of the texture coordinates inside the shader graph and wiring them into the texture sampling interface. That puts the burden on artists, is inefficient, and is error prone, because the derivative subgraph can easily drift out of sync with the forward evaluation.
  • Another option is to embed derivative code in each shader segment, alongside the forward implementation of every graph node. This is a one-time cost for the engine author and removes the burden from artists, but hand authoring those rules is tedious. Worse, each node’s segment is compiled in isolation before the graph is connected, so it must conservatively propagate derivatives for every output. At that stage, the compiler cannot know which variables will eventually feed a texture sampler. That worst-case preparation is usually wasted. A compiler might strip unused derivative code later, but relying on dead code elimination to fix an overly pessimistic design is not an elegant solution.

In practice, neither path is something we would expect in a real renderer with shader graphs.

Numerical Differentiation

Numerical differentiation follows directly from the mathematical definition of derivatives. Starting from

$$ \tag{6} f'(x) = \lim_{\delta \rightarrow 0} \dfrac{f(x+\delta) - f(x)}{\delta} $$

in plain terms, if we shift the input by a very small amount, we can evaluate the function at both the shifted and original positions. Dividing their difference by the shift size yields the derivative as $\delta$ approaches zero.

Applied to equation 3, the partial derivative is

$$ \tag{7} \frac{\partial f}{\partial x_1} = \lim_{\delta \rightarrow 0} \dfrac{f(x_1 + \delta, x_2) - f(x_1, x_2)}{\delta} $$

As the name suggests, numerical differentiation replaces the limit with a finite (but small) value of $\delta$, giving the finite difference approximation

$$ \tag{8} \frac{\partial f}{\partial x_1} \approx \dfrac{f(x_1 + \delta, x_2) - f(x_1, x_2)}{\delta} $$

This approach is simple and does not require knowing a derivative rule in advance. However, it has important trade-offs.

  • First, evaluating $N$ partial derivatives, one per input variable, requires $O(N)$ function evaluations. That cost is mild in this post, where we only need two partial derivatives per texture coordinate (with respect to screen space $x$ and $y$). It becomes prohibitive in other settings, such as deep learning, where $N$ can reach millions.
  • Second, choosing $\delta$ is not always straightforward. If $\delta$ is too large, the approximation deviates from the true derivative. This error is called truncation error. If $\delta$ is too small, floating-point rounding error dominates and the estimate becomes unstable.

Despite these drawbacks, this method is widely used in practice. In real time rendering, GPUs effectively rely on numerical differentiation for implicit derivatives. Shader programs are evaluated across many threads in parallel, and pixel shaders operate over 2x2 quads that exist in a warp/wavefront. Within a warp/wavefront, threads are synchronized, which enables hardware to obtain numerical estimates for differentiable expressions in the shader kernels literally for free, including texture coordinates. When texture sampling occurs without an explicit mip selection, the GPU numerically estimates the texture coordinate derivatives and feeds them to the texture sampling unit for mipmap level selection. In most cases, this happens transparently to the programmer.

Symbolic Differentiation

Symbolic differentiation is closely related to manual differentiation. Given an explicit formula for a texture coordinate, the algorithm derives a closed-form expression for its derivatives automatically, rather than by hand. A correct implementation applied to equation 3 produces a derivative expression equivalent to equation 4, and evaluates to the same value.

This can happen at shader compile time, which avoids extra work for artists. The compiler evaluates the generated expression at runtime with the same inputs as the forward pass. For a shader graph workflow, this may look appealing because derivatives can be generated under the hood without artist intervention.

However, symbolic differentiation does not handle conditionals, loops, or recursion cleanly. Recursion is unlikely to matter for UV derivative evaluation, but conditionals and loops are common in shader graphs. We cannot discard them this early when designing the solution.

What makes it a lot less compelling is that symbolic differentiation suffers from expression swell, just like manual differentiation. The compact form in equation 4 relies on introducing $u$. If we feed equation 3 to a symbolic differentiator literally, the subexpression $\frac{x_1 x_2}{x_1 + 1}$ appears three times and the output grows accordingly before any simplification pass. Even when a single node has a modest derivative, composed nonlinearities can be worse. Soft ReLU, a standard activation in deep learning[18], is a good example of it.

$$ \tag{9} f(x) = \log\!\left(1 + e^{wx+b}\right) $$

On its own, its derivative is modest. But once such terms are composed, the closed-form result can grow quickly. For example, if one activation feeds into another.

$$ \tag{10} f(x) = \log\!\left(1 + e^{b_2 + w_2 \log\left(1 + e^{b_1 + w_1 x}\right)}\right) $$

then

$$ \tag{11} f'(x) = \frac{w_1 w_2 \, e^{b_1 + w_1 x} \, e^{b_2 + w_2 \log\left[1 + e^{b_1 + w_1 x}\right]}}{\left(1 + e^{b_1 + w_1 x}\right)\left(1 + e^{b_2 + w_2 \log\left[1 + e^{b_1 + w_1 x}\right]}\right)} $$

In the worst case, symbolic differentiation can produce expressions far larger than the original, sometimes exponentially so, which directly affects derivative evaluation cost in SSL.

Automatic Differentiation

Automatic differentiation, sometimes also called algorithmic or computational differentiation, is another approach. Unlike symbolic differentiation, it operates on a program and produces numeric derivative values, not expanded formulas.

There are two main variants, forward mode and reverse mode.

  • Forward mode propagates derivative values alongside the forward pass. Its cost scales with the number of input variables you differentiate with respect to. That makes it a natural fit here. A shader graph may sample many textures and produce many intermediate values, but for mip selection we only need texture coordinate derivatives with respect to screen space $x$ and $y$. Those are the only two input directions that matter.
  • Reverse mode propagates adjoints backward through the program. Its cost scales with the number of outputs you differentiate. It excels when there are many inputs and few outputs, as in neural network training. In our setting the number of texture samples can still be large, making it a bad fit for our problem. And reverse mode runs a backward pass through the graph after the forward pass, which does not match the natural flow of shader evaluation as cleanly. It is especially awkward when a texture sample appears midway through the shader. We need derivatives at that point before the forward pass continues, which can force the runtime to alternate between forward and reverse passes within a single execution. And this can happen many times in a shader kernel.

For brevity, this post focuses on forward mode automatic differentiation. Reverse mode is a poor fit for the UV derivative problem described here.

To see how forward mode differentiation works, we first break the expression 3 into a sequence of elementary operations. Most programmers would not evaluate it in a single statement anyway. They use intermediate variables so shared subexpressions are computed once.

 1void eval(float x1, float x2, out float o)
 2{
 3    const float v0 = x1 + 1.f;
 4    const float v1 = x1 * x2;
 5    const float v2 = v1 / v0;
 6    const float v3 = sin(v2);
 7    const float v4 = log(1.f + v2);
 8    const float v5 = exp(x2);
 9    const float v6 = tanh(x2);
10    const float v7 = v3 + v4;
11    const float v8 = v7 - v5;
12    const float v9 = v2 - v6;
13    const float v10 = v8 * v9;
14
15    o = v10;
16}

The code above is written in SSL. Most authors would not write it exactly this way, especially with these variable names, but it is a correct decomposition of equation 3. I use this style to make the forward mode walkthrough easier to follow. Even when the source uses compound expressions, automatic differentiation still applies. The compiler lowers each line to primitive operations in the generated code.

Below is the computation graph for the implementation above.

The computation graph makes the dependencies between intermediate variables explicit. Another common representation is an evaluation trace, a sequence of variable definitions evaluated in order. Below is the trace at $(x_1, x_2) = (1.5, 0.5)$.

Variable Assignment Value
$x_1$ 1.5000
$x_2$ 0.5000
$v_0$ $x_1 + 1$ 2.5000
$v_1$ $x_1 x_2$ 0.7500
$v_2$ $v_1 / v_0$ 0.3000
$v_3$ $\sin(v_2)$ 0.2955
$v_4$ $\log(1 + v_2)$ 0.2624
$v_5$ $e^{x_2}$ 1.6487
$v_6$ $\tanh(x_2)$ 0.4621
$v_7$ $v_3 + v_4$ 0.5579
$v_8$ $v_7 - v_5$ −1.0908
$v_9$ $v_2 - v_6$ −0.1621
$v_{10}$ $v_8 v_9$ 0.1768
$o$ $v_{10}$ 0.1768

A key observation is that although every variable, intermediate or final, ultimately depends on the inputs $(x_1, x_2)$, each variable’s value depends only on its immediate operands during evaluation. For example, once $v_3$ and $v_4$ are known, $v_7 = v_3 + v_4$ follows from those values alone, we need not revisit the original inputs. This locality is exactly what makes forward mode automatic differentiation work. It holds for values, and it holds for derivatives as well. By the chain rule, $\partial v_7 / \partial x_1 = \partial v_3 / \partial x_1 + \partial v_4 / \partial x_1$, which depends on the derivatives of $v_3$ and $v_4$ only.

To evaluate $\partial o / \partial x_1$, we need $\partial v_{10} / \partial x_1$. From the product $v_{10} = v_8 v_9$, that requires $v_8$, $v_9$, $\partial v_8 / \partial x_1$, and $\partial v_9 / \partial x_1$. Note that $v_8$ and $v_9$ are already computed during forward evaluation, the extra work is computing the two partial derivatives. This decomposes recursively, $\partial v_8 / \partial x_1$ needs $\partial v_5 / \partial x_1$ and $\partial v_7 / \partial x_1$ (along with $v_5$ and $v_7$), and similarly for $\partial v_9 / \partial x_1$. In other words, we reason about the graph backward, from the output toward the inputs, until we reach seeds such as $\partial x_1 / \partial x_1 = 1$ and $\partial x_2 / \partial x_1 = 0$.

The reasoning runs backward, but the implementation runs forward. For each intermediate variable in the value pass, we allocate a companion variable that stores its derivative with respect to the input of interest. Derivatives are computed in the same order as values. In practice, the derivative update for a variable usually sits right next to the code that evaluates that variable, as if derivatives were part of the same forward pass. That is why this is called forward mode automatic differentiation. For notational convenience, let’s denote $\dot{v}_i = \partial v_i / \partial x_1$ for each intermediate variable $v_i$’s derivative. Below is the same evaluation trace at $(x_1, x_2) = (1.5, 0.5)$, extended on the right with the corresponding derivative variables, their update rules.

Variable Assignment Value Derivative Assignment Value
$x_1$ 1.5000 $\dot{x}_1$ 1.0000
$x_2$ 0.5000 $\dot{x}_2$ 0.0000
$v_0$ $x_1 + 1$ 2.5000 $\dot{v}_0$ $\dot{x}_1$ 1.0000
$v_1$ $x_1 x_2$ 0.7500 $\dot{v}_1$ $\dot{x}_1 x_2 + x_1 \dot{x}_2$ 0.5000
$v_2$ $v_1 / v_0$ 0.3000 $\dot{v}_2$ $(\dot{v}_1 v_0 - v_1 \dot{v}_0) / v_0^2$ 0.0800
$v_3$ $\sin(v_2)$ 0.2955 $\dot{v}_3$ $\cos(v_2),\dot{v}_2$ 0.0764
$v_4$ $\log(1 + v_2)$ 0.2624 $\dot{v}_4$ $\dot{v}_2 / (1 + v_2)$ 0.0615
$v_5$ $e^{x_2}$ 1.6487 $\dot{v}_5$ $v_5 \dot{x}_2$ 0.0000
$v_6$ $\tanh(x_2)$ 0.4621 $\dot{v}_6$ $(1 - v_6^2),\dot{x}_2$ 0.0000
$v_7$ $v_3 + v_4$ 0.5579 $\dot{v}_7$ $\dot{v}_3 + \dot{v}_4$ 0.1380
$v_8$ $v_7 - v_5$ −1.0908 $\dot{v}_8$ $\dot{v}_7 - \dot{v}_5$ 0.1380
$v_9$ $v_2 - v_6$ −0.1621 $\dot{v}_9$ $\dot{v}_2 - \dot{v}_6$ 0.0800
$v_{10}$ $v_8 v_9$ 0.1768 $\dot{v}_{10}$ $\dot{v}_8 v_9 + v_8 \dot{v}_9$ −0.1096
$o$ $v_{10}$ 0.1768 $\dot{o}$ $\dot{v}_{10}$ −0.1096

So $\dot{o} = \partial o / \partial x_1 \approx -0.1096$ at this point.

The update rules above map directly onto code. Below is an interleaved version of eval that computes each $v_i$ and its companion $\dot{v}_i = \partial v_i / \partial x_1$ in the same forward pass, using the seeds $\dot{x}_1 = 1$ and $\dot{x}_2 = 0$.

 1void eval_updated(float x1, float x2, out float o, out float dot_o)
 2{
 3    const float dot_x1 = 1.f;
 4    const float dot_x2 = 0.f;
 5
 6    const float v0 = x1 + 1.f;
 7    const float dot_v0 = dot_x1;
 8
 9    const float v1 = x1 * x2;
10    const float dot_v1 = dot_x1 * x2 + x1 * dot_x2;
11
12    const float v2 = v1 / v0;
13    const float dot_v2 = (dot_v1 * v0 - v1 * dot_v0) / (v0 * v0);
14
15    const float v3 = sin(v2);
16    const float dot_v3 = cos(v2) * dot_v2;
17
18    const float v4 = log(1.f + v2);
19    const float dot_v4 = dot_v2 / (1.f + v2);
20
21    const float v5 = exp(x2);
22    const float dot_v5 = v5 * dot_x2;
23
24    const float v6 = tanh(x2);
25    const float dot_v6 = (1.f - v6 * v6) * dot_x2;
26
27    const float v7 = v3 + v4;
28    const float dot_v7 = dot_v3 + dot_v4;
29
30    const float v8 = v7 - v5;
31    const float dot_v8 = dot_v7 - dot_v5;
32
33    const float v9 = v2 - v6;
34    const float dot_v9 = dot_v2 - dot_v6;
35
36    const float v10 = v8 * v9;
37    const float dot_v10 = dot_v8 * v9 + v8 * dot_v9;
38
39    o = v10;
40    dot_o = dot_v10;
41}

Each dot_v line is the code form of the corresponding derivative assignment in the table. The value lines are unchanged from the original eval. Of course, this update function only evaluates derivative with regard to $x_1$, if derivatives with regard to other inputs are needed, we can insert more instructions to make it happen.

Hand-writing the interleaved program is workable for a toy eval, but it does not scale to a full shading language with control flow and large graphs. That is where automatic differentiation comes in, the compiler emits the dot_v updates from the value code.

From Pencil and Paper to the Compiler

Now that we know the theoretical solutions to the derivative problem, it is time to get our hands dirty implementing them in the compiler. The goal is straightforward. For a practical implementation in SSL, the compiler should provide derivatives with respect to screen space $x$ and $y$ whenever a texture sample needs them.

Let’s use equation 3 as a concrete example. In practice, a shader author would write it the way they would in any other language, reuse the shared ratio once, pick local names that make sense, and move on. Below is what that might look like, the same math as the eval function above, but written as two compound expressions that read more naturally.

1float eval_practical(float x1, float x2)
2{
3    const float u = x1 * x2 / (x1 + 1.f);
4    return (sin(u) + log(1.f + u) - exp(x2)) * (u - tanh(x2));
5}

Imagine feeding that function into a texture coordinate inside an SSL shader entry.

 1texture2d g_albedo;
 2shader shader_entry(out closure output)
 3{
 4    const float3 fake_normal = vector(0.0f, 1.0f, 0.0f);
 5    const float3 global_input_pos = global_value<position>;
 6    const float x = global_input_pos.x;
 7    const float y = global_input_pos.y;
 8    const float z = global_input_pos.z;
 9    const float v = eval_practical(y, z);
10    const color basecolor = texture2d_sample<g_albedo>(x, v);
11    output = make_closure<Lambert>(basecolor, fake_normal);
12}

Clearly, this is not a sensible way to compute texture coordinates in production. It is only an example of the kind of procedural math SSL must be able to differentiate. We can ignore the odd UV mapping that results. If anything, the contrived coordinate arithmetic is representative of how arbitrary texture coordinate formulas can be in real shader source.

A few constructs in the snippet above are specific to SSL and may look unfamiliar if you are used to OSL or RSL. They are not meant as a general template for the language.

  • make_closure<Lambert> allocates a node in the closure tree that the renderer evaluates later.
  • global_value<position> reads a field from the SSL global block, the CPU-side data structure filled in before each shader execution.
  • texture2d_sample<g_albedo> samples the texture bound to the global handle g_albedo. That call crosses from JIT’d SSL into SORT’s C++ texture path, which is how the renderer can own sampling while the shader still runs as a single kernel.

Explaining SSL’s full language design is outside the scope of this post. For background, see my earlier blog post.

The compiler would derive the derivatives of eval_practical from the expressions inside it and pass them into texture2d_sample, so the renderer can select the correct mip level. That is the goal of this post, every texture sample should receive not only the texture coordinates, but also their partial derivatives with respect to screen space $x$ and $y$.

With that goal in mind, the next step is to anticipate the engineering problems that come up when theory meets a real compiler. What immediately jumped out to my mind was the following.

  • How do I deal with derivatives that cross function boundaries?
  • What if the user passes constant data that does not depend on SSL global inputs for a texture coordinate?
  • Where do I store derivative data as the shader executes?
  • For which variables do I need to track derivatives? Or do I track derivatives for all variables and rely on LLVM to eliminate dead code?
  • How do I handle conditionals, loops, or even recursive calls?
  • Would the extra instructions hurt performance?

The list is not exhaustive, but it captures the questions I had to answer before committing to an implementation. Bringing derivatives into SSL turned out to be a substantial project. It was not something I could solve by asking an AI to do in one step. It took quite a long time before I landed on a workable approach.

Implementing Derivative Lanes with SIMD

It took me about half a year to build this solution, a detour in hindsight. I had not yet studied automatic differentiation, so given my real time rendering background, it is not surprising that I tried this path first.

As explained earlier, GPU threads run in synchronized groups, warps or wavefronts. The hardware exploits that synchronization to estimate derivatives of arbitrary expressions at negligible cost.

Naturally, I looked for a CPU analogue, extra SIMD lanes in SSL, not to speed up forward evaluation, but to carry helper values so derivatives could be approximated numerically within the same execution, the way a GPU quad does. At first this looked attractive, virtually every CPU supports SSE2[21], and the extra lanes seemed like they would cost almost nothing.

I started by duplicating all data in SSL, global inputs, locals, structure members, array entries, and so on. With some initial success, I could approximate derivatives numerically. It looked promising until more and more fragile design choices surfaced.

  • The approach attaches derivative lanes to all variables, but only a tiny subset ever feeds a texture sampler. That is a large memory waste. On a GPU the comparison is different, helper lanes usually correspond to neighboring pixels, each with a real forward value. Helper lanes may be wasted on sub pixel-sized triangles, which is one reason we should avoid them. But in SSL the second and third lanes exist only for derivatives. Even in the case most variables need derivatives, it still wastes at least 25% of memory comsuption since the forth lane is totally useless.
  • Divergence was painful. Control-flow divergence happens when lanes in a SIMD group disagree on branches. A GPU typically serializes the taken paths, running each with lane masks until the warp reconverges. Data divergence happens when lanes access different addresses, for example, when each lane indexes an array differently, forcing gathers or scatters instead of one uniform load. In either case the hardware keeps every lane live at extra cost, which assumes each lane matters equally. In SSL, re-evaluating divergent paths for derivative lanes adds little value, but skipping them can produce confusing, hard to debug behavior.
  • I eventually stopped evaluating the second and third lanes altogether. That forced a new qualifier, primary, and only primary-qualified values may drive array indices, conditionals, or anything else that can diverge. The scheme worked to a degree, but it burdened shader authoring, every assignment had to preserve qualifier rules, and mixing qualifiers required careful typing (for example, a primary result could only be built from primary operands). The extra ceremony was another reason I soured on the approach.
  • DreamWorks’ Moonray renderer features a vectorized path tracer[22] that splits work into small kernels so each can be SIMD-parallelized across rays. I experimented with that direction too. Repurposing SIMD inside SSL for derivative lanes conflicts with vectorizing SSL itself, if shading is batched for throughput, the lanes are already spoken for.

With all of these problems in view, I abandoned the SIMD implementation. It was an unfortunate decision. However, it was a poor fit for this problem, which pushed me to look for alternatives that could sidestep these limits.

Implementing Forward Mode AD in SSL

After studying forward mode automatic differentiation, it became clear that this was the right tool for SSL. The overall implementation strategy is straightforward. The AST from the existing scanner and parser stays as is, and what changes is codegen. Instead of emitting only the primary evaluation path, the compiler also emits instructions that maintain derivative shadows ($\partial/\partial x$ and $\partial/\partial y$) wherever demand requires them. Unless stated otherwise, derivatives in the rest of this section mean screen space partials with respect to pixel $x$ and $y$.

Expanding SSL Global Structure

As mentioned earlier, the SSL global is the block the host fills before SSL shader kernel runs. Authors declare its layout in C++ with a small set of macros. In SORT, the hit-point payload looks like this.

1BEGIN_SSLGLOBAL_STRUCT(SSLHitGlobal)
2    SSLGLOBAL_PARAMETER(SSL_float3, uvw)
3    SSLGLOBAL_PARAMETER(SSL_float3, position)
4    SSLGLOBAL_PARAMETER(SSL_float3, normal)
5    SSLGLOBAL_PARAMETER(SSL_float3, gnormal)
6    SSLGLOBAL_PARAMETER(SSL_float3, I)
7    SSLGLOBAL_PARAMETER(SSL_float3, tangent)
8END_SSLGLOBAL_STRUCT()

This macro-based layout is one of several SSL cleanups I made before tackling automatic differentiation. Compared with my previous implementation[1], it reads much more clearly and follows the same pattern used in Unreal Engine 5.

Before automatic differentiation, those macros expand to a plain C++ struct.

1struct SSLHitGlobal{
2    float3 uvw;
3    float3 position;
4    float3 normal;
5    float3 gnormal;
6    float3 I;
7    float3 tangent;
8}

The macros let the C++ compiler gather every field at compile time. The SSL compiler uses that list to build a matching LLVM struct before it compiles shader code.

With that background in place, here is why the global block must grow. For convenience, I’ll repeat the earlier example here.

 1float eval_practical(float x1, float x2)
 2{
 3    const float u = x1 * x2 / (x1 + 1.f);
 4    return (sin(u) + log(1.f + u) - exp(x2)) * (u - tanh(x2));
 5}
 6
 7texture2d g_albedo;
 8shader shader_entry(out closure output)
 9{
10    const float3 fake_normal = vector(0.0f, 1.0f, 0.0f);
11    const float3 global_input_pos = global_value<position>;
12    const float x = global_input_pos.x;
13    const float y = global_input_pos.y;
14    const float z = global_input_pos.z;
15    const float v = eval_practical(y, z);
16    const color basecolor = texture2d_sample<g_albedo>(x, v);
17    output = make_closure<Lambert>(basecolor, fake_normal);
18}

Below is the computational graph for it.