📌

Optim.jlを使用した最適化

2022/01/15に公開

はじめに

この記事では、Optim.jlというJulia言語の制約なし非線形最適化パッケージの使用方法を説明する。特に、データを使用した (例えば最尤法の実行方法など) 最適化の実行方法について説明する。また、読者としてはJulia言語の仕様 (文法) を一通り知っている人を想定する。

Optim.jlとは?

Optim.jl とは、Julia言語の制約なし非線形最適化パッケージである。このパッケージでは、非線形最適化の手法である Nelder–Mead法, 共役勾配法, 準ニュートン法など関数が提供されている。また、The Journal of Open Source Softwareによって査読済みであり、関数の信頼性は高いと言える。現代科学では、最適化理論は機械学習やスケジューリング問題をはじめとする様々な分野の基盤を担っており、その重要度は日々に増加している。そこで、この記事では、Optim.jlを使用した最適化手法について説明しよう。具体的には、

  1. 関数に対する最適化
    • Nelder–Mead法
    • 共役勾配法
    • 準ニュートン法
  2. データを使用した最適化

の2項目を説明する。

関数に対する最適化

まず、通常の関数に対する最適化問題について説明しよう。ここでは、

y = (x-2)^2 + 4

という関数に対する最小化問題を例として挙げよう。この関数の最小値の理論解は、x = 2のときy=4という値をとる。この理論解を踏まえて、Optimで提供されている最適化手法を使用した数値解の導出を行おう。

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の関数 j のフォーマットが決まっている。具体的には、j(jaccobian::Vector,最適化パラメター, other independent variables...) と決まっている。そのため、下記のようにして共役勾配法を実行する。

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の関数 j のフォーマットが決まっている。具体的には、h(Hessian::Matrix,最適化パラメター, other independent variables...) と決まっている。そのため、下記のようにして共役勾配法を実行する。

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の最適化関数に対して引数を渡すことが出来る。

最後に、データを渡す最適化手法の具体例を提示しよう。ここでは、

y = (x-2)^2 + 4 + \epsilon,

where \epsilon\sim N(0,1). という式で発生させたデータセットを

\argmin_{\beta} \sum (y - x\beta)^2

を実現する\betaを推定する方法を説明しよう。

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