How to Find Derivatives in Julia
Overview1
The package is named Calculus.jl, but it does not support integration.
If automatic differentiation, as discussed in machine learning, is needed, refer to the Zygote.jl package.
Differentiation of Single Variable Function
Derivative function derivative()
It calculates the derivative of $f : \R \to \R$.
derivative(f)orderivative(f, :x): Returns the derivative $f^{\prime}$.derivative(f, a): Returns the differential coefficient $f^{\prime}(a)$.
julia> f(x) = 1 + 2x + 3x^2
f (generic function with 1 method)
julia> g(x) = sin(x)
g (generic function with 1 method)
julia> derivative(f)
#1 (generic function with 1 method)
julia> derivative(f, 1)
7.99999999996842
julia> Df = derivative(f)
#1 (generic function with 1 method)
julia> Dg = derivative(g)
#1 (generic function with 1 method)
#f'(x) = 2 + 6x
julia> Df(1)
7.99999999996842
#g'(x) = cos x
julia> Dg(pi)
-0.9999999999441258
Composite functions can also be differentiated.
#f∘g(x) = (2 + 6 sin x)cos x
julia> derivative(f∘g)
#1 (generic function with 1 method)
julia> derivative(f∘g, pi/4)
4.414213562300037
julia> (2+6sin(pi/4))cos(pi/4)
4.414213562373095
Second derivative function second_derivative()
It calculates the second derivative of $f : \R \to \R$.
The function returned by derivative() can take integers as input, but second_derivative() cannot use integer types. Irrational number types are also not allowed.
julia> derivative(f, 1)
7.99999999996842
julia> second_derivative(f, 1)
ERROR: MethodError: no method matching eps(::Type{Int64})
Closest candidates are:
eps() at float.jl:764
eps(::AbstractFloat) at float.jl:760
eps(::Type{Float16}) at float.jl:761
julia> second_derivative(f, 1.)
5.9999956492003905
julia> second_derivative(g, pi)
ERROR: MethodError: no method matching eps(::Type{Irrational{:π}})
Closest candidates are:
eps() at float.jl:764
eps(::AbstractFloat) at float.jl:760
eps(::Type{Float16}) at float.jl:761
julia> second_derivative(g, convert(Float64, pi))
-1.3553766145945872e-7
julia> second_derivative(g, 1pi)
-1.3553766145945872e-7
Differentiation of Multivariable Function
Gradient gradient()
Returns the gradient of $f : \mathbb{R}^{n} \to \mathbb{R}$.
Note that when defining a multivariable function, it should not be defined as a function with multiple real variables. It must be defined as a single-variable function that takes a vector as input. If it’s not a function accepting a vector as input, differentiation can be performed, but the value cannot be calculated. For example, it should not be defined as $f_{1}$, but rather as $f_{2}$.
julia> f₁(x,y,z) = x*y + z^2
f₁ (generic function with 1 method)
julia> Calculus.gradient(f₁)
#2 (generic function with 1 method)
julia> ∇f₁ = Calculus.gradient(f₁)
#2 (generic function with 1 method)
julia> ∇f₁(1,1,1)
ERROR: MethodError: no method matching (::Calculus.var"#2#4"{typeof(f₁), Symbol})(::Int64, ::Int64, ::Int64)
julia> ∇f₁([1,1,1])
ERROR: MethodError: no method matching f₁(::Vector{Float64})
julia> Calculus.gradient(f₁, 1,1,1)
ERROR: MethodError: no method matching gradient(::typeof(f₁), ::Int64, ::Int64, ::Int64)
julia> Calculus.gradient(f₁, [1,1,1])
ERROR: MethodError: no method matching f₁(::Vector{Float64})
julia> f₂(x) = x[1]*x[2] + x[3]^2
f₂ (generic function with 1 method)
julia> Calculus.gradient(f₂, [1,1,1])
3-element Vector{Float64}:
1.0000000000235538
1.0000000000235538
1.9999999999737708
julia> ∇f₂ = Calculus.gradient(f₂)
#2 (generic function with 1 method)
julia> ∇f₂(1,1,1)
ERROR: MethodError: no method matching (::Calculus.var"#2#4"{typeof(f₂), Symbol})(::Int64, ::Int64, ::Int64)
julia> ∇f₂([1,1,1])
3-element Vector{Float64}:
1.0000000000235538
1.0000000000235538
1.9999999999737708
Hessian hessian()
Returns the Hessian of $f : \mathbb{R}^{n} \to \mathbb{R}$.
Similar to second_derivative(), it only accepts Float data types as input. Like gradient(), it can only return values for functions that take vectors as input.
julia> hessian(f₂, [1.,1.,1.])
3×3 Matrix{Float64}:
3.88008e-7 1.0 0.0
1.0 3.88008e-7 0.0
0.0 0.0 2.0
julia> H = hessian(f₂)
#7 (generic function with 1 method)
julia> H([1.,1.,1.])
3×3 Matrix{Float64}:
3.88008e-7 1.0 0.0
1.0 3.88008e-7 0.0
0.0 0.0 2.0
Jacobian jacobian()
Returns the Jacobian of $f : \mathbb{R}^{m} \to \mathbb{R}^{n}$.
Similar to second_derivative(), it only accepts Float data types as input. Like gradient(), it can only return values for functions that take vectors as input.
Unlike other functions, using it like jacobian(f, [x, y, z]) is not possible.
julia> h(x) = [x[1], x[1]*x[2], x[1]*x[3]^2]
h (generic function with 2 methods)
julia> jacobian(h, [1.,1.,1.])
ERROR: MethodError: no method matching jacobian(::typeof(h), ::Vector{Float64})
julia> Jh = jacobian(h)
(::Calculus.var"#g#5"{typeof(h), Symbol}) (generic function with 1 method)
julia> Jh([1.,1.,1.])
3×3 Matrix{Float64}:
1.0 0.0 0.0
1.0 1.0 0.0
1.0 0.0 2.0
Symbolic Differentiation
Symbolic differentiation is also available in the SymEngine.jl package.
differentiate()
Performs symbolic differentiation.
While it returns constant terms and $x$ neatly, for cases like $ax$ or $x^{n}$, it returns in the form of the product rule. For example, differentiating $3x^{2}$ would return it as the product of $3$ and $x^{2}$, seen as $\dfrac{d 3}{dx} x^{2} + 3\dfrac{d x^{2}}{dx}$. Even $x^{2}$ is seen as the product of $1$ and $x^{2}$.
julia> differentiate("1", :x)
0
julia> differentiate("1 + x", :x)
1
julia> differentiate("x^2", :x)
:(2 * 1 * x ^ (2 - 1))
julia> differentiate("x^3", :x)
:(3 * 1 * x ^ (3 - 1))
julia> differentiate("1 + x + x^2", :x)
:(1 + 2 * 1 * x ^ (2 - 1))
julia> differentiate("1 + 2x + 3x^2 + 4x^3", :x)
:((0 * x + 2 * 1) + (0 * x ^ 2 + 3 * (2 * 1 * x ^ (2 - 1))) + (0 * x ^ 3 + 4 * (3 * 1 * x ^ (3 - 1))))
julia> differentiate("x^2 * sin(x) + exp(x) * cos(x)", :x)
:(((2 * 1 * x ^ (2 - 1)) * sin(x) + x ^ 2 * (1 * cos(x))) + ((1 * exp(x)) * cos(x) + exp(x) * (1 * -(sin(x)))))
Characters not specified as symbols are treated as constants, and if more than one symbol is input, it returns the differentiation for each symbol. However, writing "3yx" treats xy as a single variable, so it’s important to explicitly denote multiplication by writing it as 3x*y.
julia> differentiate("1 + x + 3yx + y^2", :x)
:(1 + (0yx + 3 * 0))
julia> differentiate("1 + x + 3y*x + y^2", :x)
:(1 + ((0y + 3 * 0) * x + (3y) * 1))
julia> differentiate("1 + x + 3yx + y^2", [:x, :y])
2-element Vector{Any}:
:(1 + (0yx + 3 * 0))
:((0yx + 3 * 0) + 2 * 1 * y ^ (2 - 1))
julia> differentiate("1 + x + 3x*y + y^2", [:x, :y])
2-element Vector{Any}:
:(1 + ((0 * x + 3 * 1) * y + (3x) * 0))
:(((0 * x + 3 * 0) * y + (3x) * 1) + 2 * 1 * y ^ (2 - 1))
simplify()
The return of differentiate() lacks readability, but simplify() neatly organizes it. However, it doesn’t perform properly if vectors are used as input.
julia> simplify(differentiate("1 + 2x + 3x^2 + 4x^3", :x))
:(2 + 3 * (2x) + 4 * (3 * x ^ 2))
julia> simplify(differentiate("1 + x + 3x*y + y^2", :x))
:(1 + 3y)
julia> simplify(differentiate("1 + x + 3x*y + y^2", :y))
:(3x + 2y)
julia> simplify(differentiate("1 + x + 3x*y + y^2", [:x, :y]))
2-element Vector{Any}:
:(1 + ((0 * x + 3 * 1) * y + (3x) * 0))
:(((0 * x + 3 * 0) * y + (3x) * 1) + 2 * 1 * y ^ (2 - 1))
deparse()
Turns the return of symbolic differentiation into a string.
julia> a = differentiate("1 + x + 3x*y + y^2", :x)
:(1 + ((0 * x + 3 * 1) * y + (3x) * 0))
julia> deparse(a)
"1 + ((0 * x + 3 * 1) * y + (3 * x) * 0)"
julia> deparse(simplify(a))
"1 + 3 * y"
Verification
check_derivative()
One can check how much the derivative obtained by derivative() differs from the actual derivative. It’s implemented for the four types of derivatives excluding jacobian().
check_derivative(f, Df, a): Returns the absolute value ofderivative(f, a)-Df(a).
julia> f(x) = 1 + x^2
f (generic function with 1 method)
julia> Df(x) = 2x
Df (generic function with 1 method)
julia> Calculus.check_derivative(f, Df, 1)
2.6229241001374248e-11
Application
Polynomials.jl
Although derivative is implemented in Polynomials.jl itself, the derivatives can also be calculated with Calculus.derivative().
julia> p = Polynomial([1,2,4,1])
Polynomial(1 + 2*x + 4*x^2 + x^3)
julia> Polynomials.derivative(p)
Polynomial(2 + 8*x + 3*x^2)
julia> Calculus.derivative(p)
#1 (generic function with 1 method)
julia> Polynomials.derivative(p,2)
Polynomial(8 + 6*x)
julia> Calculus.second_derivative(p)
#6 (generic function with 1 method)
Environment
- OS: Windows10
- Version: Julia 1.6.2, Calculus 0.5.1