Dual Numbers for Automatic Differentiation

25 September 2026 · Download PDF · Download code

1. Scalar derivatives

Dual numbers are a key part of forward-mode differentiation. However, I find the introduction of them unintuitive. Usually, as done in [1], they are introduced such that

where ${𝜀}^2=0$. This is not incorrect. However, I think it is more natural to think of dual numbers, at least in the context of automatic differentiation, as a number type that simply tracks both $𝑓\left(𝑥\right)$ and ${𝑓}^{\prime }\left(𝑥\right)$ at the same time by introducing an $𝜀$. This is sort of similar to how a complex number can be seen to track an $𝑎$ and $𝑏$ at the same time by introducing an $𝑖$. However, for complex numbers the aim is to be able to solve ${𝑥}^2=-1$. For dual numbers the intent, as noted before, is to track $𝑓\left(𝑥\right)$ and ${𝑓}^{\prime }\left(𝑥\right)$ at the same time. We therefore define how a function acts on a dual number by the composition rule

of which (1) is the special case $𝑔\left(𝑥\right)+{𝑔}^{\prime }\left(𝑥\right)𝜀=𝑥+𝜀$. In the same way we define multiplication and division of dual numbers by the product and quotient rules

and addition and subtraction by the sum and difference rules

Now an initial complaint would be: What if ${𝑓}^{\prime }\left(𝑥\right)$ (or ${𝑔}^{\prime }\left(𝑥\right)$) is not known—which is usually the case when forward-mode differentiation is applied to a general function $𝑓$? The answer is that we don’t need to know them as long as the derivatives of the basic operations that make up $𝑓$ are known. That is, we can simply evaluate $𝑓\left(𝑥+𝜀\right)$, which by (2) should give us the derivative of $𝑓$ at $𝑥$ as the coefficient of $𝜀$. In the code example below we define the following composition rules for the elementary functions $sin$ and $exp$ as

Note that in general $𝑏$ does not have to come from the derivative of $𝑎$ with respect to some variable.

1
# A small, deliberately incomplete forward-mode AD implementation.
2
struct Dual{T<:Real}
3
    value::T
4
    deriv::T
5
end
6
# Promote mixed inputs such as Dual(2, 1.0) to a common type
7
Dual(value::Real, deriv::Real) = Dual(promote(value, deriv)...)
8
9
# Product rule: (f * g)' = f' * g + f * g'
10
Base.:*(f::Dual, g::Dual) = Dual(
11
    f.value * g.value,
12
    f.deriv * g.value + f.value * g.deriv,
13
)
14
Base.:*(f::Dual, g::Real) = f * Dual(g, zero(g))
15
Base.:*(f::Real, g::Dual) = Dual(f, zero(f)) * g
16
17
# Quotient rule: (f / g)' = (f' * g - f * g') / g^2
18
Base.:/(f::Dual, g::Dual) = Dual(
19
    f.value / g.value,
20
    (f.deriv * g.value - f.value * g.deriv) / g.value^2,
21
)
22
23
# Composition rule: (f ∘ g)' = (f' ∘ g) * g'
24
Base.sin(f::Dual) = Dual(sin(f.value), cos(f.value) * f.deriv)
25
Base.exp(f::Dual) = Dual(exp(f.value), exp(f.value) * f.deriv)
26
27
# Addition rules: (f + g)' = f' + g', (f + c)' = f', (c + f)' = f'
28
Base.:+(f::Dual, g::Dual) = Dual(f.value + g.value, f.deriv + g.deriv)
29
Base.:+(f::Dual, c::Real) = f + Dual(c, zero(c))
30
Base.:+(c::Real, g::Dual) = Dual(c, zero(c)) + g
31
32
# Subtraction rules: (-f)' = -f', (f - g)' = f' - g'
33
Base.:-(f::Dual) = Dual(-f.value, -f.deriv)
34
Base.:-(f::Dual, g::Dual) = f + (-g)
35
36
# The derivative is the epsilon coefficient of f(x + 1*epsilon)
37
derivative(f, x) = f(Dual(x, one(x))).deriv
38
39
# Dual is not a Number, so Base's generic ^ does not apply to it and
40
# f_scalar(x) = sin(x^2) + x would throw a MethodError
41
f_scalar(x) = sin(x * x) + x
42
@assert isapprox(derivative(f_scalar, 2.0), 4cos(4.0) + 1)
43
44
println("f'(2) = ", derivative(f_scalar, 2.0))

