SciPy、最適化

SciPy (サイパイと発音) は、Numpy Python 拡張機能に基づく数学アプリケーション パッケージです。 SciPy を使用すると、対話型 Python セッションが、MATLAB、IDL、Octave、R-Lab、SciLab と同じ完全なデータ サイエンスおよび複雑なシステム プロトタイピング環境になります。 今日は、scipy.optimize パッケージでよく知られている最適化アルゴリズムを使用する方法について簡単に説明したいと思います。 関数の使用に関するより詳細な最新のヘルプは、help() コマンドまたは Shift+Tab を使用するといつでも取得できます。

導入

あなた自身と読者が一次情報源を検索したり読んだりする手間を省くために、メソッドの説明へのリンクは主に Wikipedia に掲載されます。 原則として、この情報は、メソッドの一般的な用語とその適用条件を理解するのに十分です。 数学的手法の本質を理解するには、各記事の最後またはお気に入りの検索エンジンにある、より信頼できる出版物へのリンクをたどってください。

したがって、scipy.optimize モジュールには次のプロシージャの実装が含まれています。

1. さまざまなアルゴリズム (ネルダー・ミード単体、BFGS、ニュートン共役勾配、 コビラ и SLSQP)
2. グローバル最適化 (例: 盆地巡り, diff_evolution)
3. 残差の最小化 MNC (least_squares) および非線形最小二乗法を使用した曲線近似アルゴリズム (curve_fit)
4. XNUMX つの変数のスカラー関数の最小化 (minim_scalar) と根の検索 (root_scalar)
5. さまざまなアルゴリズム (ハイブリッド パウエル、 レーベンバーグ＝マルカルト または次のような大規模な方法 ニュートン・クリロフ).

この記事では、このリスト全体から最初の項目のみを取り上げます。

複数の変数のスカラー関数の無条件最小化

scipy.optimize パッケージの minim 関数は、複数の変数のスカラー関数の条件付きおよび無条件の最小化問題を解決するための一般的なインターフェイスを提供します。 それがどのように機能するかを示すには、さまざまな方法で最小化するいくつかの変数の適切な関数が必要です。

これらの目的には、次の形式を持つ N 変数の Rosenbrock 関数が最適です。

Rosenbrock 関数とそのヤコビ行列とヘッセ行列 (それぞれ XNUMX 次導関数と XNUMX 次導関数) はすでに scipy.optimize パッケージで定義されていますが、これを自分たちで定義します。

import numpy as np

def rosen(x):
    """The Rosenbrock function"""
    return np.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0, axis=0)

わかりやすくするために、3 つの変数のローゼンブロック関数の値を XNUMXD で描画してみましょう。

描画コード

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter

# Настраиваем 3D график
fig = plt.figure(figsize=[15, 10])
ax = fig.gca(projection='3d')

# Задаем угол обзора
ax.view_init(45, 30)

# Создаем данные для графика
X = np.arange(-2, 2, 0.1)
Y = np.arange(-1, 3, 0.1)
X, Y = np.meshgrid(X, Y)
Z = rosen(np.array([X,Y]))

# Рисуем поверхность
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm)
plt.show()

最小値が 0 であることを事前に知っておくと、 では、さまざまな scipy.optimize プロシージャを使用して Rosenbrock 関数の最小値を決定する方法の例を見てみましょう。

ネルダー・ミードシンプレックス法

0 次元空間に初期点 x5 があるとします。 アルゴリズムを使用して、それに最も近いローゼンブロック関数の最小点を見つけてみましょう ネルダー・ミードシンプレックス (アルゴリズムはメソッドパラメータの値として指定されます):

from scipy.optimize import minimize
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)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 339
         Function evaluations: 571
[1. 1. 1. 1. 1.]

シンプレックス法は、明示的に定義された非常に滑らかな関数を最小化する最も簡単な方法です。 関数の導関数を計算する必要はなく、その値のみを指定するだけで十分です。 ネルダー・ミード法は、単純な最小化問題には適しています。 ただし、勾配推定を使用しないため、最小値を見つけるのに時間がかかる可能性があります。

パウエル法

関数値のみを計算する別の最適化アルゴリズムは、 パウエル方式。 これを使用するには、minim 関数で method = 'powell' を設定する必要があります。

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

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 19
         Function evaluations: 1622
[1. 1. 1. 1. 1.]

