optimize-sd-cg2d-linesearch.py

プログラムの動作

optimize-sd-cg2d-linesearch.py は、2次元の連続関数における最小値を探索するためのPythonプログラムです。このプログラムは、最急降下法（Steepest Descent, SD）と共役勾配法（Conjugate Gradient, CG）という2つの主要な最適化アルゴリズムを実装しており、それぞれ異なる直線探索（Line Search）手法と組み合わせて使用できます。

主な機能は以下の通りです。

2次元関数の最適化: ユーザーが定義した（またはプログラムに組み込まれた）2次元関数の局所的最小点を探索します。
複数の最適化アルゴリズム: 最急降下法と共役勾配法の選択が可能です。
多様な直線探索手法: ステップサイズを決定するための「固定ステップ（one）」、「単純探索（simple）」、「厳密解法（exact）」、「黄金分割法（golden）」、「Armijo条件（armijo）」の5つの手法をサポートします。
関数の切り替え: 楕円体関数（ellipsoid）とより複雑な一般関数（general）のいずれかを最適化対象として選択できます。
リアルタイム可視化: 最適化の進行状況を等高線図と3Dサーフェス図でリアルタイムにプロットし、探索経路を視覚的に追跡できます。

このプログラムは、多変数関数の最小値探索における基本的な最適化手法の動作原理を理解し、異なるアルゴリズムや直線探索手法が収束挙動に与える影響を比較・評価するのに役立ちます。

原理

このプログラムは、2変数関数 $f(\mathbf{x})$ の最小値を探索する非制約最適化問題に取り組んでいます。ここで $\mathbf{x} = (x_0, x_1)^T \in \mathbb{R}^2$ です。最適化は以下の反復プロセスで行われます。

\[\mathbf{x}_{k+1} = \mathbf{x}_k + \alpha_k \mathbf{d}_k\]

ここで、$\mathbf{x}_k$ は $k$ 回目の反復における現在の点、$\mathbf{d}_k$ は探索方向、$\alpha_k$ はステップサイズ（学習率）です。

勾配（Gradient）

関数の勾配 $\nabla f(\mathbf{x})$ は、各変数に関する偏導関数ベクトルとして定義されます。

\[\nabla f(\mathbf{x}) = \left( \frac{\partial f}{\partial x_0}(\mathbf{x}), \frac{\partial f}{\partial x_1}(\mathbf{x}) \right)^T\]

このプログラムでは、各変数 $x_i$ に関する1次偏導関数 diff1(i, x) を用いて勾配を計算します。

最適化アルゴリズム

最急降下法 (Steepest Descent, SD) 最急降下法は、関数の値が最も急峻に減少する方向（負の勾配方向）を探索方向とします。

\[\mathbf{d}_k = -\nabla f(\mathbf{x}_k)\]
共役勾配法 (Conjugate Gradient, CG) 共役勾配法は、各探索方向が互いに共役であるように選択することで、最急降下法よりも効率的な探索を行います。探索方向 $\mathbf{d}_k$ は、現在の負の勾配と前回の探索方向 $\mathbf{d}_{k-1}$ の線形結合として定義されます。

\[\mathbf{d}_k = -\nabla f(\mathbf{x}_k) + \beta_k \mathbf{d}_{k-1}\]

ここで、$\beta_k$ は共役性係数であり、このプログラムではHestenes-Stiefel (HS) の公式を用いて計算されます。

\[\beta_k = \frac{\nabla f(\mathbf{x}_k)^T (\nabla f(\mathbf{x}_k) - \nabla f(\mathbf{x}_{k-1}))}{\mathbf{d}_{k-1}^T (\nabla f(\mathbf{x}_k) - \nabla f(\mathbf{x}_{k-1}))}\]

コードの実装では、gradfk が $\nabla f(\mathbf{x}_k)$、gradfkm が $\nabla f(\mathbf{x}_{k-1})$、dkm が $\mathbf{d}_{k-1}$ に対応します。 yk = gradfk - gradfkm inn1 = np.inner(gradfk, yk) inn2 = np.inner(dkm, yk) そして $\beta_k = \text{inn1} / \text{inn2}$ として計算されています。 icg カウンタがパラメータの次元数（ここでは2）に達すると、共役勾配方向はリセットされ、最急降下方向から再開します。

直線探索 (Line Search) 手法

直線探索は、与えられた探索方向 $\mathbf{d}_k$ に対して最適なステップサイズ $\alpha_k$ を決定するプロセスです。目的は、1次元関数 $\phi(\alpha) = f(\mathbf{x}_k + \alpha \mathbf{d}_k)$ を最小化する $\alpha > 0$ を見つけることです。

