optimize.py 技術ドキュメント

プログラムの動作

optimize.py は、勾配降下法、共役勾配法、ニュートン法、シンプレックス法といった様々な数値最適化アルゴリズムを用いて、2変数の目的関数に対する最小値を探索するPythonプログラムです。主に、与えられた初期値から出発して関数の最小点を見つけることを目的としています。

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

  • 多様な最適化アルゴリズムのサポート: 最急降下法 (sd)、共役勾配法 (cg)、ニュートン法 (newton)、準ニュートン法 (broyden, dfp, bfgs)、およびシンプレックス法 (simplex) を選択できます。

  • 複数の目的関数の選択: 楕円体関数 (ellipsoid, ellipsoid2)、円関数 (circle)、またはより複雑な4次関数 (general) から選択できます。

  • ラインサーチ方式の指定: 各反復で最適なステップサイズを見つけるためのラインサーチ手法 (例: armijo, golden) を設定できます。

  • リアルタイムでの最適化過程の可視化: 目的関数の等高線プロットと3Dワイヤーフレーム上に、探索経路をリアルタイムで表示し、最適化の進行状況を視覚的に確認できます。

  • コマンドライン引数による柔軟な設定: 初期値、アルゴリズム、目的関数、グラフの表示設定などをコマンドライン引数で指定できます。

  • 最適化結果の保存: 最終的な最適化結果(最小点と最小値)をファイルに保存します。

このプログラムは、多変数関数の最小化問題に取り組む際の強力なツールとして機能し、異なる最適化手法の挙動を比較検討するのに役立ちます。

原理

このプログラムは、2変数関数 \(f(x_0, x_1)\) の最小値を探索します。最適化アルゴリズムの多くは、関数の勾配 \(\nabla f(\mathbf{x})\) やヘッセ行列 \(\mathbf{H} f(\mathbf{x})\) を利用して探索方向とステップサイズを決定します。

勾配ベクトルは次のように定義されます。 $\( \nabla f(\mathbf{x}) = \begin{pmatrix} \frac{\partial f}{\partial x_0} \\ \frac{\partial f}{\partial x_1} \end{pmatrix} \)$

