optimize-simplex.py 技術ドキュメント

プログラム optimize-simplex.py は、Nelder-Mead法（シンプレックス法）を用いて多変数関数の最小値を探索するPythonプログラムです。ただし、現在のコードにはPythonの構文エラーや潜在的な実行時エラーが含まれており、そのままでは期待通りに動作しません。本ドキュメントでは、コードの意図されている機能と、発生しうる問題を詳細に説明します。

プログラムの動作

このプログラムの主な目的は、与えられた多変数関数の最小値を見つけることです。具体的には、目的関数 func(x) がコード内で定義されており、その関数値を最小にするような変数 x の値を探索します。

現在の func(x) は $f(x_0, x_1) = (x_0 - 1.0)^2 + (x_1 - 3.0)^2$ と定義されており、この関数の最小値は $x_0=1.0, x_1=3.0$ のときに $0.0$ となります。プログラムは初期値 x とステップサイズ dx から探索を開始し、指定された許容誤差 tolx (変数の変化の許容誤差) および tolf (関数値の変化の許容誤差) に達するか、最大反復回数 nmaxiter に達するまで探索を続けます。

ただし、提供されたソースコードは、Pythonの構文エラー（連続する else 文）を含んでおり、また、C言語の printf 関数がPythonの print に適切に置換されていないため、そのままでは実行できません。さらに、配列のインデックスアクセスに起因する潜在的な IndexError を引き起こす可能性があります。

原理

プログラムは、Nelder-Mead法（シンプレックス法）と呼ばれる数値最適化アルゴリズムを実装しています。これは、勾配情報を使用しない直接探索法の一種であり、以下の基本的なステップを繰り返すことで関数の最小値を探索します。

シンプレックスの初期化: $N$ 変数の最適化問題に対して、$N+1$ 個の頂点を持つ単体（シンプレックス）を初期化します。各頂点には変数 x の値とそれに対応する関数値 f(x) が格納されます。
点の分類: シンプレックス内の頂点を関数値の大小で分類します。
- ih: 最も関数値が高い（悪い）頂点のインデックス
- il: 最も関数値が低い（良い）頂点のインデックス
- fs: 二番目に高い関数値
収束判定: シンプレックスのサイズ（頂点間の距離）や関数値の変化が、それぞれ tolx および tolf の許容誤差を下回った場合、収束したと判断し、探索を終了します。
幾何学的操作: 収束しない場合、シンプレックスの形状を変化させる以下のいずれかの操作を行います。
- 重心の計算: 最も関数値が高い点 x_h を除く残りの $N$ 個の点の重心 $\bar{x}$ を計算します。 $$ \bar{x}_j = \frac{1}{N} \sum_{i \neq h} x_{ij} $$
- 反射 (Reflection): 最も悪い点 x_h を重心 $\bar{x}$ に対して反射し、新しい点 $x_r$ を生成します。コードでは係数 alp1 を用いて次のように計算されます（インデックスはコードの記述に従います）： $$ x_{r,j} = \text{alp1} \cdot \bar{x}_{j+1} - x_{h,j+1} $$ ここで `alp1 = 2.0` です。これは標準的な反射式 $x_r = \bar{x} + \alpha(\bar{x} - x_h)$ で $\alpha=1$ とした場合と同じ効果を持ちます。
- 拡張 (Expansion): 反射点 $x_r$ が現在の最良点 $x_l$ よりも良い場合、シンプレックスをさらに $x_r$ の方向に拡張し、新しい点 $x_e$ を生成します。コードでは係数 gam1 を用いて次のように計算されます： $$ x_{e,j+1} = 2.0 \cdot x_{r,j-1} + \text{gam1} \cdot x_{j+1} $$ ここで `gam1 = -1.0` です。これは標準的な拡張式 $x_e = \bar{x} + \gamma(x_r - \bar{x})$ で $\gamma=2$ とした場合と同じ効果を持ちます。
- 収縮 (Contraction): 反射点 $x_r$ が最良点ほど良くない場合や、最悪点よりも悪い場合、シンプレックスを収縮させ、新しい点 $x_c$ を生成します。コードでは係数 bet1 を用いて次のように計算されます： $$ x_{c,j} = 0.5 \cdot x_{h,j} + \text{bet1} \cdot x_j $$ ここで `bet1 = 0.5` です。これは $x_{c,j} = (x_{h,j} + \bar{x}_j)/2$ という、最悪点と重心の中点を取る操作に相当します。
- 縮小 (Reduction): 収縮操作でも改善が見られない場合、最良点 x_l を中心としてシンプレックス全体を縮小します。 $$ x_{i,j} = (x_{i,j} + x_{l,j}) / 2.0 $$
反復: 収束条件が満たされるか、最大反復回数に達するまでステップ2から4を繰り返します。

