gp_simulation_animation.py 技術ドキュメント

プログラムの動作

gp_simulation_animation.py は、ガウス過程回帰 (Gaussian Process Regression, GPR) を用いたベイズ最適化の様々な探索戦略を視覚的にシミュレーションするPythonプログラムです。指定された回数だけ観測データを追加し、そのたびにガウス過程モデルの予測を更新し、その過程をアニメーションとして表示します。最終的に、このアニメーションをGIFファイルとして保存することも可能です。

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

真の関数: 正弦波と余弦波を組み合わせた模擬的な真の関数を定義し、その関数からのノイズを含んだ観測データを取得します。
ガウス過程回帰: 観測データに基づいて、RBFカーネルを用いたガウス過程モデルを構築し、未知の点における関数の予測平均と予測分散を計算します。
探索戦略のシミュレーション: 次に観測する点 (候補点) を選択するための5つの異なる戦略（獲得関数）を実装しています。
- random: 探索範囲内からランダムに点を選択します。
- std: 予測の標準偏差（不確実性）が最大となる点を選択し、探索 (exploration) を重視します。
- max: 予測平均が最大となる点を選択し、活かし (exploitation) を重視します。
- ucb (Upper Confidence Bound): 予測平均と予測標準偏差の線形結合が最大となる点を選択し、探索と活かしのバランスを取ります。
- ei (Expected Improvement): 現在の最良値からの改善が期待される量を最大化する点を選択します。
アニメーション表示と保存: 各ステップで新しいデータが追加され、GPモデルが更新される様子をリアルタイムでプロットし、最終的にこの一連の動作をアニメーションとして表示し、GIFファイルとして保存します。

このプログラムは、ベイズ最適化の基本的な動作原理と、異なる獲得関数がモデルの学習と探索にどのように影響するかを直感的に理解するための教育ツールとして機能します。

原理

このプログラムは、ガウス過程回帰 (GPR) とベイズ最適化 (Bayesian Optimization) の原理に基づいています。

ガウス過程回帰 (GPR)

ガウス過程は、関数の値を多変量ガウス分布でモデル化するノンパラメトリックな確率モデルです。GPRでは、与えられた観測データ $D = \{(x_i, y_i)\}_{i=1}^N$ をもとに、未知の入力 $x_*$ における関数の値 $f(x_*)$ の事後分布を予測します。この事後分布は平均 $\mu(x_*)$ と分散 $\sigma^2(x_*)$ を持ち、それぞれ関数の予測値と不確実性を表します。

プログラムでは、RBF (Radial Basis Function) カーネル（またはSquared Exponentialカーネル）を使用しています。RBFカーネルは、2点間の距離が近いほど類似度が高いと仮定する一般的なカーネル関数です。

\[ k(x_i, x_j) = \sigma_f^2 \exp\left(-\frac{||x_i - x_j||^2}{2l^2}\right) \]

ここで、$\sigma_f$ は信号の標準偏差（sigma_f）、$l$ は特徴量スケール（length_scale）を表します。

観測ノイズ $\sigma_n^2$ が存在する場合、観測データ $y$ が与えられたときの新しい入力 $x_*$ における関数の予測平均 $\mu(x_*)$ と予測分散 $\sigma^2(x_*)$ は、以下の数式で計算されます。

\[ \mu(x_*) = K_*^T (K + \sigma_n^2 I)^{-1} y \]

\[ \Sigma(x_*) = K_{**} - K_*^T (K + \sigma_n^2 I)^{-1} K_* \]

ここで、

$K$ は訓練データ $x_{train}$ 間のカーネル行列です。$K_{ij} = k(x_i, x_j)$。
$K_*$ は訓練データ $x_{train}$ とテストデータ $x_*$ 間のカーネル行列です。$(K_*)_{ij} = k(x_i, x_j^*)$。
$K_{**}$ はテストデータ $x_*$ 間のカーネル行列です。$(K_{**})_{ij} = k(x_i^*, x_j^*)$。
$I$ は単位行列です。
$y$ は観測データ $y_{train}$ です。

プログラムの gp_posterior 関数では、np.linalg.cholesky を用いて、$(K + \sigma_n^2 I)^{-1}$ の計算を数値的に安定かつ効率的に行っています。

ベイズ最適化と獲得関数

ベイズ最適化は、GPRを用いて評価関数をモデル化し、次に評価すべき点を選択するための「獲得関数 (Acquisition Function)」を最大化することで、評価回数が限られている高価な関数の最適化を目指します。本プログラムでは以下の獲得関数が実装されています。

random (ランダム探索): 探索範囲内の点を完全にランダムに選択します。これは比較対象となるベースラインです。
std (標準偏差最大): 予測標準偏差 $\sigma(x)$ が最大となる点を選択します。これにより、モデルが不確実性の高い領域を優先的に探索し、探索 (exploration) を促進します。
max (予測平均最大): 予測平均 $\mu(x)$ が最大となる点を選択します。これは、現在最も良いと予測される点を活かす (exploitation) 戦略です。
ucb (Upper Confidence Bound): 以下の式を最大化する点を選択します。 $$ x_{next} = \arg\max_x \left( \mu(x) + \kappa \sigma(x) \right) $$ ここで、$\mu(x)$ は予測平均、$\sigma(x)$ は予測標準偏差、$\kappa$ は探索の強さを調整するパラメータです。$\kappa$ を大きくすると探索が促進され、小さくすると活かしが促進されます。プログラムでは $\kappa=2.0$ に設定されています。
ei (Expected Improvement): 現在の観測データで得られた最良値 $y_{best}$ から、どれだけ改善が期待できるかを最大化する点を選択します。 $$ x_{next} = \arg\max_x \left( (\mu(x) - y_{best}) \Phi(Z) + \sigma(x) \phi(Z) \right) $$ ここで、
- $y_{best} = \max(y_{train})$ はこれまでの観測値の最大値です。
- $\mu(x)$ と $\sigma(x)$ はそれぞれ予測平均と予測標準偏差です。
- $Z = \frac{\mu(x) - y_{best}}{\sigma(x) + \epsilon}$ で、$\epsilon$ は数値安定性のための小さな値 ($10^{-9}$ を使用) です。
- $\Phi(Z)$ は標準正規分布の累積分布関数 (CDF) です。
- $\phi(Z)$ は標準正規分布の確率密度関数 (PDF) です。 ei は探索と活かしのバランスを効率的に取るとされています。

