SciPyによる最適化(scipy.optimize)①|SciPy入門 #5

f:id:lib-arts:20190225182115g:plain

SciPyについて色々と話題になり面白そうだったので公式チュートリアルを元にまとめています。

SciPy Tutorial — SciPy v1.2.1 Reference Guide
#3、#4ではscipy.statsを元に確率分布や検定に関して取り扱いました。

 #4でscipy.statsのチュートリアルは一通り終了したので、#5からは最適化に関する機能であるscipy.optimizeについて取り扱っていきます。

Optimization (scipy.optimize) — SciPy v1.2.1 Reference Guide
以下目次になります。

1. scipy.optimizeの概要(冒頭部)
2. Unconstrained minimization of multivariate scalar functions (minimize)
2.1 Nelder-Mead Simplex algorithm
2.2 Broyden-Fletcher-Goldfarb-Shanno algorithm
2.3 Newton-Conjugate-Gradient algorithm
3. まとめ

 

1. scipy.optimizeの概要(冒頭部)
冒頭部にはscipy.optimizeの概要がまとまっており、scipy.optimizeはいくつかの最適化でよく用いられるメソッドが実装されているとされています。

f:id:lib-arts:20190301114752p:plain
上記の5点がまとまっており、中でも良く用いるのが3の二乗誤差最小化(Least-squares minimization)と曲線フィッティング(curve fitting)なのではないかと思います。
とはいえ、機能も多いため#5では曲線フィッティングまでは取り扱わず、多変数関数の制約なしの最適化(最小化)までを取り扱えればと思います。2節ではこちらについて諸々まとめたいと思います。


2. Unconstrained minimization of multivariate scalar functions (minimize)
scipy.optimize.minimizeが制約なし/制約ありの双方の最適化アルゴリズムに関するインターフェースを提供しているとされています。問題を具体的に取り扱うにあたって、Rosenbrock functionの最小化が例として取り上げられています。

f:id:lib-arts:20190301120042p:plain
実際に最小化するRosenbrock functionは上記のように立式されています。この問題を実際にscipy.optimizeを用いて解いていきます。


2.1 Nelder-Mead Simplex algorithm
まず最初にNelder-Mead simplex algorithmについて解説されています。

import numpy as np
from scipy.optimize import minimize
def rosen(x):
    """The Rosenbrock function"""
    return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)

x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
res = minimize(rosen, x0, method='nelder-mead',
                options={'xtol': 1e-8, 'disp': True})

print(res.x)

上記を実行すると最適解として[1. 1. 1. 1. 1.]を得ることができます。scipy.minimizeの使用に関しては、引数に(1)最適化対象の関数、(2)初期値を与え、あとはメソッドとオプションを指定して実行する形式になっています。結果はres.xで得ることができます。


2.2 Broyden-Fletcher-Goldfarb-Shanno algorithm
次にNelder-Mead Simplexよりも高速な手法として、BFGS(Broyden-Fletcher-Goldfarb-Shanno)アルゴリズムについて解説されています。BFGSでは一次微分(first-differences)を用いて解の導出を行います。

def rosen_der(x):
    xm = x[1:-1]
    xm_m1 = x[:-2]
    xm_p1 = x[2:]
    der = np.zeros_like(x)
    der[1:-1] = 200*(xm-xm_m1**2) - 400*(xm_p1 - xm**2)*xm - 2*(1-xm)
    der[0] = -400*x[0]*(x[1]-x[0]**2) - 2*(1-x[0])
    der[-1] = 200*(x[-1]-x[-2]**2)
    return der
res = minimize(rosen, x0, method='BFGS', jac=rosen_der,
                    options={'disp': True})
print(res.x)

上記において、関数rosen_derは元のRosenbrock functionの微分にあたるのですが、scipy.optimize.minimizeのメソッドないではjacというパラメータで関数の一次微分の数式を取り扱います。また、2.1同様関数自体はminimizeを使用して今すが、methodでBFGSを指定しているところが違います。
上記を実行することで、最適化を実行し実行結果を得ることができます。


2.3 Newton-Conjugate-Gradient algorithm
次にニュートン-CG法(Newton-Conjugate-Gradient algorithm)についてまとめられています。

def rosen_hess(x):
    x = np.asarray(x)
    H = np.diag(-400*x[:-1],1) - np.diag(400*x[:-1],-1)
    diagonal = np.zeros_like(x)
    diagonal[0] = 1200*x[0]**2-400*x[1]+2
    diagonal[-1] = 200
    diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:]
    H = H + np.diag(diagonal)
    return H
res = minimize(rosen, x0, method='Newton-CG',
                jac=rosen_der, hess=rosen_hess,
                options={'xtol': 1e-8, 'disp': True})
print(res.x)

上記のようにヘッセ行列を作成し、jacとして勾配(rosen_der)、hessとしてヘッセ行列(rosen_hess)を引数として与え、methodに'Newton-CG'を与え実行を行います。
上記を実行することで、最適化を実行し実行結果を得ることができます。


3. まとめ
上記でRosenbrock functionを題材にした制約なしの多変量スカラー関数の最適化について色々と見てきました。
#6では制約があるケースの取り扱いや最小二乗法などを見ていければと思います。