Choosing the Right Derivative: A Practical

Comparison of Julia’s Automatic Differentiation

Strategies
Jason liu
Master Student Software Engineering
University of Antwerp
Antwerp, Belgium
[Link]@[Link]

Abstract—Automatic differentiation (AD) is the core mecha­ known derivative. By applying the chain rule at every step, the
nism behind differentiable programming, enabling gradients to derivative of the whole program is built up incrementally [1].
be computed through arbitrary programs. This paper compares Forward mode AD attaches a derivative component to every
four AD strategies in Julia, a self-implemented forward mode
value in the computation, requiring n passes for a function
using dual numbers, a self-implemented reverse mode using
a tape, [Link], and [Link], across three dimensions: cor­ with n inputs. Reverse mode AD records the computation
rectness, runtime performance, and implementation complexity. during a forward pass and recovers all n gradients in a single
Experiments are conducted on a parameter fitting problem using backward sweep, making it the preferred choice when n is
the Lotka-Volterra predator-prey ODE system. Results show that large and the output is a single scalar. [1], [2]
Enzyme is the fastest by a large margin due to LLVM-level The Lotka-Volterra predator-prey model [3] serves as the
differentiation, forward mode is the simplest to implement, and experimental vehicle, a system of two coupled ODEs describ­
Zygote suffers significant overhead on scalar simulation code due
to heap allocation per operation.
ing population dynamics of prey and predators. We frame it
Index Terms—Julia, Automatic differentiation, AD, Forward, as a parameter fitting problem: given observed trajectories,
Reverse, Enzyme, Zygote find the parameters minimizing a least-squares loss, and use
AD to compute the gradient of that loss.
I. Introduction
III. Experiment & Results
Differentiable programming is a paradigm in which pro­
grams are written such that derivatives can be computed Our experiment is divided into four sections. First, we ver­
through them automatically. This enables gradient-based ify that all differentiation techniques produce correct gradients
optimization, the backbone of machine learning, scientific by comparing them against finite differences (correctness).
parameter estimation, and numerical simulation, to be applied Next, we measure the time each technique takes to compute
to arbitrary code rather than only to hand-derived formulas. those gradients (timing). Then, we scale the problem by
Automatic differentiation (AD) is the computational tech­ increasing the number of parameters to observe how each
nique that makes this possible. Unlike symbolic differentia­ method behaves under a heavier workload.
tion, which manipulates mathematical expressions, or numer­ The experiment is based on the Lotka-Volterra predator-
ical differentiation, which approximates derivatives by finite prey model [3], a system of two coupled ordinary differential
differences, AD computes exact derivatives by applying equations that describes the population dynamics of two inter­
the chain rule systematically through every operation in a acting species, prey (rabbits) and predators (foxes). Rather
program. There are two fundamental modes: forward mode, than solving the equations directly, we frame the problem as
which propagates derivatives alongside values in a single parameter fitting: given observed population trajectories gen­
forward pass, and reverse mode, which records the compu­ erated from known true parameters, we compute the gradient
tation and recovers all gradients in a single backward pass. [1] of a least-squares loss with respect to all model parameters.
In practice, multiple AD systems exist with different design This gives us a realistic many-inputs-to-one-output problem,
philosophies and performance characteristics. This paper eval­ which is precisely the setting where the difference between
uates four strategies available in Julia and examines how they forward and reverse mode automatic differentiation becomes
differ in correctness, speed, and the burden they place on the most apparent.
programmer. A. Methods

