ジュリアでの微分の求め方
概要1
名前はCalculus.jlだけど、積分はサポートしない。
機械学習などで話される自動微分が必要ならZygote.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
二階導関数 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