Scipyの最適化に導関数をどう与えるか：Autogradとの組み合わせ

はじめに

scipy.optimize.minimizeのmethod引数には10種類を超えるアルゴリズムを指定でき、どれを選んでもたいてい何かしらの答えは返ってきます。ただ、必要とする情報はアルゴリズムごとに違い、関数値だけで探索するものから、勾配やヘッセ行列まで要求するものまであります。この違いが収束までの関数評価回数、つまり計算コストを左右します。以下では、Scipyの最適化手法を導関数の要否で整理したうえで、導関数を手で求めるのが難しい場合にAutogradの自動微分をどう組み合わせるかを実装で確認します。

対象読者:

Pythonを用いて最適化問題を解きたいと考えているエンジニアや研究者
Scipyのoptimizationパッケージの基本的な使い方を学びたい方
自動微分(Autograd)に興味があり、最適化問題への応用を検討している方

記事のポイント:

Scipyのoptimizationパッケージに含まれる主要な最適化アルゴリズムの紹介と使い分け
導関数(1次、2次)の情報を使うと最適化がどう効率化されるか
Autogradによる自動微分の導入と使いどころ
複数パラメータの関数を最適化するときの注意点と初期値の影響

Scipyの最適化アルゴリズム

Scipyのoptimizationパッケージには、さまざまな最適化アルゴリズムが実装されています。いずれも目的関数の情報(値、勾配、ヘッセ行列など)をもとにパラメータを探索しますが、必要とする情報はアルゴリズムごとに異なります。

制約なし最適化で使える主なアルゴリズムを、必要な情報とあわせて整理します。

アルゴリズム	必要情報	特徴
Nelder-Mead	関数値のみ	導関数が不要。パラメータ数が多いと計算コストが高くなる可能性がある。
Powell	関数値のみ	導関数が不要。パラメータ数が多いと計算コストが高くなる可能性がある。
CG	関数値、1次導関数(勾配)	勾配情報を利用。
BFGS	関数値、1次導関数(勾配)	勾配情報を利用。準ニュートン法の一種。
Newton-CG	関数値、1次導関数、2次導関数(ヘッセ行列)	ヘッセ行列を利用。パラメータ数が多いとヘッセ行列の計算コストが高くなる可能性がある。
dogleg	関数値、1次導関数、2次導関数(ヘッセ行列)	ヘッセ行列を利用。信頼領域法の一種。
trust-ncg	関数値、1次導関数、2次導関数(ヘッセ行列)	ヘッセ行列を利用。信頼領域法の一種。
trust-krylov	関数値、1次導関数、2次導関数(ヘッセ行列)	ヘッセ行列を利用。信頼領域法の一種。
trust-exact	関数値、1次導関数、2次導関数(ヘッセ行列)	ヘッセ行列を利用。信頼領域法の一種。

与える情報が多いほど(2次導関数まで使うほど)、収束までのステップ数は少なくなる傾向があります。解析的な導関数が手に入るなら、それを使うアルゴリズムのほうが、少ない関数評価回数で解にたどり着きます。

導関数の利用とアルゴリズムの選択

どれを選ぶかは、計算コストと精度のつり合いで決まります。Nelder-MeadやPowellは導関数なしで動きますが、関数評価回数がステップごとにパラメータ数の数倍になるため、パラメータ数が多いと計算コストがかさみます。

逆にNewton-CGのように2次導関数を使うアルゴリズムでは、パラメータ数が増えるとヘッセ行列の計算そのものが重くなります。結局のところ、パラメータ数と導関数の計算コストを見比べて選ぶことになります。

簡単な例：4次関数の最小化

パラメータが $x$ ひとつだけの関数 $f(x) = x^4$ の最小化問題を例に、各アルゴリズムの動作を確認します。この関数の最小値は $x = 0$ のとき $f(x) = 0$ です。

Nelder-Mead法

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$ が得られました。

BFGS法（1次導関数を使用）

1次導関数(勾配)をjac引数に指定します。

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])

Newton-CG法（2次導関数を使用）

2次導関数をhess引数に指定します。

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による自動微分

複雑な関数になると、導関数を手で導出するのは現実的でなくなります。そこで使えるのが自動微分で、Autogradはそのための代表的なPythonライブラリです。

Autogradのインストール

pip install autograd

Autogradの基本的な使い方

Autogradを使うと、解析的に導関数を計算しなくても、数値的に勾配やヘッセ行列を計算できます。

import autograd.numpy as np
from autograd import grad

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

grad(tanh)(1.0)

# 出力
0.41997434161402614

grad(tanh)は、tanh関数の微分(導関数)を計算する関数を返します。grad(tanh)(1.0)は、 $x=1.0$ におけるtanh関数の微係数を計算します。

複数パラメータ関数の最適化

複数パラメータの関数 $f(x, y, z) = x^2 + (y^2 + z^2 + 1) * (2 + sin(y))$ の最適化を考えます。この関数の偏導関数は解析的に計算可能ですが、ここではAutogradを使って勾配とヘッセ行列を計算します。この関数は、 $x = z = 0$ のとき、 $(y^2 + 1) * (2 + sin(y))$ が最小になり、最小値として $1.8578$ をとります。

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]]

Autogradで計算したjacとhessを使って、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$ が得られました。ただし、この関数は極小値を多数持つため、出発点によっては別の局所最適解に収束します。実際に、初期値を[1.0, 1.0, 1.0]から[1.0, 5.0, 1.0]に変えると、局所最適解に収束してそこから動きません。

大域的最適化アルゴリズム(例えばscipy.optimize.dual_annealing)で局所解を避ける手もありますが、計算コストは大きくなります。多くの場合は、関数の形状について分かっていることを初期値に反映するほうが、少ない計算量で良い解に届きます。

まとめ

Scipyのminimizeは、導関数なしで動くNelder-MeadやPowellから、勾配とヘッセ行列まで使うNewton-CG系まで、同じインターフェースで呼び分けられます。導関数を渡せる問題なら渡したほうが少ない関数評価で収束し、手で導出できない関数でも、Autogradのjacobianを2回かければ勾配とヘッセ行列をjacとhessにそのまま渡せます。どのアルゴリズムも保証するのは局所解までなので、極小値が多い関数では初期値の置き方が結果を左右します。

手元の最適化問題に取り組むときは、まずパラメータ数と導関数の入手しやすさでアルゴリズムを絞り、導関数が書けなければ自動微分で補ってください。そのうえで、関数の形状について知っていることを初期値に反映し、余裕があれば初期値を数点変えて同じ解に収束するかを確かめると、局所解に落ちていないかの確認になります。