Python

Scipyによる最適化計算と自動微分

Scipyによる最適化計算

Pythonで最適化計算を行う場合、Scipyのoptimizationパッケージを使うことになると思う。これには、いろいろな最適化アルゴリズム(Nelder-Mead, Powell, BFGSなど)が実装されており、アルゴリズムによって、2次導関数(ヘッセ行列)まで必要であったり、あるいは、1次導関数(勾配ベクトル)すら必要なかったり、いろいろである。次に示すのは、制約なし最適化で使えるアルゴリズム。

関数のみが必要:
– Nelder-Mead
– Powell

1次導関数が必要:
– CG
– BFGS

2次導関数まで必要:
– Newton-CG
– dogleg
– trust-ncg
– trust-krylov
– trust-exact

アルゴリズムに与えることのできる情報が多いほど、収束までのステップ数が減る傾向がある。特に、Nelder-MeadとPowellは関数の評価回数がステップごとにパラメータ数倍になるので、解析的な導関数が用意できるなら、それを活用できるアルゴリズムの方が良い。ただし、2次導関数を使うアルゴリズムに関しては、パラメータの数が増えてくると、その数の二乗で計算量が増えるヘッセ行列の計算に時間がかかってしまうこともあり、勾配ベクトルのみを使うアルゴリズムを選択すべきこともある。また、導関数を解析的に導出するのが面倒なことも多い。

さて、手始めに、パラメータが x ひとつだけの、f(x)=x^4の最小化問題を解いてみることにする。解はf(x=0)=0である。

import scipy.optimize
import numpy as np

def f(x):
    return x**4

x0 = np.array([1.0])
scipy.optimize.minimize(f, x0, method='Nelder-Mead')

# 出力
 final_simplex: (array([[-8.8817842e-16], [ 9.7656250e-05]]), array([6.22301528e-61, 9.09494702e-17]))
           fun: 6.223015277861142e-61
       message: 'Optimization terminated successfully.'
          nfev: 34
           nit: 17
        status: 0
       success: True
             x: array([-8.8817842e-16])

簡単に、x≈0が得られた。

せっかくなので、導関数を与える方法も試してみる。jacに1次導関数、つまり、勾配を指定する。

import scipy.optimize
import numpy as np

def f(x):
    return x**4

def dfdx(x):
    return 4*x**3

x0 = np.array([1.0])
scipy.optimize.minimize(f, x0, jac=dfdx, method='BFGS')

# 出力
      fun: 1.0000000000000035e-08
 hess_inv: array([[1]])
      jac: array([-4.e-06])
  message: 'Optimization terminated successfully.'
     nfev: 2
      nit: 1
     njev: 2
   status: 0
  success: True
        x: array([-0.01])

さらに、hessに2次導関数を指定してみる。

import scipy.optimize
import numpy as np

def f(x):
   return x**4

def dfdx(x):
   return 4*x**3

def ddfdxdx(x):
   return 12*x**2

x0 = np.array([1.0])
scipy.optimize.minimize(f, x0, jac=dfdx, hess=ddfdxdx, method='Newton-CG')

# 出力
     fun: array([6.97034909e-10])
     jac: array([5.42626361e-07])
 message: 'Optimization terminated successfully.'
    nfev: 14
    nhev: 14
     nit: 14
    njev: 14
  status: 0
 success: True
       x: array([0.00513823])

Autogradによる自動微分

pip install autogradでインストールできる。GitHubからサンプルを引用すると次のようなスタイルで使える。これを使えば、解析的に導関数を導出しなくても、勾配やヘッセ行列を得ることができる。自動微分の原理については、こちらのビデオを見るとよいらしい(私はまだ見ていないので、あしからず)。

import autograd.numpy as np
from autograd import grad

def tanh(x):
    return (np.exp(x) - np.exp(-x)) / (np.exp(x) + np.exp(-x))

grad(tanh) # -> 0.4199

ここで、いまさらながら複数パラメータでないとヘッセ行列が行列っぽくないと思ったので、複数パラメータの関数を考える。
f(x, y, z) = x^2 + (y^2 + z^2 + 1) * (2 + sin(y))

偏導関数は、解析的に計算できるが、これを明示的に記述せず、Autogradによって計算するようにしてみる。
式から明らかに、x=z=0のとき、(y^2+1)*(2+sin(y))が最小になるときがこれの最小である。

import autograd.numpy as np
from autograd import grad, jacobian

def f(c):
    x, y, z = c
    return x**2 + (y**2 + z**2 + 1) * (2 + np.sin(y))

jac = jacobian(f)
hess = jacobian(jac)

c = np.array([1.0, 1.0, 1.0])
print(jac(c))
print(hess(c))

# 出力
[2.         7.30384889 5.68294197]

[[2.         0.         0.        ]
 [0.         5.31973824 1.08060461]
 [0.         1.08060461 5.68294197]]

このjacとhessを使えば、2次導関数まで必要な、Newton-CG法を次のように実行できる。

import scipy.optimize

c0 = np.array([1.0, 1.0, 1.0])
scipy.optimize.minimize(f, c0, jac=jac, hess=hess, method='Newton-CG')

# 出力
     fun: 1.857815580169501
     jac: array([ 1.31720930e-07, -1.30621411e-08, -3.76504435e-07])
 message: 'Optimization terminated successfully.'
    nfev: 8
    nhev: 7
     nit: 7
    njev: 8
  status: 0
 success: True
       x: array([ 2.66610421e-11, -3.07249594e-01, -3.62583333e-12])

最小値1.8578が求まる。ただし、この関数は極小値をたくさん持つので、今回は初期値の取り方が良かった。実際に、yの初期値を大きくずらすと、見事に極小値にはまる。(y^2+1)*(2+sin(y))のグラフを書くと次のようになる。

大域的最適化など、いくつかやりようはあるものの、初期値を適切に選ぶのは重要で、初期値をその関数に対して何の知識もない状態で、うまく決めるのは難しい。

RECOMMENDこちらの記事も人気です。