ブロイデン・フレッチャー・ゴールドファーブ・シャンノ (BFGS) アルゴリズム

解への収束を速くするには、次の手順を実行します。 BFGS 目的関数の勾配を使用します。 勾配は関数として指定することも、一次差分を使用して計算することもできます。 いずれの場合も、BFGS メソッドは通常、シンプレックス メソッドよりも必要な関数呼び出しの数が少なくなります。

Rosenbrock 関数の導関数を解析形式で求めてみましょう。

この式は、最初と最後の変数を除くすべての変数の導関数に対して有効であり、次のように定義されます。

この勾配を計算する Python 関数を見てみましょう。

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

勾配計算関数は、以下のようにminim関数のjacパラメータの値として指定します。

res = minimize(rosen, x0, method='BFGS', jac=rosen_der, options={'disp': True})
print(res.x)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 25
         Function evaluations: 30
         Gradient evaluations: 30
[1.00000004 1.0000001  1.00000021 1.00000044 1.00000092]

共役勾配アルゴリズム (Newton)

アルゴリズム ニュートンの共役勾配 修正ニュートン法です。
ニュートンの方法は、局所領域の関数を XNUMX 次の多項式で近似することに基づいています。

どこ は XNUMX 次導関数の行列 (ヘッセ行列、ヘッセ行列) です。
ヘッセ行列が正定値の場合、二次形式のゼロ勾配をゼロとみなすことによって、この関数の極小値を見つけることができます。 結果は次の式になります。

逆ヘッセ行列は、共役勾配法を使用して計算されます。 この方法を使用してローゼンブロック関数を最小化する例を以下に示します。 Newton-CG 法を使用するには、ヘッセ行列を計算する関数を指定する必要があります。
分析形式のローゼンブロック関数のヘッセ行列は次と等しくなります。

どこ и 、行列を定義します .

行列の残りの非ゼロ要素は次と等しくなります。

たとえば、5 次元空間 N = XNUMX では、Rosenbrock 関数のヘッセ行列はバンドの形式になります。

このヘッセ行列を計算するコードと、共役勾配 (ニュートン) 法を使用してローゼンブロック関数を最小化するコード:

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)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 24
         Function evaluations: 33
         Gradient evaluations: 56
         Hessian evaluations: 24
[1.         1.         1.         0.99999999 0.99999999]

ヘッセ行列と任意のベクトルの積関数の定義の例

現実の問題では、ヘッセ行列全体を計算して保存するには、かなりの時間とメモリ リソースが必要になる場合があります。 この場合、実際にはヘッセ行列自体を指定する必要はありません。 最小化手順には、ヘッセ行列と別の任意のベクトルの積に等しいベクトルのみが必要です。 したがって、計算の観点からは、ヘッセ行列と任意のベクトルの積の結果を返す関数をすぐに定義することが非常に望ましいです。

hess 関数を考えてみましょう。この関数は、最小化ベクトルを最初の引数として、任意のベクトルを XNUMX 番目の引数として (最小化される関数の他の引数とともに) 取ります。 私たちの場合、Rosenbrock 関数のヘッセ行列と任意のベクトルの積を計算することはそれほど難しくありません。 もし p が任意のベクトルの場合、積は の形式は次のとおりです。

ヘッセ行列と任意のベクトルの積を計算する関数は、hessp 引数の値として minimum 関数に渡されます。

def rosen_hess_p(x, p):
    x = np.asarray(x)
    Hp = np.zeros_like(x)
    Hp[0] = (1200*x[0]**2 - 400*x[1] + 2)*p[0] - 400*x[0]*p[1]
    Hp[1:-1] = -400*x[:-2]*p[:-2]+(202+1200*x[1:-1]**2-400*x[2:])*p[1:-1] 
    -400*x[1:-1]*p[2:]
    Hp[-1] = -400*x[-2]*p[-2] + 200*p[-1]
    return Hp

res = minimize(rosen, x0, method='Newton-CG',
               jac=rosen_der, hessp=rosen_hess_p,
               options={'xtol': 1e-8, 'disp': True})

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 24
         Function evaluations: 33
         Gradient evaluations: 56
         Hessian evaluations: 66

共役勾配信頼領域アルゴリズム (Newton)