II. Background Here we will describe the four methods that we choose to
benchmark, where two of them are custom made.
Automatic differentiation exploits the fact that any program a) Forward mode (self implemented):
is ultimately composed of primitive operations, each with a Forward mode automatic differentiation is implemented
using dual numbers. A dual number carries two components
simultaneously: the actual value of the computation and its
derivative with respect to one input. By overloading all arith­
metic operators to propagate both components according to
the chain rule, any existing function can be differentiated
without modification.
struct Dual
value::Float64
deriv::Float64
end
To compute the gradient of a function with n inputs, forward
mode must perform n separate passes, one for each input
parameter, seeding that input’s derivative component to 1.0
and all others to 0.0. This is the fundamental cost of forward
mode: the number of passes scales linearly with the number
of inputs. [1], [4], [5]
Note: to keep the implementation even simpler, we only
overloaded the basic things and some other things we needed.
To make this fully functional, you need to overload way more
than what we did. Fig. 1. Zygote workflow
b) Reverse mode (self implemented):
Phase 1 — Forward pass: You run the computation nor­ At its core, Zygote uses reverse-mode automatic differenti­
mally, but secretly record every single operation onto a “tape” ation, which is ideal for deep learning scenarios with many
as you go. You’re not computing any derivatives yet, you’re input parameters and a single scalar output like a loss function.
just logging what happened and saving the input values of [7], [8]
each operation, because you’ll need those values later. d) Enzyme:
Phase 2 — Backward pass: You seed the output with a [Link] is a high-performance automatic differentiation
gradient of 1.0 (meaning “the output is 100% sensitive to (AD) system for Julia that operates by binding to the
itself”), then walk the tape backwards. At each recorded LLVM-based Enzyme differentiator, allowing it to compute
operation you apply the local derivative rule and push the derivatives of Julia functions by analyzing and transforming
gradient back to the inputs of that operation. By the time LLVM intermediate representation (IR) rather than source
you’ve walked the entire tape, every input has accumulated code. [9], [10]
its [Link] crucial difference from forward mode is that
one backward pass gives you all gradients simultaneously.
You don’t re-run the simulation for each parameter, you run
it once, record it, then read all gradients off in one backward
sweep. [1], [6]
Note: again we overload only what is needed to keep the
implementation at minimum. Fig. 2. Zygote workflow
c) Zygote:
[Link] operates using source-to-source transformation, Enzyme differs from Zygote, where we can specify ourself
where it parses existing Julia code and generates new Julia what differentiation we want to apply. We can both choose
code specifically designed to compute gradients, rather than forward and reverse. In our experiment, we will make use of
building a computational graph or relying on operator over­ reverse because we want to see how the best differentiation
loading. compares to others.
B. Experiments
As mentioned before we will run three different types of
experiments.
a) Setup:
First let us discuss a bit about the minor setup we do for
the benchmarking tests.
We choose to opt to Euler integration, mainly because this
is the simplest way to solve ODE. [11], [12]
const DT = 0.01
const STEPS = 400

const TRUE_α = 1.0

const TRUE_β = 0.1
...
 Where DT is what we use for the Euler method, how TABLE II
