optimize-newton-raphson2d.py テクニカルドキュメント

プログラムの動作

optimize-newton-raphson2d.py は、2変数の非線形関数 $f(x_0, x_1)$ の極値（最小値または最大値）をニュートン・ラプソン法を用いて探索するPythonプログラムです。

このプログラムの主な機能は以下の通りです。

極値探索: ソースコードにハードコードされた特定の2変数関数について、その1階および2階偏導関数（ヘッセ行列の要素）を用いてニュートン・ラプソン法による数値最適化を実行します。
リアルタイム可視化: 探索の進行状況をMatplotlibライブラリを用いてリアルタイムで可視化します。具体的には、関数の2D等高線図と3Dワイヤーフレーム図を表示し、探索経路（$x_0, x_1$ の遷移）をプロットします。
初期値の指定: 探索の開始点となる$x_0, x_1$ の初期値をコマンドライン引数で指定できます。
収束判定と終了条件: 更新ステップの最大絶対値が事前に設定された閾値 eps を下回るか、最大反復回数 nmaxiter に達した場合に探索を終了します。

本プログラムは、多変数関数の数値最適化、特に勾配情報（1階および2階導関数）が利用可能な場合の極値探索という課題を解決します。

目的関数 $f(x_0, x_1)$ は以下の式で定義されています。

\[f(x_0, x_1) = -3 - 10x_0 - 30x_0^2 + 1.5x_0^3 + 3x_0^4 + 30x_1 - 30x_1^2 + 3x_1^4 + 3x_0x_1^2\]

1階偏導関数 $\nabla f(x) = \begin{pmatrix} \frac{\partial f}{\partial x_0} \\ \frac{\partial f}{\partial x_1} \end{pmatrix}$ の各要素は、diff1 関数によって計算され、以下の式に対応します。

\[\frac{\partial f}{\partial x_0} = -10 - 60x_0 + 4.5x_0^2 + 12x_0^3 + 3x_1^2\]

\[\frac{\partial f}{\partial x_1} = 30 - 60x_1 + 12x_1^3 + 6x_0x_1\]

2階偏導関数（ヘッセ行列 $H(x)$ の要素）は、diff2 関数によって計算され、以下の式に対応します。

\[\begin{split}H(x) = \begin{pmatrix} \frac{\partial^2 f}{\partial x_0^2} & \frac{\partial^2 f}{\partial x_0 \partial x_1} \\ \frac{\partial^2 f}{\partial x_1 \partial x_0} & \frac{\partial^2 f}{\partial x_1^2} \end{pmatrix}\end{split}\]

ただし、コード内の diff2 関数は次のように定義されています。

\[\frac{\partial^2 f}{\partial x_0^2} = -60 + 9x_0 + 36x_0^2 + 3x_1^2\]

\[\frac{\partial^2 f}{\partial x_0 \partial x_1} = 6x_0x_1\]

\[\frac{\partial^2 f}{\partial x_1 \partial x_0} = 6x_1\]

\[\frac{\partial^2 f}{\partial x_1^2} = -60 + 36x_1^2 + 6x_0\]

原理

本プログラムは、非線形最適化手法の一つであるニュートン・ラプソン法を2変数の関数に適用しています。

ニュートン・ラプソン法は、テイラー展開を用いて関数を近似し、その近似関数の極値を求めることで、元の関数の極値に近づいていく反復法です。現在の探索点 $x_k$ から次の探索点 $x_{k+1}$ を更新する式は以下の通りです。

\[x_{k+1} = x_k - H(x_k)^{-1} \nabla f(x_k)\]

ここで:

$x_k = \begin{pmatrix} x_{k,0} \\ x_{k,1} \end{pmatrix}$ は現在の探索点における変数のベクトルです。
$\nabla f(x_k)$ は現在の探索点における勾配ベクトル（1階偏導関数ベクトル）です。 $$\nabla f(x_k) = \begin{pmatrix} \frac{\partial f}{\partial x_0}(x_k) \\ \frac{\partial f}{\partial x_1}(x_k) \end{pmatrix}$$ プログラムでは Si 変数がこれに対応します。
$H(x_k)$ は現在の探索点におけるヘッセ行列（2階偏導関数行列）です。ヘッセ行列は対称行列であり、各要素は2階偏導関数で構成されます。 $$H(x_k) = \begin{pmatrix} \frac{\partial^2 f}{\partial x_0^2}(x_k) & \frac{\partial^2 f}{\partial x_0 \partial x_1}(x_k) \\ \frac{\partial^2 f}{\partial x_1 \partial x_0}(x_k) & \frac{\partial^2 f}{\partial x_1^2}(x_k) \end{pmatrix}$$ プログラムでは Sij 変数がこれに対応します。
$H(x_k)^{-1}$ はヘッセ行列の逆行列です。
$H(x_k)^{-1} \nabla f(x_k)$ は探索方向とステップ幅を示すベクトルであり、プログラムでは -dx に対応します。

プログラムは、この更新式を反復的に適用し、$\|dx\|$ が非常に小さくなった（収束した）場合に探索を終了します。