linesearch_one (固定ステップ): ステップサイズ $\alpha_k$ は、初期設定された alpha 値に固定されます。 $$\alpha_k = \text{alpha}$$
linesearch_simple (単純探索): 初期ステップサイズ alpha から開始し、それを倍々 (i * alpha) にしながら関数値が減少する限り探索を続けます。関数値が増加に転じたら、一つ前のステップサイズを採用します。
linesearch_exact (厳密解法): ステップサイズ $\alpha_k$ を決定するために、1次元関数 $\phi(\alpha)$ を2次関数で近似します。最適な $\alpha$ は、近似関数の導関数がゼロになる点で見つけることができます。近似導関数 $\phi'(\alpha) \approx \phi'(0) + \phi''(0)\alpha$ から $\alpha = -\phi'(0)/\phi''(0)$ を導出します。プログラムでは、中心差分近似を用いて $\phi'(0)$ と $\phi''(0)$ を計算します。 $f_1 = \phi'(0) \approx \frac{f(\mathbf{x}_k + h \mathbf{sdir}/||\mathbf{sdir}||) - f(\mathbf{x}_k - h \mathbf{sdir}/||\mathbf{sdir}||)}{2h}$ $f_2 = \phi''(0) \approx \frac{f(\mathbf{x}_k + h \mathbf{sdir}/||\mathbf{sdir}||) - 2f(\mathbf{x}_k) + f(\mathbf{x}_k - h \mathbf{sdir}/||\mathbf{sdir}||)}{h^2}$ として、ステップサイズ $\alpha_2 = f_1 / (|f_2| + \text{dump})$ を更新に利用します。
linesearch_golden (黄金分割法): 単峰性の関数を前提とし、探索区間を黄金比に基づいて効率的に縮小していく手法です。初期探索範囲 xrange を持ち、$\eta = (\sqrt{5}-1)/(\sqrt{5}+1) \approx 0.381966$ を用いて内側の点を選択し、関数値に基づいて区間を更新します。
linesearch_armijo (Armijo条件): Armijo条件は、関数値の十分な減少を保証するための基準です。初期ステップサイズ alpha から開始し、条件 $f(\mathbf{x}_k + \alpha \mathbf{d}_k) \le f(\mathbf{x}_k) - \xi \alpha ||\mathbf{d}_k||^2$ が満たされるまで $\alpha$ を減少させます。ここで、$\xi$ はArmijoパラメータ（コードでは xi）、$||\mathbf{d}_k||^2$ は探索方向 $\mathbf{d}_k$ のL2ノルムの2乗です。

対象関数

このプログラムには、以下の2つの2次元関数が組み込まれています。

ellipsoid 関数 (efunc): $$f(x_0, x_1) = x_0^2 + 9x_1^2$$ 勾配 `ediff1`: $$\frac{\partial f}{\partial x_0} = 2x_0$$\frac{\partial f}{\partial x_1} = 18x_1$$ この関数は $(0,0)$ で最小値 $0$ を持ちます。
general 関数 (gfunc): $$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$$ 勾配 `gdiff1`: $$\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$$ この関数はより複雑な形状を持ち、複数の局所的最小値を持つ可能性があります。

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

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

NumPy: 数値計算（特に配列操作や線形代数）に使用されます。
Matplotlib: データのプロットと可視化に使用されます。3Dプロットのためには mpl_toolkits.mplot3d も必要ですが、これは通常 Matplotlib に含まれています。

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

pip install numpy matplotlib

必要な入力ファイル

optimize-sd-cg2d-linesearch.py は、外部の入力ファイルを必要としません。すべての設定（初期点、アルゴリズム、直線探索モード、対象関数など）は、プログラム内の定数として設定されているか、コマンドライン引数として直接渡されます。

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

このプログラムは、ディスク上にいかなるファイルも生成しません。実行結果は以下の形式で出力されます。

標準出力:
- 最適化プロセスの概要（選択されたアルゴリズム、直線探索モード、対象関数）。
- 各反復ステップ（イテレーション）における以下の情報:
  - 反復回数（iter）
  - 直線探索の内部イテレーション回数（insiter）
  - 現在の点の座標 $(x_0, x_1)$
  - 前回のステップからの移動量 $dx = \max(|dx_0|, |dx_1|)$
  - 現在の関数値 $f(x_0, x_1)$
