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次導関数を使うアルゴリズムに関しては、パラメータの数が増えてくると、その数の2乗で計算量が増えるヘッセ行列の計算に時間がかかってしまうこともあり、勾配ベクトルのみを使うアルゴリズムを選択すべきこともあります。また、導関数を解析的に導出するのが面倒なことも多いので、あえて1次導関数を使うものにとどめることもあります。

さて、手始めに、パラメータが 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])

どちらも、x≈0が得られましたね。

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))のグラフを書くと次のようになります。

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