Optim.jlを使用した最適化
はじめに
この記事では、Optim.jlというJulia言語の制約なし非線形最適化パッケージの使用方法を説明する。特に、データを使用した (例えば最尤法の実行方法など) 最適化の実行方法について説明する。また、読者としてはJulia言語の仕様 (文法) を一通り知っている人を想定する。
Optim.jlとは？
Optim.jl とは、Julia言語の制約なし非線形最適化パッケージである。このパッケージでは、非線形最適化の手法である Nelder–Mead法, 共役勾配法, 準ニュートン法など関数が提供されている。また、The Journal of Open Source Softwareによって査読済みであり、関数の信頼性は高いと言える。現代科学では、最適化理論は機械学習やスケジューリング問題をはじめとする様々な分野の基盤を担っており、その重要度は日々に増加している。そこで、この記事では、Optim.jlを使用した最適化手法について説明しよう。具体的には、
- 関数に対する最適化
- Nelder–Mead法
- 共役勾配法
- 準ニュートン法
- データを使用した最適化
の2項目を説明する。
関数に対する最適化
まず、通常の関数に対する最適化問題について説明しよう。ここでは、
という関数に対する最小化問題を例として挙げよう。この関数の最小値の理論解は、
Nelder-Mead法
Nelder-Mead法とは、無制約非線形最適化問題を解くためのアルゴリズムである。このアルゴリズムでは、目的関数さえ定義すれば実行可能な最適化アルゴリズムである。そのため、後に紹介する方法と比較すると精度・収束速度ともに遅い。
それでは、Nelder-Mead法の実行方法について説明しよう。Optimでは、"optimize"という関数を使用して最適化を実行する。具体的には、
using Optim
# target variable
f(x) = (x[1]-2)^2 + 4
#template: optmize(目的関数, 初期値, method = NelderMead())
optimize(f, [1.], method = NelderMead())
共役勾配法
共役勾配法とは、Nelder-Mead法と異なり勾配 (Jacobbian) を必要とする。その代わり、Nelder-Mead方と比較して精度・収束速度ともにより良いパフォーマンスを持つ。
それでは、共役勾配法の実行方法について説明しよう。Optimで使用可能なJacobbianの関数
using Optim
# target variable
f(x) = (x[1] - 2)^2 + 4
function j(J::Vector, x)
J[1] = 2 *(x[1] - 2)
end
#template: optmize(目的関数, jacobbian, 初期値, method = ConjugateGradient())
optimize(f,j , [1.], method = ConjugateGradient())
準ニュートン法
準ニュートン法では、Hessian, Jacobbianの2つを入力する必要がある。最もよく使われている最適化手法であり、精度・収束速度ともに優れている。
それでは、準ニュートン法の実行方法について説明しよう。Optimで使用可能なHessianの関数
using Optim
# target variable
f(x) = (x[1] - 2)^2 + 4
function j(J, x)
J[1] = 2 *(x[1] - 2)
end
function h(H::Matrix, x)
H[1,1] = 2
end
#template: optmize(目的関数, jacobbian, 初期値, method = ConjugateGradient())
optimize(f,j , H, [1.], method = Newton())
データを使用した最適化
Optim.jlでは、python言語のscipy.optimizeで提供されているようなkwarg={"x":x}のようなフォーマットを使用したデータを渡すための引数が用意されていない。そのため、Optimでデータを使用した最適化を行うためには、function-like objectを使用する必要がある。function-like object とは、「Julia上で定義した任意のオブジェクトは、関数の機能を同時に保有することが出来るという性質」のことである。例えば、構造体 struct に対してfunction-like objectの性質を同時にもたせる方法を説明しよう。
# I strongly recommend you to define the type of struct obviously.
struct test{Type}
x::Type
end
s = test(2.0)
# add function onto test object.
function (s::test)(x)
# test型の構造体 sが参照可能
# 今回の例だと、s.x = 2.0を参照できる
return s.x^(x)
end
このようにして、function-like objectを定義することが出来る。その最大の特性は、objectの状態を関数内で引数として取ることが出来るという特性だ。たとえば、今回定義した「test」という構造体 "s" の挙動について例をあげよう。ここで定義された function-like objectである"s"というオブジェクトは、関数内で宣言時の状態 "s.x = 2"という状態が参照可能である。そのため、定義されたs(x)という関数は"s.x^x = 2^x"という値を返す。この性質を利用すると、Optim.jlの最適化関数に対して引数を渡すことが出来る。
最後に、データを渡す最適化手法の具体例を提示しよう。ここでは、
where
を実現する
struct target
y::Vector
x::Matrix
end
x = collect(-10:0.1:100)
y = (x .-2).^2 .+ 4 .+ randn(length(x))
X = ones(length(x),3)
X[:,2] .= x
X[:,3] .= x.^2
F = target(y,X)
function (F::target)(β)
return sum( (F.y .- F.x *β).^2 )
end
JJ = target(y,X)
function (JJ::target)(G,β)
G[1] = sum( -2 .*(JJ.y - JJ.x * β) .* JJ.x[:,1] )
G[2] = sum( -2 .*(JJ.y - JJ.x * β) .* JJ.x[:,2] )
G[3] = sum( -2 .*(JJ.y - JJ.x * β) .* JJ.x[:,3] )
end
HH = target(y,X)
function (HH::target)(H::Matrix,β)
H[1,1] = sum( -2 .*(HH.x[:,1] .* HH.x[:,1]) )
H[1,2] = sum( 2 .* (HH.x[:,2] .* HH.x[:,1]) )
H[1,3] = sum( 2 .* (HH.x[:,3] .* HH.x[:,1]) )
H[2,1] = sum( 2 .* (HH.x[:,1] .* HH.x[:,2]) )
H[2,2] = sum( -2 .*(HH.x[:,2] .* HH.x[:,2]) )
H[2,3] = sum( 2 .*(HH.x[:,3] .* HH.x[:,2]) )
H[3,1] = sum( 2 .* (HH.x[:,1] .* HH.x[:,3]) )
H[3,2] = sum( 2 .* (HH.x[:,2] .* HH.x[:,3]) )
H[3,3] = sum( -2 .*(HH.x[:,3] .* HH.x[:,3]) )
end
opt = optimize(F,JJ,HH,[1.,1.,1.], method= Newton())
Optim.minimizer(opt)
この手続きを踏めば、データを使用した最適化手法が使用可能である。
Discussion