long the distance is we travel for every step. Next, we have Gradient computation time for the 6-parameter Lotka-Volterra
problem.
our amount of steps we will take before we conclude our
gradients. The TRUE_* are the ground truth values we just Method Median Min Std
assume that we have observed.
Next, because we don’t have real observed data, we just Forward (Dual) 15.3 μs 13.6 μs 15.9 μs
generate this and state that this is our real observed data:
const OBS_X, OBS_Y = Reverse (Tape) 1.767 ms 393.4 μs 3.797 ms
generate_observations(TRUE_α, ...)
const PARAMS = [1.1, 0.09, 1.4, 0.08, Zygote 5.364 ms 1.922 ms 6.685 ms
10.5, 4.8]
This method of experimenting is what we call synthetic Enzyme 5.5 μs 4.2 μs 7.36 μs
experiments and will suffice for this small benchmarking.
The PARAMS will represent our permuted data, and is what Table 2 reports the median gradient computation time
we will use to find (eventually) our ground truth. for each method on the standard 6-parameter Lotka-Volterra
And at last we need to make the simulation functions for problem, measured using BenchmarkTools over 200 runs to
each method (some share the same function) so we can finally avoid JIT compilation noise. [13]
start experimenting. There are three versions of the same Enzyme is the fastest method at 5.5 μs, nearly three times
simulation for three different API’s. faster than forward mode and roughly 1000 times faster than
b) Correctness: Zygote. This is expected, Enzyme operates at the LLVM IR
Before measuring performance, we verify that all four level, below Julia’s object system entirely, meaning interme­
methods produce correct gradients. Each method is evaluated diate values live in CPU registers rather than on the heap. No
at the perturbed parameter vector and compared against a memory allocation occurs during the backward pass. [8], [10]
central finite difference baseline computed as: Forward mode at 15.3 μs is the second fastest despite being
𝜕𝐿 𝐿(𝜃𝑖 + 𝜀) − 𝐿(𝜃𝑖 − 𝜀) a simpler implementation. For 6 parameters it runs 6 separate
≈ (1)
𝜕𝜃𝑖 2𝜀 simulation passes, but each pass is lightweight, it carries only
with epsilon = 10−5 . We report both the maximum absolute two floating point values per operation and allocates nothing
error and maximum relative error across all six parameters. extra beyond the Dual numbers themselves.
The manual reverse mode tape is significantly slower at
TABLE I 1.767 ms. While it correctly computes all gradients in a
Gradient correctness across all four AD methods compared to single backward pass, the tape infrastructure introduces over­
central finite differences. head: every operation during the forward pass allocates a
Method Max Absolute Error Max Relative Error
TapeEntry struct onto the heap, and the backward pass
must walk through thousands of these entries.
Forward (Dual) 0.161 7.24e-8 Zygote is the slowest at 5.364 ms, despite also being
reverse mode. As discussed by the Julia community [14],
Reverse (Tape) 0.161 7.24e-8 Zygote generates the backward pass as a chain of closures,
each of which is a heap allocation. For a simulation with 400
Zygote 0.161 7.24e-8 timesteps and roughly 10 scalar operations per step, this pro­
duces thousands of heap objects per gradient call, confirmed
Enzyme 0.161 7.24e-8 by the high standard deviation of 6.685 ms indicating frequent
garbage collector interference. This behaviour is well-suited
All four methods return identical relative errors, confirming to machine learning workloads dominated by large matrix
that they agree with each other to machine precision. This operations, but is a poor fit for tight scalar simulation loops.
validates that the timing comparisons in the following sections
are between four correct implementations of the same compu­
tation.
c) Timing:
Next, we check how fast each method performs and how
good it is for the basic Lotka-Volterra experiment.
 Next, we have the tape simulation. This is broadly O(1)