2. Jacobian-vector products

For a multivariate function $𝑓:{ℝ}^{𝑛}\to {ℝ}^{𝑚}$ the simplest extension keeps a single $𝜀$ and seeds each input ${𝑥}_{𝑗}$ with a direction ${𝑣}_{𝑗}$. The dual numbers then track the function value and a directional derivative (a Jacobian-vector product) at the same time, i.e.,

where ${𝐽}_{𝑓}\left(𝑥\right)$ is the Jacobian. This definition follows directly by looking at the $𝑖$th entry as

A thing to note here is that the sum is computed automatically by the dual number arithmetic, and hence we never form the Jacobian. From a scaling perspective, this is a big win as the Jacobian is $𝑚\times 𝑛$ in size, and hence can be very large.

The price is that one evaluation gives one product ${𝐽}_{𝑓}\left(𝑥\right)𝑣$. Its cost is a small multiple of the cost of evaluating $𝑓$, but the full Jacobian requires $𝑛$ evaluations, one for each unit vector $𝑣={𝑒}_{𝑗}$. Forward mode is therefore well suited when $𝑛$ is small compared to $𝑚$. In the opposite case reverse mode (backpropagation) is preferred. It computes one vector-Jacobian product ${𝑤}^{\top }{𝐽}_{𝑓}\left(𝑥\right)$ per evaluation, also at a small multiple of the cost of evaluating $𝑓$, so the full Jacobian requires $𝑚$ evaluations, one for each unit vector $𝑤={𝑒}_{𝑖}$. The price is that the intermediate values of $𝑓$ must be stored for the backward pass. The common case of a scalar loss of many parameters ($𝑚=1$) needs only a single evaluation; see Backpropagation as a Linear Solve.¹¹ In the notation of that post, the Jacobian of a network is $𝑃{𝐿}^{-1}{𝑀}^{\top }$, where $𝐿$ is block lower triangular. Forward mode multiplies $𝑣$ from the right, $𝑃\left({𝐿}^{-1}\left({𝑀}^{\top }𝑣\right)\right)$, so the solve with $𝐿$ is a forward substitution that follows the layers in order, exactly as the dual numbers do. Reverse mode multiplies ${𝑤}^{\top }$ from the left, $\left(\left({𝑤}^{\top }𝑃\right){𝐿}^{-1}\right){𝑀}^{\top }$, so the solve becomes the backward substitution ${𝐿}^{\top }𝜆={𝑃}^{\top }𝑤$.

1
@isdefined(Dual) || include("dual_numbers.jl")
2
3
function jvp(f, x, v)
4
    @assert length(x) == length(v)
5
    seeded = [Dual(xi, vi) for (xi, vi) in zip(x, v)]
6
    result = f(seeded)
7
    return map(z -> z.value, result), map(z -> z.deriv, result)
8
end
9
10
f_vector(x) = (x[1] * x[2] + sin(x[1]), exp(x[2]))
11
x, v = [2.0, 3.0], [1.0, -1.0]
12
value, deriv = jvp(f_vector, x, v)
13
@assert all(isapprox.(value, (6 + sin(2.0), exp(3.0))))
14
@assert all(isapprox.(deriv, (1 + cos(2.0), -exp(3.0))))
15
16
println("f($x) = ", value)
17
println("Jf($x)$v = ", deriv)

@misc{paltorp2026dualnumbers,
  author = {Paltorp, Mikkel},
  title = {{Dual Numbers for Automatic Differentiation}},
  year = {2026},
  month = sep,
  url = {https://mipals.github.io/pdf/2026/09-25-dual-numbers/dual-numbers.pdf},
  note = {Blog post}
}