このページは、nAGライブラリのJupyterノートブックExampleの日本語翻訳版です。オリジナルのノートブックはインタラクティブに操作することができます。
一般化されたRosenbrock関数の最小化
境界制約付き最適化を使用して「一般化されたRosenbrock」関数を最小化したいとします。
nAGライブラリへのインターフェースはnaginterfaces.libraryサブパッケージで提供されています
https://www.nag.com/numeric/py/nagdoc_latest/naginterfaces.library.htmlのHTMLドキュメントから、(局所)最適化に関連するアルゴリズムのサブモジュールがoptであることがわかります。
opt機能インデックスを調べると、呼び出すべき関連する最適化ソルバーがbounds_quasi_func_easyであることが確認できます。このソルバーのHTMLドキュメントは
https://www.nag.com/numeric/py/nagdoc_latest/naginterfaces.library.opt.bounds_quasi_func_easy.html
にあります。
必要に応じて、最適化ソルバーを直接インポートすることができます
from naginterfaces.library.opt import bounds_quasi_func_easyここで最適化問題を定義します：まず、一般化されたRosenbrock問題の目的関数です。nAGルーチンの文書化された署名から
help(bounds_quasi_func_easy)naginterfaces.library.optモジュールのbounds_quasi_func_easy関数に関するヘルプ:
bounds_quasi_func_easy(ibound, funct1, bl, bu, x, liw=None, lw=None, data=None) 関数値のみを使用する準ニュートン法アルゴリズムによる、境界制約付き最小化（使いやすいバージョン）。
`bounds_quasi_func_easy`は、独立変数x_1,x_2,...,x_nの固定された上限と下限の制約の下で、
関数F(x_1,x_2,...,x_n)の最小値を見つけるための、使いやすい準ニュートン法アルゴリズムです。
関数値のみを使用します。
この関数は、連続で、連続な一次および二次導関数を持つ関数を対象としています
（ただし、導関数に時折不連続性がある場合でも通常は機能します）。
詳細については、e04jyに関するnAGライブラリドキュメントを参照してください
https://www.nag.com/numeric/nl/nagdoc_27.1/flhtml/e04/e04jyf.html
パラメータ
----------
ibound : int
特殊な形式の境界を扱う機能を使用するかどうかを示します。
以下のいずれかの値に設定する必要があります：
`ibound` = 0
すべてのl_jとu_jを個別に提供する場合。
`ibound` = 1
x_jに境界がない場合。
`ibound` = 2
すべての境界が0 <= x_jの形式の場合。
`ibound` = 3
l_1 = l_2 = ... = l_nかつu_1 = u_2 = ... = u_nの場合。
funct1 : callable fc = funct1(xc, data=None)
任意の点xにおける関数F(x)の値を計算するために`funct1`を提供する必要があります。
`bounds_quasi_func_easy`で使用する前に、個別にテストする必要があります
（[E04イントロダクション]を参照）。
パラメータ
~~~~~~~~~~
xc : float, ndarray, shape (n)
関数値が必要な点x。
data : 任意, オプション, その場で変更可能
コールバック関数用のユーザー通信データ。
戻り値
~~~~~~~
fc : float
現在の点xにおける関数Fの値。
bl : float, array-like, shape (n)
下限l_j。
`ibound`が0に設定されている場合、j = 1,2,...,nに対して`bl`[j-1]をl_jに設定する必要があります。
（特定のx_jに下限が指定されていない場合、対応する`bl`[j-1]を-10^6に設定する必要があります。）
`ibound`が3に設定されている場合、`bl`[0]をl_1に設定する必要があります；
`bounds_quasi_func_easy`は`bl`の残りの要素を`bl`[0]と等しく設定します。
bu : float, array-like, shape (n)
上限u_j。
`ibound`が0に設定されている場合、j = 1,2,...,nに対して`bu`[j-1]をu_jに設定する必要があります。
（特定のx_jに上限が指定されていない場合、対応する`bu`[j-1]を10^6に設定する必要があります。）
`ibound`が3に設定されている場合、`bu`[0]をu_1に設定する必要があります；
`bounds_quasi_func_easy`は`bu`の残りの要素を`bu`[0]と等しく設定します。
x : float, array-like, shape (n)
j = 1,2,...,nに対して、`x`[j-1]を最小値の位置のj番目の成分の推定値に設定する必要があります。
liw : None または int, オプション注意: この引数がNoneの場合、デフォルト値が使用されます。
以下のように決定されます: n+2。
配列 `iw` の次元。
lw : None または int, オプション
注意: この引数がNoneの場合、デフォルト値が使用されます。
以下のように決定されます: max(n*(n-1)/2+12*n, 13)。
配列 `w` の次元。
data : 任意, オプション
コールバック関数用のユーザー通信データ。
戻り値
-------
bl : float, ndarray, 形状 (n)
``bounds_quasi_func_easy`` で実際に使用される下限。
bu : float, ndarray, 形状 (n)
``bounds_quasi_func_easy`` で実際に使用される上限。
x : float, ndarray, 形状 (n)
計算中に見つかった最小点。したがって、終了時に例外や警告が発生しない場合、
`x`[j-1] は最小値の位置のj番目の成分です。
f : float
`x` に格納された最終点に対応する F(x) の値。
iw : int, ndarray, 形状 (`liw`)
関数が正常に終了するか、`errno` が 3 または 5 の場合、
`iw` の最初の n 要素には、どの変数が現在境界上にあり、
どの変数が自由であるかについての情報が含まれます。
具体的には、x_i が:
- 上限に固定されている場合、`iw`[i-1] は -1 です。
- 下限に固定されている場合、`iw`[i-1] は -2 です。
- 実質的に定数の場合（つまり、l_j = u_j）、`iw`[i-1] は -3 です。
- 自由な場合、`iw`[i-1] は自由変数のシーケンスにおけるその位置を示します。
さらに、`iw`[n] には自由変数の数（つまり、n_z）が含まれます。
配列の残りの部分はワークスペースとして使用されます。
w : float, ndarray, 形状 (`lw`)
関数が正常に終了するか、`errno` が 3 または 5 の場合、`w`[i-1] には
射影された勾配ベクトル g_z の i 番目の要素の有限差分近似が含まれます（i = 1,2,...,n）。
さらに、`w`[n] には射影されたヘッセ行列の条件数の推定値（つまり、k）が含まれます。
配列の残りの部分はワークスペースとして使用されます。
例外
------
NagValueError
(`errno` 1)
入力時、n = *<値>*。
制約: n >= 1。
(`errno` 1)
入力時、`ibound` = *<値>*。
制約: 0 <= `ibound` <= 3。
(`errno` 1)
入力時、`lw` = *<値>*。
制約: `lw` >= *<値>*。
(`errno` 1)
入力時、`liw` = *<値>*。
制約: `liw` >= *<値>*。
(`errno` 1)
入力時、`ibound` = 0 で一部の j について `bl`[j-1] > `bu`[j-1]、
または `ibound` = 3 で `bl`[0] > `bu`[0]。
(`errno` 2)
`funct1` への呼び出しが 400*n 回行われました。
(`errno` 4)
計算中にオーバーフローが発生しました。
(`errno` 9)
変数の絶対値が非常に大きくなりました。`funct1` にミスがあるか、
問題に有限の解がないか、問題のスケール変更が必要かもしれません。
(`errno` 10)前方差分の1つが負でした。
警告
-----
NagAlgorithmicWarning
(`errno` 3)
最小値の条件がすべて満たされていませんが、
より低い点が見つからず、アルゴリズムが失敗しました。
(`errno` 5)
局所最小値が見つかった可能性が高いですが、
保証することはできません。
(`errno` 6)
局所最小値が見つかった可能性が高いですが、
保証することはできません。
(`errno` 7)
局所最小値が見つかった可能性は低いです。
(`errno` 8)
局所最小値が見つかった可能性は非常に低いです。
注意
-----
このルーチンに相当する従来のCインターフェースは
nAGライブラリには存在しません。
``bounds_quasi_func_easy``は以下の形式の問題に適用可能です：
F(x_1,x_2,...,x_n)を最小化する。ただし、l_j <= x_j <= u_j, j =
1,2,...,nの制約条件がある
F(x)の導関数が利用できない場合。
x_jに実際に制限がない問題、非負制約のみの問題、
l_1 = l_2 = ... = l_nかつu_1 = u_2 = ... = u_nの問題に対して
特別な対応がなされています。
任意の点xにおけるF(x)の値を計算する関数を提供する必要があります。
あなたが提供した開始点から、F(x)の勾配と曲率の推定に基づいて、
制約付き関数の局所最小値に収束することを意図した
実行可能な点の列が生成されます。
最終点が最小値であることを検証する試みがなされます。
典型的な反復は、現在の点xから始まります。ここでn_z個（仮に）の
変数が両方の境界から自由です。
自由変数に関するF(x)の導関数の有限差分近似である
射影勾配ベクトルg_zが既知です。
単位下三角行列Lと対角行列D（どちらもn_z次元）も保持されており、
LDL^Tは自由変数に関する2次導関数行列（つまり、射影ヘッセ行列）の
正定値近似です。
以下の方程式を解いて探索方向p_zを求めます：
LDL^Tp_z = -g_z
p_zは適切な0要素を挿入してn次元ベクトルpに拡張されます。
次に、F(x+alpha p)がalphaに関して（固定境界に従って）
おおよそ最小になるようなalphaを見つけます。xはx+alpha pに置き換えられ、
行列LとDは、ステップalpha pによって推定勾配に生じた変化と
一致するように更新されます。
pに沿った探索中に実際に変数が境界に達した場合、
その変数は固定され、次の反復でn_zが減少します。
ほとんどの反復ではg_zを前方差分で計算しますが、
必要と思われる場合は中心差分を使用します。
2セットの収束基準があります - 弱いものと強いものです。
弱い基準が満たされるたびに、すべてのアクティブな制約に対して
ラグランジュ乗数が推定されます。
ラグランジュ乗数の推定値が著しく負の場合、
負のラグランジュ乗数に関連する変数の1つが推定値は境界から解放され、次の探索方向が
拡張された部分空間で計算されます（つまり、n_zが増加します）。
そうでない場合、実行可能である限り、現在の部分空間での最小化が続行されます。
それが実行不可能な場合、または強力な収束基準がすでに
満たされている場合、1つ以上のラグランジュ乗数の
推定値がゼロに近い場合、対応する変数の値に
わずかな摂動が加えられ、より低い関数値が
得られるまで順次行われます。
その後、摂動を加えた点から通常のアルゴリズムが再開されます。
鞍点が疑われる場合、鞍点から
離れるための局所探索が実行されます。
制約付き最小値と考えられる点が見つかった場合にも
局所探索が実行されます。
参考文献
----------
Gill, P E and Murray, W, 1976, `変数に境界がある場合の最小化`,
NPL Report NAC 72, National Physical Laboratoryパラメータ funct1 は、関数に渡す通信データ
（data）がない場合、 ラムダ式として指定される可能性があると
推測できます。
rosen = lambda x: (sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1.0-x[:-1])**2.0))以下はRosenbrock関数の2次元等高線プロットです
# Jupyter の表示バックエンドを選択:
%matplotlib nbaggimport matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import numpy as np
delta = 0.05
x = np.arange(-3., 3., delta)
y = np.arange(-3., 3, delta)
X, Y = np.meshgrid(x, y)
Z = np.empty((len(X), len(Y)))
for i, x_i in enumerate(x):
for j, y_j in enumerate(y):
Z[j, i] = rosen(np.array([x_i, y_j]))
fig = plt.figure()
ax = Axes3D(fig, auto_add_to_figure=False)
ax.grid(False)
ax.plot_wireframe(X, Y, Z, color='red', linewidths=0.4)
ax.contour(X, Y, Z, levels=[5, 25, 50, 100, 250, 500, 1000, 1500, 2000, 2500], offset=-3000.0, cmap=cm.jet)
ax.set_title('The 2D Rosenbrock Function')
ax.set_xlabel(r'$\mathit{x}$')
ax.set_ylabel(r'$\mathit{y}$')
ax.set_zlim3d(-1.0, 10000.0)
ax.azim = -20
ax.elev = 20
fig.add_axes(ax)
plt.show()
2次元のRosenbrock関数の最小値を見つけることは比較的簡単なので、10次元のRosenbrock関数の最小値を見つけることにします（これはプロットできません！）。
まず、最適化の初期推測を定義します。
naginterfaces.library サブパッケージでは、上記の引数
x の float, array-like, shape (n)
仕様で述べたように、入力配列データは任意の ‘配列のような’
コンテナで提供できます。1次元の x の場合、
これはデータのシーケンスであれば適切なコンテナとなり、
numpy.ndarray も同様です。
（多次元データを扱う関数では、ネストされたシーケンスや
numpy.ndarray のインスタンスも有効です。）さらに、 提供する
x の形状（長さ）によって、問題の n
の（推論された）値が決まります。
選択した開始点は \((0., 0., ..., 0.,
0.)\) であり、したがって ‘配列のような’ ベクトル x
を提供するには以下のいずれかを使用できます：
listとして
x = [0.]*10- 「タプル」として
x = (0., 0., 0., 0., 0., 0., 0., 0., 0., 0.)- 「ndarray」として
x = np.array([0.]*10)次に問題のボックス制約を定義します
n = len(x)
bl, bu = [0.0]*n, [2.0]*n
ibound = 0問題を最小化する
opt_soln = bounds_quasi_func_easy(ibound, rosen, bl, bu, x)nAGルーチンによって返される引数は、戻り値タプルの名前付きフィールドとしてアクセスできます。
結果を表示する
print('Function value at lowest point found is {:.5f}.'.format(opt_soln.f))
print('The corresponding x is (' + ', '.join(['{:.4f}'] * n).format(*opt_soln.x) + ').')見つかった最低点での関数値は0.00000です。 対応するxは(1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000)です。