ヘッセ行列の条件付けが不十分で検索方向が間違っていると、ニュートンの共役勾配アルゴリズムが無効になる可能性があります。 このような場合には、以下が優先されます。 トラストリージョン方式 (trust-region) 共役ニュートン勾配。

ヘッセ行列の定義の例:

res = minimize(rosen, x0, method='trust-ncg',
               jac=rosen_der, hess=rosen_hess,
               options={'gtol': 1e-8, 'disp': True})
print(res.x)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 20
         Function evaluations: 21
         Gradient evaluations: 20
         Hessian evaluations: 19
[1. 1. 1. 1. 1.]

ヘッセ行列と任意のベクトルの積関数の例:

res = minimize(rosen, x0, method='trust-ncg', 
                jac=rosen_der, hessp=rosen_hess_p, 
                options={'gtol': 1e-8, 'disp': True})
print(res.x)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 20
         Function evaluations: 21
         Gradient evaluations: 20
         Hessian evaluations: 0
[1. 1. 1. 1. 1.]

クリロフ型メソッド

trust-ncg 手法と同様、クリロフ型手法は行列ベクトルの積のみを使用するため、大規模な問題を解決するのに適しています。 それらの本質は、切り詰められたクリロフ部分空間によって制限された信頼領域内で問題を解決することです。 不確実な問題の場合は、trust-ncg 法と比較して部分問題ごとの行列ベクトル積の数が少ないため、使用する非線形反復の回数が少なくなるこの方法を使用することをお勧めします。 さらに、二次部分問題の解は、trust-ncg 法を使用するよりも正確に見つかります。
ヘッセ行列の定義の例:

res = minimize(rosen, x0, method='trust-krylov',
               jac=rosen_der, hess=rosen_hess,
               options={'gtol': 1e-8, 'disp': True})

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 19
         Function evaluations: 20
         Gradient evaluations: 20
         Hessian evaluations: 18

print(res.x)

    [1. 1. 1. 1. 1.]

ヘッセ行列と任意のベクトルの積関数の例:

res = minimize(rosen, x0, method='trust-krylov',
               jac=rosen_der, hessp=rosen_hess_p,
               options={'gtol': 1e-8, 'disp': True})

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 19
         Function evaluations: 20
         Gradient evaluations: 20
         Hessian evaluations: 0

print(res.x)

    [1. 1. 1. 1. 1.]

信頼領域における近似解のアルゴリズム

すべての手法 (Newton-CG、trust-ncg、trust-krylov) は、大規模な問題 (数千の変数を含む) を解決するのに適しています。 これは、基礎となる共役勾配アルゴリズムが逆ヘッセ行列の近似決定を意味するという事実によるものです。 解は、ヘッセ行列を明示的に展開することなく、反復的に求められます。 ヘッセ行列と任意のベクトルの積の関数を定義するだけでよいため、このアルゴリズムはスパース (バンド対角) 行列を扱うのに特に適しています。 これにより、メモリのコストが低くなり、時間が大幅に節約されます。

中規模の問題の場合、ヘッセ行列の保存と因数分解のコストは重要ではありません。 これは、より少ない反復で解を取得し、信頼領域の部分問題をほぼ正確に解決できることを意味します。 これを行うには、二次部分問題ごとにいくつかの非線形方程式を反復的に解きます。 このような解決策には、通常、ヘッセ行列の 3 つまたは 4 つのコレスキー分解が必要です。 その結果、この方法はより少ない反復回数で収束し、他の実装された信頼領域方法よりも必要な目的関数の計算が少なくなります。 このアルゴリズムには完全なヘッセ行列の決定のみが含まれており、ヘッセ行列と任意のベクトルの積関数を使用する機能はサポートされていません。

Rosenbrock 関数を最小化した例:

res = minimize(rosen, x0, method='trust-exact',
               jac=rosen_der, hess=rosen_hess,
               options={'gtol': 1e-8, 'disp': True})
res.x

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 13
         Function evaluations: 14
         Gradient evaluations: 13
         Hessian evaluations: 14

array([1., 1., 1., 1., 1.])

おそらくそこで止まるでしょう。 次の記事では、条件付き最小化、近似問題を解く際の最小化の応用、XNUMX 変数の関数の最小化、任意のミニマイザー、および scipy.optimize を使用した方程式系の根の検索について、最も興味深いことを説明しようと思います。パッケージ。

出所： https://docs.scipy.org/doc/scipy/reference/

出所： habr.com