🎒
Optim.jl入門
Optim.jlの使い方を簡単に解説します.
インストール
REPLまたはノートブック上でusing Pkg; Pkg.add("Optim")を実行するか, パッケージモードでadd Optimを実行してもインストールできます. ノートブック上ではusing Optimを宣言します. 可視化に使用するため, 同様にPlots.jlをインストール・宣言しておきます. 最後にusing Printfも使いますが, こちらは標準ライブラリなのでインストールは不要です.
# using Pkg
# Pkg.add("Optim")
using Optim
# using Pkg
# Pkg.add("Plots")
using Plots
using Printf
Hello World!
Rosenbrock関数を最小化する例です. デフォルトではNelder-Mead法が採用されていることがわかります.
# using Optim
f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
x0 = [0.0, 0.0]
optimize(f, x0)
* Status: success
* Candidate solution
Minimizer: [1.00e+00, 1.00e+00]
Minimum: 3.525527e-09
* Found with
Algorithm: Nelder-Mead
Initial Point: [0.00e+00, 0.00e+00]
* Convergence measures
√(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 60
f(x) calls: 118
特定の値を取り出す
結果の中から特定の値を取り出すには, res = optimize(f, x0)のように戻り値を受け取り, Optim.summary(), Optim.minimizer(), Optim.minimum()などに戻り値resを渡すと取り出せます. 詳細はこちら.
# using Optim
f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
x0 = [0.0, 0.0]
res = optimize(f, x0)
println("手法 ：", Optim.summary(res))
println("最小値 ：", Optim.minimum(res))
println("最小点 ：", Optim.minimizer(res))
println("反復回数：", Optim.iterations(res))
println("評価総数：", Optim.f_calls(res))
println("収束判定：", Optim.converged(res))
手法 ：Nelder-Mead
最小値 ：3.5255270584829996e-9
最小点 ：[0.9999634355313174, 0.9999315506115275]
反復回数：60
評価総数：118
収束判定：true
手法を指定する
手法はoptimize(f, x0, method=NelderMead())のように指定できます. デフォルトはNelderMead()なので, ここではLBFGS()を試します. 反復回数, 評価回数どちらもLBFGS()の方が少なくなっているように見え, 結果も正確です.
# using Optim
f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
x0 = [0.0, 0.0]
res = optimize(f, x0, method=LBFGS())
println("手法 ：", Optim.summary(res))
println("最小値 ：", Optim.minimum(res))
println("最小点 ：", Optim.minimizer(res))
println("反復回数：", Optim.iterations(res))
println("評価総数：", Optim.f_calls(res))
println("収束判定：", Optim.converged(res))
手法 ：L-BFGS
最小値 ：5.378388330692143e-17
最小点 ：[0.9999999926662504, 0.9999999853325008]
反復回数：24
評価総数：67
収束判定：true
評価回数についての注意
Optim.f_calls(res)は正確な呼び出し回数をカウントしているわけではないので注意してください.
# using Optim
function f(x)
global calls += 1
return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end
x0 = [0.0, 0.0]
calls = 0
res = optimize(f, x0, method=NelderMead())
println("手法 ：", Optim.summary(res))
println("最小値 ：", Optim.minimum(res))
println("最小点 ：", Optim.minimizer(res))
println("反復回数：", Optim.iterations(res))
println("評価総数：", Optim.f_calls(res))
println("実際の数：", calls)
calls = 0
res = optimize(f, x0, method=LBFGS())
println("手法 ：", Optim.summary(res))
println("最小値 ：", Optim.minimum(res))
println("最小点 ：", Optim.minimizer(res))
println("反復回数：", Optim.iterations(res))
println("評価総数：", Optim.f_calls(res))
println("実際の数：", calls)
手法 ：Nelder-Mead
最小値 ：3.5255270584829996e-9
最小点 ：[0.9999634355313174, 0.9999315506115275]
反復回数：60
評価総数：118
実際の数：118
手法 ：L-BFGS
最小値 ：5.378388330692143e-17
最小点 ：[0.9999999926662504, 0.9999999853325008]
反復回数：24
評価総数：67
実際の数：335
１変数関数
デフォルトのNelder-Mead法は2変数以上を前提とした手法なので, 1変数関数を扱う場合は別の手法を指定する必要があります. 1変数であったとしても, 関数の引数や初期条件を配列として定義します.
f(x) = (x[1] + 2.0)^2 + 1
x0 = [5.0]
optimize(f, x0, method=GradientDescent())
* Status: success
* Candidate solution
Minimizer: [-2.00e+00]
Minimum: 1.000000e+00
* Found with
Algorithm: Gradient Descent
Initial Point: [5.00e+00]
* Convergence measures
|x - x'| = 7.00e+00 ≰ 0.0e+00
|x - x'|/|x'| = 3.50e+00 ≰ 0.0e+00
|f(x) - f(x')| = 4.90e+01 ≰ 0.0e+00
|f(x) - f(x')|/|f(x')| = 4.90e+01 ≰ 0.0e+00
|g(x)| = 9.17e-11 ≤ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 1
f(x) calls: 3
∇f(x) calls: 3
収束の様子
可視化にはPlots.jlを使用します. store_trace=trueを指定すると, Optim.f_trace()で反復毎の関数値（評価毎の関数値ではない）を格納した配列が受け取れます. なお, シンプレックスの履歴はOptim.simplex_trace()で取り出せますが, バグがあり, 現在対応中のようです.
f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
x0 = [0.0, 0.0]
res = optimize(f, x0, store_trace=true, extended_trace=true, method=ConjugateGradient())
Y = Optim.f_trace(res)
plot(Y, label="", xlabel="number of iteration", ylabel="f(x,y)") |> display
X = [x[1] for x in Optim.x_trace(res)]
Y = [x[2] for x in Optim.x_trace(res)]
heatmap(-1.5:0.01:1.5, -0.5:0.01:2.0, (x,y)->f([x,y]), c=:thermal, clim=(0,100), zscale=:log10, xlabel="x", ylabel="y")
plot!(X, Y, label="", lc=:red)
scatter!(X[end:end], Y[end:end], label="", mws=0, mc=:red) |> display
収束の速さ
Optim.iterations()を用いると反復回数, Optim.f_calls()を用いると評価回数を取り出すことができます. 横軸に反復回数を取ると, LBFGS()はNelderMead()より速く収束することが確かめられます. ただし, 関数の評価回数としてはLBFGS()の方が多いので注意してください.
# using Optim
f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
x0 = [0.0, 0.0]
function myOptimize(f, x0, method)
res = optimize(f, x0, method=method, store_trace=true)
println("手法 ：", Optim.summary(res))
println("最小値 ：", Optim.minimum(res))
println("最小点 ：", Optim.minimizer(res))
println("反復回数：", Optim.iterations(res))
println("評価総数：", Optim.f_calls(res))
return plot!(Optim.f_trace(res), label=Optim.summary(res))
end
plot(xlabel="number of iteration", ylabel="f(x,y)")
myOptimize(f, x0, NelderMead())
myOptimize(f, x0, LBFGS())
手法 ：Nelder-Mead
最小値 ：3.5255270584829996e-9
最小点 ：[0.9999634355313174, 0.9999315506115275]
反復回数：60
評価総数：118
手法 ：L-BFGS
最小値 ：5.378388330692143e-17
最小点 ：[0.9999999926662504, 0.9999999853325008]
反復回数：24
評価総数：67
関数の評価回数を自前の関数で数えて, 評価回数を横軸に, 評価値を縦軸に取るとNelder-Mead法が一番早く収束していることがわかる(デフォルト設定でRosenbrock関数に適用した場合).
function myOptimize(f,x0,method;lw=1)
function ff(x)
y = f(x)
push!(Y,y)
return y
end
global Y = []
res = optimize(ff, x0, method=method, store_trace=true, iterations=10000)
println("手法 ：", Optim.summary(res))
println("最小値 ：", Optim.minimum(res))
println("最小点 ：", Optim.minimizer(res))
println("反復回数：", Optim.iterations(res))
println("評価総数：", Optim.f_calls(res))
println("収束判定：", Optim.converged(res))
println()
#plot!(Optim.f_trace(res), label=Optim.summary(res))
plot!(Y, label="$(Optim.summary(res)) : $(Optim.iterations(res)==10000 ? "<" : "")$(length(Y))", lw=lw)
end
f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
x0 = [0.0, 0.0]
Y = []
plot(xlabel="Number of Evaluation", ylabel="\$f(x,y)\$", xlim=(0,200), ylim=(-0.04,1.54))
myOptimize(f,x0,ConjugateGradient())
myOptimize(f,x0,GradientDescent())
myOptimize(f,x0,LBFGS())
myOptimize(f,x0,SimulatedAnnealing())
myOptimize(f,x0,NelderMead(),lw=3)
手法 ：Conjugate Gradient
最小値 ：5.4324993781297417e-17
最小点 ：[0.9999999926294509, 0.9999999852589598]
反復回数：21
評価総数：59
収束判定：true
手法 ：Gradient Descent
最小値 ：1.7140627552964986e-13
最小点 ：[0.9999995872120353, 0.9999991712424621]
反復回数：10000
評価総数：25314
収束判定：false
手法 ：L-BFGS
最小値 ：5.378388330692143e-17
最小点 ：[0.9999999926662504, 0.9999999853325008]
反復回数：24
評価総数：67
収束判定：true
手法 ：Simulated Annealing
最小値 ：0.0026031581168682403
最小点 ：[1.0359966119697626, 1.069673178519981]
反復回数：10000
評価総数：10001
収束判定：false
手法 ：Nelder-Mead
最小値 ：3.5255270584829996e-9
最小点 ：[0.9999634355313174, 0.9999315506115275]
反復回数：60
評価総数：118
収束判定：true
閾値の指定
収束判定の閾値は手法毎に異なります. Nelder-Mead法の場合はg_tolで指定され, 判定に使われる値の記録はOptim.g_norm_trace(res)で取り出せます. 片対数グラフにすると収束がわかりやすいです.
# using Optim
f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
x0 = [0.0, 0.0]
res = optimize(f, x0, store_trace=true, g_tol=1e-8)
Y = Optim.g_norm_trace(res)
plot(Y, label="trace", xlabel="number of iteration", ylabel="", yscale=:log10)
plot!(x->1e-8, label="g_tol")
引用方法
Optim.jlのソースコードはGitHubにあります. 研究で使うときはこの論文を引用しましょう. bibtexはGitHubのReadmeにあります.
動作環境
versioninfo()
Julia Version 1.6.2
Commit 1b93d53fc4 (2021-07-14 15:36 UTC)
Platform Info:
OS: Windows (x86_64-w64-mingw32)
CPU: Intel(R) Core(TM) i7-4650U CPU @ 1.70GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-11.0.1 (ORCJIT, haswell)
参考文献
この記事の元になったJupyter Notebookのデータは下記のリンクにあります.
Discussion