Differentiating Through a Simulator

A simulator is a function iterated. One step maps state to state, $x_{t+1}=f(x_t,\theta)$, with parameters $\theta$ (controls, masses, gains) held fixed across the rollout, and a scalar loss $L(x_T)$ reads off the final state. To improve $\theta$ you want $\mathrm{d}L/\mathrm{d}\theta$. Running reverse-mode autodiff through the unrolled trajectory is not a separate trick: it is exactly the adjoint method, and it returns the full gradient in a single backward pass whose cost does not depend on how many parameters $\theta$ holds.

The unrolled trajectory

A rollout of $T$ steps is the composition $x_T=f(\cdot,\theta)\circ\cdots\circ f(x_0,\theta)$. The loss depends on $\theta$ through every step, so the chain rule threads a product of per-step Jacobians:

Forming those products left to right (forward mode) costs one pass per input direction, so $\theta\in\mathbb{R}^p$ needs $p$ passes. The adjoint reorders the sum to sweep once from the right.

Reverse mode is the adjoint

Define the adjoint $\lambda_t=\partial L/\partial x_t$, the sensitivity of the loss to the state at time $t$. Seed it at the end with the loss gradient, then propagate backward with the transposed state Jacobian, accumulating a parameter term at each step:

One backward traversal produces every entry of the gradient at once, so the cost is independent of $p$: a million-parameter controller and a one-parameter one cost the same backward sweep. That asymmetry, cheap gradient regardless of input dimension, is why reverse mode drives differentiable physics¹.

A scalar sim, by hand

Take the smallest nontrivial case: a scalar linear sim $x_{t+1}=a\,x_t+\theta$ with loss $L=(x_T-x^\star)^2$. The two partials are constant,

so the seed is $\lambda_T=2(x_T-x^\star)$, the recurrence is $\lambda_t=a\,\lambda_{t+1}$, and every accumulation term is just $\lambda_{t+1}$. Summing the geometric chain gives a closed form:

Now plug in $a=0.9$, $T=3$, $x_0=0$, $\theta=1$, $x^\star=2$. The forward rollout is

with $L=(2.71-2)^2=0.5041$. The backward sweep seeds $\lambda_3=2(2.71-2)=1.42$, then

and the parameter gradient is the running sum

The closed form agrees: $2(0.71)(1+0.9+0.81)=1.42\times2.71=3.8482$. And so does a central finite difference, which never sees the adjoint at all:

The demo runs the same recurrence for $T=8$: drag $a$ and $\theta$ to set the forward trajectory, then step the adjoint backward and watch each $\lambda_t$ get folded into $\mathrm{d}L/\mathrm{d}\theta$, ending at the finite-difference value.

The continuous limit and solver layers

As the step size $h\to0$ the backward recurrence becomes an ODE for $\lambda(t)$ integrated in reverse time, the continuous adjoint behind neural ordinary differential equations². And when a step is not an explicit update but an equation solve, say a contact quadratic program, you need not unroll the solver at all: the implicit function theorem differentiates its solution directly from the optimality conditions, giving an exact gradient through a fixed number of solver iterations³.

One caveat carries over from the contact-cliff view of differentiable simulation: when contact makes $f$ nonsmooth, $\partial f/\partial x_t$ jumps or spikes, and the adjoint faithfully propagates a gradient that is biased or high-variance. The machinery is exact for the model it is given; it cannot repair a discontinuous $f$.