勾配計算と最適化問題
初めに
多くの物理の問題は、あるパラメータに依存したコスト関数を考えて、その最大・最小化問題 (もしくは、停留値問題）を解くことに帰着されます。この章では、このような最適化問題の数値解法を学びます。多くの数値的な最小アルゴリズムは、コスト関数の勾配の情報を必要としますので、まずは数値微分の方法について学びます。
以下のコードは適当なプロジェクト環境を作って、その中で作業してみましょう。
プロジェクト環境の準備
Pkgライブラリは使わず、プロジェクト環境を作りたいディレクトリ内で、activateすることで、Project.tomlを作ります。
$ cd ~/workspace
$ mkdir optimizead # 今回のプロジェクト用
$ cd optimizead
$ julia
(@v1.9) pkg> activate .
Activating new project at `~/workspace/optimizead`
(optimizead) pkg> add Optim, Zygote, ForwardDiff, ReverseDiff
Gitレポジトリの初期化とProject.tomlの登録
$ git init
Initialized empty Git repository ~/workspace/optimizead/.git/
$ git add Project.toml
$ git commit -m "Initial commit"
勾配の計算
差分公式
中心差分法を紹介しましょう。関数
と近似されます。ここで、
数
-
h O(h^2) -
h f(\vec{x} + h \vec{e}_i) - f(\vec{x} - h \vec{e}_i)
このように、適切な
using Optim
# 最小化したい関数
f(x::Array{Float64}) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
# 初期値
initial_x = [0.0, 0.0]
# 数値微分（差分法）で勾配を計算
function g!(G::Array{Float64}, x::Array{Float64})
h = 1e-8
for i in 1:length(x)
tmp = x[i]
x[i] = tmp + h
f1 = f(x)
x[i] = tmp - h
f2 = f(x)
G[i] = (f1 - f2) / (2h)
x[i] = tmp
end
end
# 最適化問題を解く
result = optimize(f, g!, initial_x, LBFGS(), Optim.Options(show_trace=true))
println("The minimum is at: ", Optim.minimizer(result))
最適化の過程を表示するためにOptim.Options(show_trace=true)を使っています。
勾配を計算する関数g!の!は、第一引数Gの内容を関数内で変更することが一目で分かるように付けています (Juliaの慣習)。
勾配の計算結果を配列としてreturnするような関数を定義したい場合もあります。
このような方法は直感的に分かりやすいですが、多くの配列をallocateするので、効率性では劣ります。
その場合には、以下の様にかけます (注目点: similar, inplace)。
# 数値微分（差分法）で勾配を計算
function g(x::Array{Float64})
G = similar(x)
h = 1e-8
for i in 1:length(x)
tmp = x[i]
x[i] = tmp + h
f1 = f(x)
x[i] = tmp - h
f2 = f(x)
G[i] = (f1 - f2) / (2h)
x[i] = tmp
end
return G
end
# 最適化問題を解く
result = optimize(f, g, initial_x, LBFGS(), Optim.Options(show_trace=true), inplace=false)
println("The minimum is at: ", Optim.minimizer(result))
実際には、Optimには数値微分（差分法）を使わなくても、勾配を計算する機能があります。
result = optimize(f, initial_x, LBFGS(), Optim.Options(show_trace=true), autodiff = :finite)
演習問題
- 刻み幅
h
自動微分
自動微分の一般論については、Wikipediaの記事が参考になります。Juliaでは自動微分を提供するライブラリが複数存在しています。プロジェクト環境に導入した上で、利用してみましょう。一般に、関数の出力の要素数 < 入力の要素数の場合、トップダウン型自動微分 (リバースモード自動微分とも呼ばれる)が効率的です。
using Optim
# 最小化したい関数 (xの型指定をしていない！)
f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
# 数値微分（差分法）で勾配を計算
function g(x::Array{Float64})
G = similar(x)
h = 1e-8
for i in 1:length(x)
tmp = x[i]
x[i] = tmp + h
f1 = f(x)
x[i] = tmp - h
f2 = f(x)
G[i] = (f1 - f2) / (2h)
x[i] = tmp
end
return G
end
# 初期値
initial_x = [0.0, 0.0]
差分公式で計算した勾配には「大きな」誤差がある。
# Finite difference method
g(initial_x)
2-element Vector{Float64}:
-1.9999999933961732
0.0
一方、自動微分の場合には、ほぼ厳密な値が計算可能。
# Reverse-mode AD/back propagation using Zygote
using Zygote
gradient(f, initial_x)
([-2.0, 0.0],)
自動微分による勾配の計算結果を返す関数g_zygoteを作って、optimizeに渡してみましょう。
g_zygote = x -> gradient(f, x)[1]
g_zygote(initial_x)
2-element Vector{Float64}:
-2.0
0.0
# 最適化問題を解く
result = optimize(f, g_zygote, initial_x, LBFGS(), Optim.Options(show_trace=true), inplace=false)
実際には、forward-mode ADはOptimに内蔵されています。
# もしくは、Optim内蔵のForwardDiffを使う
result = optimize(f, initial_x, LBFGS(), Optim.Options(show_trace=true), autodiff = :forward)
演習問題
- 差分公式と自動微分で計算した勾配の精度を比較しましょう。
- [発展] 最小自乗法を実装してみましょう。
Discussion