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

はじめに

この記事では、PythonのScipyライブラリを用いて最適化問題を解く方法について解説します。Scipyのoptimizationパッケージには多様な最適化アルゴリズムが実装されており、問題の特性に応じて最適なアルゴリズムを選択することが重要です。また、Autogradライブラリを用いた自動微分を利用することで、複雑な関数の最適化も効率的に行うことができます。

この記事を書いたひと

デジタルリアクタ合同会社 代表
望月 優輝 / Yuki Mochizuki   

機械学習・統計、数値計算などの領域を軸としたソフトウェアエンジニアリングを専門としています。スタートアップからグローバル企業まで、さまざまなスケールの企業にて、事業価値に直結する計算システムを設計・構築してきました。主に機械学習の応用分野において、出版・発表、特許取得等の実績があり、また、IT戦略やデータベース技術等の情報処理に関する専門領域の国家資格を複数保有しています。九州大学理学部卒業。日本ITストラテジスト協会正会員。

対象読者:

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

記事のポイント:

• Scipyのoptimizationパッケージに含まれる主要な最適化アルゴリズムの紹介と使い分け
• 導関数(1次、2次)の情報を活用することによる最適化の効率化
• Autogradを用いた自動微分の導入とそのメリット
• 複数パラメータを持つ関数の最適化における注意点と初期値の重要性

Scipyの最適化アルゴリズム

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

以下は、制約なし最適化で利用できるアルゴリズムとその特徴をまとめたテーブルです。

一般的に、アルゴリズムに与える情報が多いほど(例えば、2次導関数まで利用するなど)、収束までのステップ数は減少する傾向があります。解析的な導関数が利用できる場合は、それを利用できるアルゴリズムの方が、より少ない関数評価回数で効率的に最適解を探索できます。

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

アルゴリズムの選択は、計算コストと精度のバランスを考慮して行う必要があります。Nelder-MeadやPowellは導関数を必要としませんが、関数評価回数がステップごとにパラメータ数の数倍になるため、パラメータ数が多い場合には計算コストが高くなる可能性があります。

一方で、Newton-CGのように2次導関数を利用するアルゴリズムは、パラメータ数が増えるとヘッセ行列(2次導関数)の計算に時間がかかってしまうことがあります。そのため、問題の規模や特性(パラメータ数、導関数の計算コストなど)に応じて、適切なアルゴリズムを選択することが重要です。

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

パラメータが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のoptimizationパッケージを用いた最適化計算と、Autogradによる自動微分について解説しました。

• 問題の特性(パラメータ数、導関数の計算コストなど)に応じて、Scipyのoptimizationパッケージから適切なアルゴリズムを選択する。
• 解析的に導関数が計算できる場合は、それを利用することで最適化の効率が向上する。
• 導関数の計算が困難な場合は、Autogradなどの自動微分ライブラリを利用する。
• 初期値の設定は最適化の結果に大きく影響するため、関数の形状に関する事前知識に基づいて慎重に決定する。

これらの点を踏まえ、ScipyとAutogradを効果的に活用することで、さまざまな最適化問題を効率的に解くことができるようになります。

関連投稿:

確率的最適化アルゴリズムの実装と比較：遺伝的アルゴリズムと勾配降下法 3次元回転の最適化計算 Graph Neural Networksを用いた異常検知入門 任意の確率密度関数を持つ乱数の生成