with respect to the number of parameters as expected, though
the tape overhead grows slightly because a larger system
means more operations recorded per step. Notably at n=4 it is
actually slower than forward mode, for very small parameter
counts the tape allocation and backward pass overhead out­
weighs the cost of simply running a few extra forward passes.
The crossover point where reverse mode becomes faster than
forward mode occurs around n=8 parameters.
Forward mode grows roughly linearly. This is consistent
with its O(n) cost: each additional parameter requires one full
additional simulation pass.
Zygote scales linearly like forward mode but starts at a
Fig. 3. Bar diagram timing much higher baseline, because its per-operation heap alloca­
tion cost dominates regardless of parameter count.
d) Scaling: e) Code complexity:
To observe how each method behaves as the problem grows, Beyond performance, the four methods differ significantly
we extend the Lotka-Volterra system to n interacting species in how much code they require the programmer to write.
arranged in a ring, giving 2n parameters total, n growth rates Forward mode requires no changes to the simulation func­
and n interaction coefficients. We sweep from n=4 to n=24 tion at all. The plain Julia code runs unchanged because
parameters and measure the median gradient computation operator overloading on the Dual struct propagates derivatives
time over 50 samples at each size. automatically through every arithmetic operation. Computing
the gradient is a single call:
TABLE III
Median gradient computation time as the number of parameters dual_gradient(p) do args...
increases. simulate_lv_scalar(args...)
end
Parameters Forward Reverse Zygote Enzyme Zygote is equally simple from the user’s perspective, the
same unmodified simulation function is passed directly to
4 1.07 ms 2.14 ms 79.02 ms 0.10 ms [Link]. The entire gradient computation is one line:
[Link](simulate_lv_scalar, p...)
8 5.36 ms 3.80 ms 140.18 ms 0.20 ms
The simplicity of both forward mode and Zygote comes
from the same underlying reason: the simulation is written in
12 13.29 ms 4.32 ms 219.33 ms 0.23 ms
plain Julia and the AD system handles differentiation invisibly.
16 18.12 ms 7.15 ms 264.94 ms 0.27 ms Enzyme requires a small but explicit change. The simu­
lation must be rewritten to accept a Vector{Float64}
20 29.51 ms 10.02 ms 324.76 ms 0.23 ms instead of individual scalar arguments, and the caller must
manually allocate a shadow array and annotate each argument
24 36.01 ms 11.04 ms 405.70 ms 0.40 ms with Active, Const, or Duplicated:
d = zeros(Float64, length(p))
[Link](ReverseWithPrimal,
simulate_lv_vector, Active, Duplicated(p,
d))
Manual reverse mode is by far the most complex. The
simulation cannot be written in plain Julia, every single
arithmetic operation must be manually routed through the
tape, replacing natural expressions like 𝛼 ∗ 𝑥 − 𝛽 ∗ 𝑥 ∗ 𝑦 with
explicit tape calls:
s_αx = taped_mul(tape, s_α, s_x)
s_βxy = taped_mul(tape, taped_mul(tape,
s_β, s_x), s_y)
s_dx = taped_sub(tape, s_αx, s_βxy)
What is two lines of readable simulation code becomes ten
lines of tape manipulation. This verbosity scales with every
Fig. 4. Scaling graph plot operation in the simulation, making manual reverse mode
impractical for larger problems despite its correctness.
We can see that Enzyme stays stable, and this was to be
expected because everything is done in LLVM, so values stay
in CPU registers.
 IV. Conclusion
So overall Enzyme comes on top. It performs the best and
scales the best with many inputs. Besides that it is also fairly
easy to write.
If you want more inside on how your code works, you
could write it yourself. Depending on what the purpose is,
you might choose forward mode over reverse mode. Forward
mode is especially easy to implement and fairly fast for small
amount of inputs. But it scales fairly poorly.

References
[1] Wikipedia contributors, “Automatic Differentiation.” Wikipedia, The
Free Encyclopedia, 2025.
[2] J. Shi, “Automatic Differentiation: Forward and Reverse.” 2022.
[3] Wikipedia contributors, “Lotka–Volterra equations.” Wikipedia, The
Free Encyclopedia, 2025.
[4] J. Revels, M. Lubin, and T. Papamarkou, “Forward-Mode Automatic
Differentiation in Julia,” 2016, [Online]. Available: [Link]
abs/1607.07892
[5] JuliaDiff Contributors, “[Link]: Forward Mode Automatic Dif­
ferentiation for Julia.” [Online]. Available: [Link]
[Link]
[6] JuliaDiff Contributors, “[Link]: Reverse Mode Automatic Dif­
ferentiation for Julia.” [Online]. Available: [Link]
[Link]
[7] ApX Machine Learning, “Automatic Differentiation with [Link].”
2026.
[8] The [Link] Developers, “[Link] Documentation, Version 0.5.”
2020.
[9] The [Link] Developers, “[Link] Documentation (Stable).” 2026.
[10] DeepWiki, “[Link] – EnzymeAD/[Link] Overview.” 2025.
[11] Wikipedia contributors, “Euler Method.” Wikipedia, The Free Encyclo­
pedia, 2025.
[12] Petersd, “Euler's Method – Demo.” 2025.
[13] JuliaCI, “[Link] Documentation (Stable).” 2026.
[14] cortner, “Zygote Performance.” 2019.