Implementing Automatic Differentiation Forward Mode Using Dual Numbers in Julia

Overview

The forward mode of automatic differentiation can be easily implemented using dual numbers. This post explains how to implement forward mode in Julia. For background knowledge on dual numbers and automatic differentiation, the following posts are recommended.

• dual numbers
  • Differentiable real functions defined over the ring of dual numbers
• automatic differentiation
  • Automatic differentiation and dual numbers

Code ¹

This is an example of computing automatic differentiation for the function $y(x) = \ln (x^{2} + \sin x)$.

Defining the dual number struct

First, let us tell the computer what a dual number is. In Julia, a dual number can be defined as a struct as follows. With applications in automatic differentiation in mind, let the first component v be the value (function value), and the second component d be the derivative.

struct Dual 
    v::Float64 # (function) value
    d::Float64 # derivative
end

Now we can define the dual numbers $x = (3, 1)$ and $y = (2, 0)$.

julia> x = Dual(3, 1)
Dual(3.0, 1.0)

julia> x.v, x.d
(3.0, 1.0)

julia> y = Dual(2, 0)
Dual(2.0, 0.0)

Defining addition of dual numbers

Now let us add x and y. Of course, this raises an error. This is because the binary operation + built into Julia by default is not defined for Dual and Dual.

julia> x + y
ERROR: MethodError: no method matching +(::Dual, ::Dual)

Currently, there are 189 methods defined for the addition +.

julia> methods(+)
# 189 methods for generic function "+" from Base:

If we define + for Dual and Dual as below, we can see that one more method has been added.

julia> Base.:+(x::Dual, y::Dual) = Dual(x.v + y.v, x.d + y.d)

julia> methods(+)
# 190 methods for generic function "+" from Base:

Now we can compute x + y.

julia> x + y
Dual(5.0, 1.0)

Defining multiplication of dual numbers

In the same way as addition, we can define and compute multiplication as follows.

julia> Base.:*(x::Dual, y::Dual) = Dual(x.v*y.v, x.v*y.d + x.d*y.v)

julia> x * y
Dual(6.0, 2.0)

Defining functions over dual numbers

Since we need to compute the function $y(x) = \ln (x^{2} + \sin x)$, let us define the logarithm function and the sine function for dual numbers.

Base.sin(x::Dual) = Dual(sin(x.v), x.d * cos(x.v) )
Base.log(x::Dual) = Dual(log(x.v), x.d / x.v)

Forward mode computation

Now, by computing as follows, we can obtain the function value and the derivative simultaneously.

julia> y₁ = x*x
Dual(9.0, 6.0)

julia> y₂ = sin(x)
Dual(0.1411200080598672, -0.9899924966004454)

julia> y₃ = y₁ + y₂
Dual(9.141120008059866, 5.010007503399555)

julia> y₄ = log(y₃)
Dual(2.2127829171337874, 0.548073704204972)

julia> log(3^2 + sin(3))
2.2127829171337874

julia> (2*3 + cos(3)) / (3^2 + sin(3))
0.548073704204972

Complete Code

# 이원수 구조체 정의
struct Dual 
    v::Float64 # (function) value
    d::Float64 # derivative
end

x = Dual(3, 1)
y = Dual(2, 0)

# 이원수의 덧셈 정의
x + y

methods(+)

Base.:+(x::Dual, y::Dual) = Dual(x.v + y.v, x.d + y.d)
methods(+)

x + y

# 이원수의 곱셈 정의
Base.:*(x::Dual, y::Dual) = Dual(x.v*y.v, x.v*y.d + x.d*y.v)
x * y

# 이원수 위의 함수 정의
Base.sin(x::Dual) = Dual(sin(x.v), x.d * cos(x.v) )
Base.log(x::Dual) = Dual(log(x.v), x.d / x.v)

# 전진모드 계산
y₁ = x*x
y₂ = sin(x)
y₃ = y₁ + y₂
y₄ = log(y₃)

log(3^2 + sin(3))
(2*3 + cos(3)) / (3^2 + sin(3))

Environment

• OS: Windows11
• Version: Julia 1.7.1

See Also

• dual numbers
  • Differentiable real functions defined over the ring of dual numbers
• automatic differentiation
  • Automatic differentiation and dual numbers
  • Implementing automatic differentiation forward mode with dual numbers in Julia