avoid_duplicate 関数は、数値計算上の問題や学習効率の低下を防ぐため、既に観測された点に非常に近い点を次の候補点として選ばないように、わずかなオフセットを加える処理です。

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

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

numpy: 数値計算を効率的に行うためのライブラリです。
scipy: 科学技術計算ライブラリで、特にここではscipy.stats.norm（正規分布関連の関数）を使用します。
matplotlib: データの可視化やアニメーション作成に使用します。matplotlib.animation.FuncAnimationとmatplotlib.animation.PillowWriterが必要です。

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

pip install numpy scipy matplotlib

必要な入力ファイル

gp_simulation_animation.py は、外部からの入力ファイルを必要としません。プログラム内部で真の関数 f_true が定義され、シミュレーションに必要な観測データはプログラム実行中に生成されます。

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

プログラムは、ガウス過程回帰のシミュレーション過程をアニメーションとして保存することができます。

ファイル名: コマンドライン引数 outfile で指定されます。指定がない場合、または空文字列が指定された場合は、ファイルは保存されず、matplotlibのウィンドウでアニメーションが表示されるのみとなります。
ファイル形式: GIF形式 (.gif) のアニメーションファイルが生成されます。
内容: 各ステップで新しい観測データが追加され、ガウス過程モデルの予測平均（青線）と95%信頼区間（青い帯）が更新されていく様子が描画されます。真の関数（破線）とこれまでの観測データ（黒点）も表示されます。アニメーションのタイトルには、使用された探索モードと現在のデータ数が表示されます。

出力ファイルの例: ucb_animation.gif (ucbモードで実行した場合)

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

gp_simulation_animation.py は、以下の形式でコマンドラインから実行できます。

python gp_simulation_animation.py [mode] [nmaxiter] [outfile]

[mode]:
- 探索戦略を指定する文字列です。以下のいずれかを指定できます。
  - "random": ランダムに次の点をサンプリングします。
  - "std": 予測標準偏差が最大となる点を選択します。
  - "max": 予測平均が最大となる点を選択します。
  - "ucb": Upper Confidence Bound 戦略で点を選択します（デフォルト）。
  - "ei": Expected Improvement 戦略で点を選択します。
- 省略した場合、デフォルト値は "ucb" です。
[nmaxiter]:
- シミュレーションの繰り返し回数を指定する整数です。
- この回数だけ新しい観測点が追加され、アニメーションのフレーム数が決まります。
- 省略した場合、デフォルト値は 20 です。
[outfile]:
- 生成されるアニメーションGIFファイルのファイル名を指定する文字列です。
- 例えば "animation.gif" のように指定します。
- 省略した場合、または空文字列が指定された場合、アニメーションは保存されず、matplotlibのウィンドウで表示されるだけになります。

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

1. ランダム探索で30ステップ実行し、ファイルを保存しない場合

以下のコマンドを実行すると、random モードで30ステップのガウス過程シミュレーションが実行され、そのアニメーションが画面に表示されます。ファイルは保存されません。

python gp_simulation_animation.py random 30

実行結果の説明: matplotlibのウィンドウが開き、30ステップにわたってランダムに選択された観測点が追加され、ガウス過程の予測が更新される様子がアニメーションとして表示されます。ウィンドウのタイトルには「mode=random | データ数 n = ...」と表示されます。

2. UCB戦略で40ステップ実行し、`ucb_animation.gif` として保存する場合

以下のコマンドを実行すると、ucb モードで40ステップのガウス過程シミュレーションが実行され、アニメーションが画面に表示されるとともに、ucb_animation.gif というファイルが現在のディレクトリに保存されます。

python gp_simulation_animation.py ucb 40 ucb_animation.gif

実行結果の説明: matplotlibのウィンドウが開き、40ステップにわたってUCB戦略に基づいて選択された観測点が追加され、ガウス過程の予測が更新される様子がアニメーションとして表示されます。シミュレーション終了後、コンソールに「Saved animation to ucb_animation.gif」というメッセージが表示され、ucb_animation.gif ファイルが生成されます。

3. EI戦略で25ステップ実行し、ファイルを保存しない場合

以下のコマンドを実行すると、ei モードで25ステップのガウス過程シミュレーションが実行され、そのアニメーションが画面に表示されます。ファイルは保存されません。

python gp_simulation_animation.py ei 25

実行結果の説明: matplotlibのウィンドウが開き、25ステップにわたってEI戦略に基づいて選択された観測点が追加され、ガウス過程の予測が更新される様子がアニメーションとして表示されます。ウィンドウのタイトルには「mode=ei | データ数 n = ...」と表示されます。