必要な非標準ライブラリとインストール方法

このプログラムは以下の非標準ライブラリを使用します。

numpy: 数値計算（配列操作など）
matplotlib: グラフ描画（ただし、このプログラムのバージョンでは実際にグラフ描画は行われていません）

これらのライブラリは pip コマンドを使用してインストールできます。

pip install numpy matplotlib

必要な入力ファイル

このプログラムは外部の入力ファイルを必要としません。目的関数 func(x) はプログラムのソースコード内に直接定義されており、初期値や探索パラメータも main 関数内でハードコードされています。

生成される出力ファイル

このプログラムは、いかなる出力ファイルも生成しません。最適化の進行状況や結果は、標準出力（コンソール）に表示されます。具体的には、反復回数（ITER）、最低関数値（FL）、二番目に高い関数値（FS）、最高関数値（FH）が表示されます。また、収束時や最大反復回数に達した際には、対応するメッセージが出力されます。

コマンドラインでの使用例 (Usage)

このプログラムは引数なしで実行されます。

python optimize-simplex.py

コマンドラインでの具体的な使用例

提供されたソースコードは、複数の問題を含んでいるため、そのままでは実行できません。

SyntaxError: 154行目と155行目に else: が連続して記述されており、Pythonの構文として不正です。これにより、プログラムの実行前に SyntaxError: invalid syntax が発生します。
NameError: printf 関数がPythonの print に適切に置換されていないため、もし構文エラーを修正して実行できたとしても、NameError: name 'printf' is not defined が発生します。
IndexError: コードには、配列のインデックスアクセスに j+1 や j-1 といったオフセットが多用されており、特に n=2 のような次元の低い問題では、IndexError: index is out of bounds for axis が発生する可能性が高いです。例えば、simplex 関数の103行目や144行目などで発生しうる問題です。

これらの問題を修正し、プログラムが意図通りに動作したと仮定した場合の出力例は以下のようになります。

$ python optimize-simplex.py
ITER=  0  FL, FS, FH=2.0,5.0,10.0
ITER= 10  FL, FS, FH=0.0000216,0.0000490,0.0000854
ITER= 20  FL, FS, FH=0.0000000,0.0000000,0.0000000
(SUBR.SMPLX) CONVERGENCE AT ITER=25

実行結果の説明:

上記の例は、プログラムがエラーなく実行された場合の、最適化の進行状況を示しています。

ITER= 0 FL, FS, FH=2.0,5.0,10.0: 0回目の反復（初期状態）でのシンプレックス内の関数値の分布を示します。FL が最小値、FS が二番目に高い値、FH が最大値です。
ITER= 10 FL, FS, FH=0.0000216,0.0000490,0.0000854: 10回目の反復では、関数値が大幅に減少し、最小値 $0$ に近づいていることがわかります。
ITER= 20 FL, FS, FH=0.0000000,0.0000000,0.0000000: 20回目の反復では、関数値がほぼ $0$ に到達しています。
(SUBR.SMPLX) CONVERGENCE AT ITER=25: 25回目の反復で収束条件が満たされたことを示し、最適化が完了したことを報告しています。このとき、目的関数 $f(x_0, x_1) = (x_0 - 1.0)^2 + (x_1 - 3.0)^2$ の最小値 $0$ に近い解が見つかっていることを意味します。最終的な変数 x の値は、[1.0, 3.0] に近いものになります。