줄리아에서 도함수 구하는 법
개요1
이름은 Calculus.jl인데 적분은 지원하지 않는다.
머신러닝등에서 말하는 자동 미분이 필요하다면 Ztgote.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
심볼릭 미분
심볼릭 미분은 SymEngine.jl 패키지에서도 사용할 수 있다.
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
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