줄리아에서 도함수 구하는 법

줄리아에서 도함수 구하는 법

How to Use Calculus.jl in Julia

개요1

이름은 Calculus.jl인데 적분은 지원하지 않는다.

일변수 함수의 미분

도함수 derivative()

$f : \R \to \R$의 도함수를 구해준다.

  • derivative(f) 혹은 derivative(f, :x): 도함수 $f^{\prime}$을 리턴한다.
  • derivative(f, a): 미분계수 $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

합성함수도 미분할 수 있다.

#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

2계도함수 second_derivative()

$f : \R \to \R$의 이계도함수를 구해준다.

derivative()로 반환된 함수는 정수를 입력값으로 쓸 수 있지만, second_derivative()는 정수형 타입을 쓸 수 없다. 무리수 타입도 쓸 수 없다.

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

다변수 함수의 미분

그래디언트 gradient()

$f : \mathbb{R}^{n} \to \mathbb{R}$의 그래디언트를 리턴한다.

주의해야 할 점은 다변수 함수를 정의할 때 진짜 변수가 여러개인 함수로 정의하면 안된다는 것이다. 벡터를 입력으로 받는 일변수 함수로 정의해야 한다. 벡터를 입력으로 받는 함수가 아니면 미분은 할 수 있지만, 값을 계산할 수 없다는 것이다. 다음의 예에서 $f_{1}$과 같이 정의하면 안되고, $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()

$f : \mathbb{R}^{n} \to \mathbb{R}$의 헤시안를 리턴한다.

second_derivative()와 마찬가지로 Float 자료형만 입력을 받을 수 있다. gradient()와 마찬가지로 벡터를 입력으로 받는 함수에 대해서만 값을 리턴받을 수 있다.

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()

$f : \mathbb{R}^{m} \to \mathbb{R}^{n}$의 자코비안를 리턴한다.

second_derivative()와 마찬가지로 Float 자료형만 입력을 받을 수 있다. gradient()와 마찬가지로 벡터를 입력으로 받는 함수에 대해서만 값을 리턴받을 수 있다.

다른 함수들과는 다르게 jacobian(f, [x, y, z])와 같이 사용하는 건 불가능하다.

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

심볼릭 미분

differentiate()

심볼릭한 미분을 수행한다.

상수항과 $x$는 깔끔하게 리턴하는 반면에, $ax$나 $x^{n}$과 같은 경우에는 곱의 미분법의 꼴로 리턴한다. 가령 $3x^{2}$을 미분해보면 이를 $3$과 $x^{2}$의 곱으로 보고 $\dfrac{d 3}{dx} x^{2} + 3\dfrac{d x^{2}}{dx}$와 같이 리턴한다는 말이다. 심지어 $x^{2}$도 $1$과 $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)))))

입력한 심볼이 아닌 문자는 상수로 처리하고, 두 개 이상의 심볼을 입력한 경우에는 각 심볼에 대한 미분을 리턴한다. 이때, "3yx"라고 적으면 xy 자체를 하나의 변수로 보기 때문에 꼭 곱하기 기호를 넣어 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()

differentiate()의 리턴은 가독성이 떨어지는데, simplify()는 이를 깔끔하게 정리해준다. 단, 벡터를 입력으로 쓰면 제대로 수행하지 못한다.

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()

심볼릭 미분의 리턴을 문자열로 바꿔준다.

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"

검산

check_derivative()

derivative()로 구한 도함수가 진짜 도함수와 얼마나 차이가 나는지를 확인할 수 있다. jacobian()을 제외한 다른 4개의 도함수에 대해서 구현되어있다.

  • check_derivative(f, Df, a): derivative(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

응용

Polynomials.jl

🔒(06/08)Polynomials.jl 자체에도 derivative가 구현되어있기는 하지만 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)

환경

  • OS: Windows10
  • Version: Julia 1.6.2, Calculus 0.5.1

  1. https://github.com/JuliaMath/Calculus.jl ↩︎

댓글