また、本プログラムでは探索の安定性を向上させるためにダンピング（減衰）ファクターを導入しています。ヘッセ行列 $H(x_k)$ が正定値でない場合（特に最大値探索や鞍点近傍で不安定になる可能性がある場合）に、単位行列のdump倍をヘッセ行列に加算または減算することで、探索の安定化を図っています。これは Sij += dump * tr * np.eye(n) や Sij -= dump * np.eye(n) の部分に相当します。tr はヘッセ行列のトレース（対角成分の和）です。

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

本プログラムの実行には、以下の非標準Pythonライブラリが必要です。

numpy: 数値計算（配列操作、線形代数演算など）
matplotlib: グラフ描画（2D/3Dプロット、等高線図など）

これらのライブラリは、Pythonのパッケージマネージャーである pip を使用してインストールできます。コマンドプロンプトやターミナルで以下のコマンドを実行してください。

pip install numpy matplotlib

必要な入力ファイル

このプログラムは、特定の外部入力ファイルを必要としません。最適化対象の関数 $f(x_0, x_1)$ とその導関数は、すべてプログラムのソースコード内に直接定義されています。

探索の開始点（初期値 $x_0, x_1$）は、コマンドライン引数として与えることができます。引数が指定されない場合は、デフォルト値として $x_0 = 0.0, x_1 = 0.0$ が使用されます。

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

本プログラムは、実行中にいかなるファイルも生成したり保存したりしません。

計算の進行状況と最終的な結果は、標準出力（コンソール）にテキストとして表示されます。また、Matplotlibウィンドウに表示されるグラフは、探索のリアルタイムな可視化を提供しますが、このグラフ自体がファイルとして保存されることはありません。

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

optimize-newton-raphson2d.py は、以下の基本的な形式でコマンドラインから実行できます。初期探索点 $x_0, x_1$ はオプションで指定可能です。

python optimize-newton-raphson2d.py [initial_x0] [initial_x1]

initial_x0 (オプション): 探索を開始する $x_0$ 座標の初期値です。浮動小数点数で指定します。省略された場合、デフォルト値 0.0 が使用されます。
initial_x1 (オプション): 探索を開始する $x_1$ 座標の初期値です。浮動小数点数で指定します。省略された場合、デフォルト値 0.0 が使用されます。

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

例1: デフォルトの初期値で実行する

初期値を指定せずにプログラムを実行する場合、x0 = 0.0 と x1 = 0.0 が初期探索点として使用されます。

python optimize-newton-raphson2d.py

実行結果の例: プログラムが開始され、以下のような情報が標準出力に表示され、同時にグラフウィンドウが開きます。

Find minimum / maximum point by Newton-Raphson method

x0 = (0.0, 0.0): f = -3.0
   dx= [0.16666667 -0.5       ]
   0: x = (0.166667,    -0.5) dx =      0.5: f = -18.7847
   dx= [0.15835942 -0.2100868 ]
   1: x = (0.325026, -0.710087) dx =  0.210087: f = -25.7516
   dx= [0.09675373 -0.16911634]
   2: x = (0.421779, -0.879203) dx =  0.169116: f = -28.5375
   dx= [0.04618753 -0.12932145]
   3: x = (0.467967,    -1.00852) dx =  0.129321: f = -29.6053
   ...
Converged at x =  [0.5, -2.23606797749979]  dx = 3.65478e-06   f = -30.0
Press enter to terminate: 

グラフウィンドウでは、関数の等高線と3Dワイヤーフレームが表示され、原点 $(0.0, 0.0)$ から極値点 $(0.5, -2.236...)$ へと移動する探索経路が青い線でリアルタイムに描画されます。最終的に、収束した極値点が表示されます。

例2: 特定の初期値で実行する

$x_0 = 1.0$, $x_1 = 2.0$ を初期探索点としてプログラムを実行します。

python optimize-newton-raphson2d.py 1.0 2.0

実行結果の例: プログラムは指定された初期値から探索を開始し、以下のような出力がコンソールに表示されます。

Find minimum / maximum point by Newton-Raphson method

x0 = (1.0, 2.0): f = -27.0
   dx= [-0.23192451  0.0357774 ]
   0: x = (0.768075,   2.035777) dx =  0.231925: f = -28.7905
   dx= [-0.17846549  0.07166316]
   1: x = (0.58961 ,   2.107441) dx =  0.178465: f = -29.6209
   dx= [-0.07689945  0.10667856]
   2: x = (0.512711,   2.21412 ) dx =  0.106679: f = -29.9678
   dx= [-0.00760921  0.02194883]
   3: x = (0.505102,   2.236069) dx =  0.0219488: f = -30.0
   ...
Converged at x =  [0.5, 2.23606797749979]  dx = 3.65478e-06   f = -30.0
Press enter to terminate: 

この場合、グラフ上では $(1.0, 2.0)$ から探索が始まり、極値点 $(0.5, 2.236...)$ へと収束する経路が描画されます。表示される極値は初期値によって異なる局所極値に収束する可能性があります。