ヘッセ行列は次のように定義されます。 $\( \mathbf{H} f(\mathbf{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} \)$

プログラムで利用可能な主な最適化アルゴリズムと、それらが利用する情報の概要は以下の通りです。

  1. シンプレックス法 (simplex): 関数の勾配情報を使用せず、点の集合(シンプレックス)を反転、伸縮、収縮させて関数の谷底に向かって移動させます。

  2. 最急降下法 (sd): 現在の点における勾配の負の方向 (\(\mathbf{d} = -\nabla f(\mathbf{x})\)) に関数値が最も急峻に減少するため、この方向に沿って探索を行います。ステップサイズはラインサーチで決定されます。

  3. 共役勾配法 (cg): 最急降下法を改良したもので、探索方向が過去の勾配と「共役」になるように調整され、探索の効率を高めます。勾配情報を使用します。

  4. ニュートン法 (newton): 関数の2次近似(テイラー展開)に基づき、ヘッセ行列の逆行列を用いて探索方向 \(\mathbf{d} = -(\mathbf{H} f(\mathbf{x}))^{-1} \nabla f(\mathbf{x})\) を決定します。勾配とヘッセ行列の両方の情報が必要です。

  5. 準ニュートン法 (broyden, dfp, bfgs): ニュートン法と同様に2次近似を利用しますが、ヘッセ行列やその逆行列を直接計算する代わりに、勾配情報のみからヘッセ行列の近似を構築します。これにより、2次微分の計算コストを削減しつつ、ニュートン法に近い収束性能を目指します。

探索は通常、各ステップで最適なステップサイズ \(\alpha\) を見つけるための「ラインサーチ」と組み合わせて行われます。ラインサーチ手法には armijo (アルミホ条件), golden (黄金分割探索), exact (厳密ラインサーチ) などがあります。

プログラムで定義されている目的関数 (func) は以下の通りです。

  • ellipsoid タイプ: $\(f(x_0, x_1) = x_0^2 + 9.0x_1^2\)$

  • ellipsoid2 タイプ: 定数 \(ae2xx = 1.0\), \(ae2xy = -2.0\), \(ae2yy = 3.0\) を用いて、 $\(f(x_0, x_1) = ae2xx \cdot x_0^2 + ae2xy \cdot x_0x_1 + ae2yy \cdot x_1^2\)$

  • circle タイプ: 定数 \(acxx = 3.0\) を用いて、 $\(f(x_0, x_1) = acxx \cdot (x_0^2 + x_1^2)\)$

  • general タイプ (デフォルト): $\(f(x_0, x_1) = -3.0 - 10.0x_0 - 30.0x_0^2 + 1.5x_0^3 + 3.0x_0^4 + 30.0x_1 - 30.0x_1^2 + 3.0x_1^4 + 3.0x_0x_1^2\)$

これらの関数に対し、それぞれに対応する1次微分 (diff1) と2次微分 (diff2、ヘッセ行列要素) もプログラム内で定義されており、選択された最適化アルゴリズムに応じて利用されます。

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

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

  1. NumPy: 数値計算を効率的に行うためのPythonライブラリです。

    • インストール方法:

      pip install numpy
      
  2. tklib: このプログラムで使用されている tkargs, tksci.tkoptimize, tkgraphic.tkplot モジュールを含むカスタムライブラリです。

    • インストール方法: tklib は標準の pip でインストールできる公開ライブラリではありません。このプログラムの実行には、tklib ディレクトリがPythonのパスから利用可能であるか、またはプログラムファイルと同じディレクトリに配置されている必要があります。通常、これはプログラム作者が用意したカスタムライブラリです。別途入手するか、プログラムに同梱されていることを想定しています。

必要な入力ファイル

本プログラムの実行に必要な入力ファイルはありません。 初期パラメータはプログラム内でデフォルト値が設定されており、コマンドライン引数によって上書きすることが可能です。

生成される出力ファイル

プログラムの実行が完了すると、以下のファイルが生成されます。

  • final_params.prm: 最適化によって見つけられた最終的なパラメータ(\(x_0, x_1\))とその点での関数値 \(f(x_0, x_1)\) が保存されます。ファイル形式はキーと値のペアからなるシンプルなテキスト形式です。

    内容例:

    x: 1.23456789
    y: 9.87654321
    f: -50.1234567
    

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

optimize.py はコマンドライン引数を受け取って動作を制御します。引数は以下の順序で指定します。

python optimize.py [x0_x] [x0_y] [algorism] [lsmode] [functype] [colormap] [nlevels] [initial_simplex] [simplex_scale]

重要: このプログラムの主要な最適化処理は、コマンドライン引数が最低9個(スクリプト名を含めると sys.argv の要素数が10個)指定された場合にのみ実行されます。それ未満の引数では、一部の初期情報が表示されるだけで最適化は行われません。

  • [x0_x] (必須): 初期探索点のX座標 (\(x_0\) の初期値)。浮動小数点数。

  • [x0_y] (必須): 初期探索点のY座標 (\(x_1\) の初期値)。浮動小数点数。

  • [algorism] (必須): 使用する最適化アルゴリズム。文字列。

    • 選択肢: simplex, sd (最急降下法), cg (共役勾配法), newton, broyden, dfp, bfgs

  • [lsmode] (必須): ラインサーチのモード。文字列。

    • 選択肢: none, newton, one, simple, exact, golden, armijo

  • [functype] (必須): 最適化対象の関数タイプ。文字列。

    • 選択肢: ellipsoid, ellipsoid2, circle, general

  • [colormap] (必須): 等高線プロットのカラーマップ。文字列。

    • 選択肢: line, hsv, cool, winter, gray, gist_gray, cividis, Spectral など。

  • [nlevels] (必須): 等高線の数。整数。

  • [initial_simplex] (必須): シンプレックス法で初期シンプレックスを生成する方法 (現在はコード内で使用されていない可能性が高いが、引数として要求される)。文字列。

  • [simplex_scale] (必須): シンプレックス法の初期シンプレックスのサイズを決定するためのスケール値。浮動小数点数。

引数を省略した場合、プログラム内部で定義されたデフォルト値が使用されますが、前述の通り、引数の総数が10未満の場合には最適化処理が実行されない点に注意してください。

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

以下の例では、プログラムが最低9個の引数(スクリプト名を含めると sys.argv の要素数が10個)を必要とすることに基づいてコマンドを記述します。

例1: 一般的な4次関数を最急降下法 (SD) で最適化

初期点を \((1.0, 1.0)\) とし、general 関数を sd (最急降下法) と armijo ラインサーチで最適化します。

python optimize.py 1.0 1.0 sd armijo general cool 51 default 1.0

実行結果の説明: このコマンドを実行すると、初期点 \((1.0, 1.0)\) から開始し、general 関数を sd アルゴリズムで最小化する過程がリアルタイムでグラフ表示されます。探索の進行とともに、最適化経路が等高線グラフ上にプロットされます。最終的に、最適化されたパラメータと関数値がコンソールに出力され、final_params.prm ファイルに保存されます。

コンソール出力例 (最終部分):

Find minimum point by steepest-descend / conjugate gradient / simplex methods

x0= (1.0, 1.0)
algorism= sd
  For simplex:
    initial_simplex= default
    scale= 1.0
lsmode= armijo
functype= general
colormap= cool

Initial parameters:
  x(1):    1.000000000000000e+00  (dx=1.000000e+00)
  y(1):    1.000000000000000e+00  (dx=1.000000e+00)
Initial func:  -4.500000e+01

x0 = (1.0, 1.0): f = -45.0

Optimized parameters:
  x(1):    3.058359230554030e+00  (dx=1.000000e+00)
  y(1):   -0.000000000000000e+00  (dx=1.000000e+00)
Optimized func:  -1.393433e+02

Final resuls are stored to [final_params.prm]
Press enter to terminate:

この例では、初期点 \((1.0, 1.0)\) から最適化を開始し、最終的に \((3.058, 0.0)\) 付近の点で関数値が \(-139.34\) に達したことが示されています。

例2: 楕円体関数をシンプレックス法で最適化

初期点を \((-2.0, -2.0)\) とし、ellipsoid 関数を simplex アルゴリズムで armijo ラインサーチ(シンプレックス法ではラインサーチは直接使用されないが、引数として指定する必要がある)で最適化します。シンプレックスの初期サイズを 0.5 に設定します。

python optimize.py -2.0 -2.0 simplex armijo ellipsoid hsv 31 default 0.5

実行結果の説明: このコマンドは、初期点 \((-2.0, -2.0)\) から ellipsoid 関数に対してシンプレックス法を実行します。hsv カラーマップで31本の等高線が描画され、シンプレックスの動きがグラフ上に表示されます。ellipsoid 関数の最小値は \((0, 0)\) です。

コンソール出力例 (最終部分):

Find minimum point by steepest-descend / conjugate gradient / simplex methods

x0= (-2.0, -2.0)
algorism= simplex
  For simplex:
    initial_simplex= default
    scale= 0.5
lsmode= armijo
functype= ellipsoid
colormap= hsv

Initial parameters:
  x(1):   -2.000000000000000e+00  (dx=5.000000e-01)
  y(1):   -2.000000000000000e+00  (dx=5.000000e-01)
Initial func:   4.000000e+01

x0 = (-2.0, -2.0): f = 40.0

Optimized parameters:
  x(1):    0.000000000000000e+00  (dx=5.000000e-01)
  y(1):    0.000000000000000e+00  (dx=5.000000e-01)
Optimized func:  -3.418434e-15

Final resuls are stored to [final_params.prm]
Press enter to terminate:

この例では、最適化により \((0.0, 0.0)\) 付近の点で関数値が非常に小さい(ほぼゼロ)値に収束したことが示されています。これは ellipsoid 関数の期待される最小値です。

例3: 一般的な4次関数をBFGS法で最適化

初期点を \((0.0, 0.0)\) とし、general 関数を bfgs (準ニュートン法の一種) と armijo ラインサーチで最適化します。

python optimize.py 0.0 0.0 bfgs armijo general cividis 61 default 1.0

実行結果の説明: このコマンドは、初期点 \((0.0, 0.0)\) から general 関数に対して bfgs アルゴリズムを実行します。cividis カラーマップで61本の等高線が描画され、最適化経路が表示されます。

コンソール出力例 (最終部分):

Find minimum point by steepest-descend / conjugate gradient / simplex methods

x0= (0.0, 0.0)
algorism= bfgs
  For simplex:
    initial_simplex= default
    scale= 1.0
lsmode= armijo
functype= general
colormap= cividis

Initial parameters:
  x(1):    0.000000000000000e+00  (dx=1.000000e+00)
  y(1):    0.000000000000000e+00  (dx=1.000000e+00)
Initial func:   -3.000000e+00

x0 = (0.0, 0.0): f = -3.0

Optimized parameters:
  x(1):    3.058359230554030e+00  (dx=1.000000e+00)
  y(1):   -0.000000000000000e+00  (dx=1.000000e+00)
Optimized func:  -1.393433e+02

Final resuls are stored to [final_params.prm]
Press enter to terminate:

この例では、bfgs アルゴリズムも同様に \((3.058, 0.0)\) 付近の点で関数値が \(-139.34\) に収束したことが示されており、例1の sd アルゴリズムと同じ局所最小点に到達しました。