Largely inspired from Prof. Alan Edelman's talk Automatic Differentiation in 10 minutes with Julia
In linear algebra, the dual numbers extend the real numbers by adjoining one new element $\varepsilon$ (epsilon) with the property $\varepsilon^2 = 0$. Thus the multiplication of dual numbers is given by: $$(a+b\varepsilon )(c+d\varepsilon )=ac+(ad+bc)\varepsilon,$$ and addition is done componentwise: $$(a+b\varepsilon )+(c+d\varepsilon )=(a+c)+(b+d)\varepsilon.$$
Using Julia, we can create a data type for dual numbers as follows:
struct Dual <: Number
value::Float64
epsilon::Float64
end
We now define addition, substraction and multiplication for dual numbers. To do this, we define a new method to the Julia Base operators +,-,*
.
import Base: +, -, *
+(z::Dual, w::Dual) = Dual(z.value+w.value, z.epsilon+w.epsilon)
-(z::Dual, w::Dual) = Dual(z.value-w.value, z.epsilon-w.epsilon)
*(z::Dual, w::Dual) = Dual(z.value*w.value, z.value*w.epsilon+z.epsilon*w.value)
For convenience, we also define a new display function for dual numbers.
Base.show(io::IO,x::Dual) = print(io,x.value," + ",x.epsilon," ε")
a=Dual(2,1)
We are now ready to start playing with dual numbers.
a*a
a^2
a^7
This last two commands are rather surprising as we never define the power of a dual numbe! But it turns out that in Julia power is defined for any x
supporting *
as can be seen here or by following the link given by the following command:
@which a^7
2*a
To correct this problem, we need to:
2
into the Dual(2,0)
This is what is called conversion and promotion and here it can be done as follows:
import Base: convert, promote_rule
convert(::Type{Dual}, x::Real) = Dual(x,zero(x))
promote_rule(::Type{Dual}, ::Type{<:Number}) = Dual
2*a
a-1
b=Dual(1,1)
3*(a+b)^2
We can now get derivatives for polynomials. Consider $P(x) = p_0+p_1x+p_2x^2+\dots +p_n x^n$, note that since $\varepsilon^2 =0$ ($\varepsilon$ is nilpotent), we have $(x+\varepsilon y)^k = x^k + kx^{k-1}\varepsilon y$. Hence, we have $$ P(x+\varepsilon y) = P(x) + \left(p_1 + 2p_2 x+ \dots np_n x^{n-1}\right)y\varepsilon = P(x) + P'(x)y\varepsilon $$
More info on Wikipedia
a
3*(1+a)^2
value(z::Dual) = z.value
epsilon(z::Dual) = z.epsilon
epsilon(3*(1+a)^2)
We can now define a simple function to compute the derivative of a polynomials at a given point as follows:
function derivative(f,x::Real)
epsilon(f(Dual(x,1)))
end
derivative(x->1+x+3x^2,1)
a^(1/2)
To deal with this problem, we follow the Babylonian as described in this nice tutorial
Repeat $ t \leftarrow \frac{1}{2}\left(t+\frac{x}{t}\right)$ until $t$ converges to $\sqrt{x}$.
function Babylonian(x; N = 10)
t = (1+x)/2
for i = 2:N; t=(t+x/t)/2 end
t
end
Babylonian(2), √2 # Type \sqrt+<tab> to get the symbol
The Babylonian algorithm uses only addition and division, hence we need to define a division for dual numbers: $$\frac{a+b\varepsilon}{c+d\varepsilon}= \frac{a}{c}\left(1+\frac{b}{a}\varepsilon\right)\left(1-\frac{d}{c}\varepsilon\right)= \frac{a}{c} +\frac{bc-ad}{c^2}\varepsilon.$$
import Base:/
/(x::Dual, y::Dual) = Dual(x.value/y.value, (y.value*x.epsilon - x.value*y.epsilon)/y.value^2)
a
1/(1+a)
Babylonian(a), √2, 0.5/√2
derivative(x->Babylonian(x),2)
function dBabylonian(x; N = 10)
t = (1+x)/2
dt = 1/2
for i = 1:N;
t = (t+x/t)/2;
dt = (dt+(t-x*dt)/t^2)/2;
end
dt
end
x = 2; dBabylonian(x), .5/√x
derivative(x->Babylonian(x)^4,4)
derivative(x->(1+3*Babylonian(x))^3/Babylonian(x),2)
using Pkg; Pkg.add("ForwardDiff")
using ForwardDiff
ForwardDiff.derivative(sqrt, 2)
ForwardDiff.derivative(Babylonian, 2)
ForwardDiff.derivative(x->(1+3*sqrt(x))^3/sqrt(x),2)
For more on dual numbers, have a look at the (not active anymore) Julia package DualNumbers.jl