- 収束メッセージ、または非収束メッセージ。
GUIウィンドウ:
- 対象関数の3Dワイヤーフレーム図。
- 対象関数の等高線図。この等高線図には、最適化プロセス中に探索された点の軌跡が青い線でリアルタイムに描画されます。

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

プログラムは以下の形式でコマンドラインから実行できます。

python optimize-sd-cg2d-linesearch.py [x0_val] [x1_val] [algorithm] [linesearch_mode] [function_type]

引数は以下の通りです。

[x0_val] (オプション): 初期点の $x_0$ 座標。浮動小数点数。
[x1_val] (オプション): 初期点の $x_1$ 座標。浮動小数点数。
[algorithm] (オプション): 使用する最適化アルゴリズム。文字列 'sd' (最急降下法) または 'cg' (共役勾配法)。
[linesearch_mode] (オプション): 使用する直線探索モード。文字列 'one', 'simple', 'exact', 'golden', 'armijo' のいずれか。
[function_type] (オプション): 最適化対象の関数。文字列 'ellipsoid' または 'general'。

デフォルト値: 引数が指定されない場合、以下のデフォルト値が使用されます。

x0_val: 0.0
x1_val: 0.0
algorithm: 'cg'
linesearch_mode: 'simple'
function_type: 'general'

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

例1: 最急降下法、Armijo直線探索、ellipsoid関数

初期点を $(3.0, 3.0)$ とし、最急降下法 (sd) とArmijo直線探索 (armijo) を用いて ellipsoid 関数を最適化します。

python optimize-sd-cg2d-linesearch.py 3.0 3.0 sd armijo ellipsoid

実行結果の説明:

プログラムは、初期点 $(3.0, 3.0)$ から開始し、最急降下法とArmijo条件に基づくステップサイズで、関数 $f(x_0, x_1) = x_0^2 + 9x_1^2$ の最小値 $(0,0)$ を探索します。

標準出力には、以下のような内容が表示されます（一部抜粋）。

Find minimum point by steepest-descend / conjugate gradient methods

x0 = (3.0, 3.0): f = 90.0
   0(  4): x = (  2.625000,   2.625000) dx =   0.375000: f =   65.812500
   1(  3): x = (  2.321875,   2.321875) dx =   0.303125: f =   48.514063
   ...
 198(  3): x = ( 0.000003,  0.000000) dx =  0.000000: f =  0.000000
 199(  3): x = ( 0.000003,  0.000000) dx =  0.000000: f =  0.000000
Converged at x =  [3.03362153e-06 3.37069059e-07]  dx = 3.033621531777085e-06   f = 9.202359218228394e-12

GUIウィンドウには、$(3.0, 3.0)$ から $(0,0)$ へと向かう探索経路が等高線図上に青い線でリアルタイムに描画され、3Dワイヤーフレーム図とともに表示されます。最終的には、$(0,0)$ に非常に近い点で関数値が最小となり、収束したことが報告されます。

例2: 共役勾配法、黄金分割法、general関数

初期点を $(1.0, -1.0)$ とし、共役勾配法 (cg) と黄金分割法 (golden) を用いて general 関数を最適化します。

python optimize-sd-cg2d-linesearch.py 1.0 -1.0 cg golden general

実行結果の説明:

プログラムは、初期点 $(1.0, -1.0)$ から開始し、共役勾配法と黄金分割法に基づくステップサイズで、より複雑な general 関数の最小値を探索します。

標準出力には、以下のような内容が表示されます（一部抜粋）。

Find minimum point by steepest-descend / conjugate gradient methods

x0 = (1.0, -1.0): f = 7.5
   0(  0): x = (  0.428571,  -0.342857) dx =   0.571429: f =  -30.825816
   1(  0): x = (  0.490060,  -0.342857) dx =   0.061489: f =  -30.835497
   ...
  45(  0): x = (  0.500000,  -0.499999) dx =  0.000000: f =  -31.187500
  46(  0): x = (  0.500000,  -0.499999) dx =  0.000000: f =  -31.187500
Converged at x =  [ 0.5 -0.5]  dx = 5.679198642959604e-07   f = -31.1875

GUIウィンドウには、$(1.0, -1.0)$ から局所的な最小値へと向かう探索経路が描画され、3Dワイヤーフレーム図とともに表示されます。この general 関数には複数の局所解が存在する可能性がありますが、この実行例では約 $(0.5, -0.5)$ 付近で関数値 $-31.1875$ の最小値に収束したことが報告されます。