技術ドキュメント: bayes_gp_plain_shimizu.py

このドキュメントは、ガウス過程 (Gaussian Process: GP) を用いたベイズ最適化を PHYSBO ライブラリを通じて実行する Python プログラム bayes_gp_plain_shimizu.py の技術的な側面を解説します。

プログラムの動作

bayes_gp_plain_shimizu.py は、入力データに基づいてガウス過程モデルを構築し、ベイズ最適化を通じて目的関数の最適値を探索するプログラムです。特に、高価な実験やシミュレーションの回数を減らしつつ、効率的に最適解を見つけ出すことを目的としています。

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

  1. データ読み込みと前処理: ExcelまたはCSV形式の入力ファイルから記述子（説明変数）と目的関数（応答変数）を自動で識別し、読み込みます。ヘッダーの制御コードによって、最大化、最小化、特定値への収束などの最適化モードを指定できます。また、必要に応じて記述子の標準化を行います。

  2. ガウス過程モデルの構築: 既存の観測データ（訓練データ）を用いてガウス過程モデルを構築し、未観測の領域における目的関数の予測平均値と不確実性（標準偏差）を推定します。

  3. ベイズ最適化の実行: PHYSBO ライブラリを利用し、構築したガウス過程モデルと指定された獲得関数（Expected Improvement: EI, Probability of Improvement: PI, Thompson Sampling: TS）に基づいて、次に評価すべき最適な候補点（アクション）を提案します。

  4. 逐次最適化: ユーザーが指定した繰り返し回数に応じて、ベイズ最適化によって提案された最良の候補点を新たな訓練データとして自動的に追加し、モデルを再学習・更新する機能があります。これにより、最適化プロセスを継続的に改善し、より良い解を探索することが可能になります。

  5. 結果の出力: 予測された目的関数の平均値、標準偏差、獲得関数値などの詳細な結果を Excel 形式で出力します。また、PHYSBO の探索履歴も .npz 形式で保存されます。

  6. 結果の可視化:

    • 等高線図: 記述子が複数ある場合、予測された目的関数の平均値と標準偏差の等高線図を2次元で表示します。

    • パレート図: 複数の目的関数がある場合、目的関数間のトレードオフを示すパレート図を表示します。

    • 予測値と獲得関数のプロット: 各データ点に対する目的関数の予測値（平均値と信頼区間）と獲得関数値をプロットし、最適化の進行状況を示します。

    • インタラクティブなデータ探索: グラフ上の点をクリックすることで、対応するデータ点（記述子、予測値、スコアなど）の詳細情報をコンソールに表示する機能があります。

原理

本プログラムは、ベイズ最適化の枠組みの中でガウス過程を代理モデルとして使用しています。

1. ベイズ最適化

ベイズ最適化は、目的関数 \(f(x)\) の最適値を探索するための効率的な手法です。目的関数 \(f(x)\) が直接評価するにはコストが高い（時間、資源など）場合に特に有効です。この手法は、以下の2つの主要なコンポーネントから構成されます。

  • 代理モデル (Surrogate Model): 既知の観測データに基づいて目的関数の挙動を近似するモデル。本プログラムではガウス過程が用いられます。

  • 獲得関数 (Acquisition Function): 代理モデルの予測に基づいて、次に評価すべき入力点 \(x\) を決定するための基準。

ベイズ最適化は、以下のステップを繰り返して最適解に収束します。

  1. 初期の観測データを用いて代理モデルを構築する。

  2. 獲得関数を最大化する点 \(x_{next}\) を探索する。

  3. \(x_{next}\) で目的関数 \(f(x_{next})\) を評価し、新たな観測データを取得する。

  4. 新たな観測データを既存のデータセットに追加し、代理モデルを更新する。

  5. ステップ2に戻り、最適化を続ける。

2. ガウス過程 (Gaussian Process, GP)

ガウス過程は、任意の有限個の入力点における関数値の組が、同時ガウス分布に従う確率分布として定義されます。これにより、関数全体にわたる不確実性をモデル化することができます。

  • 平均関数 (Mean Function): 事前分布における関数値の期待値 \(m(x) = E[f(x)]\)。通常、0に設定されることが多いです。

  • 共分散関数 (Covariance Function, カーネル関数): 2つの入力点 \(x_i\) と \(x_j\) 間の類似度を測り、それに応じて出力値 \(f(x_i)\) と \(f(x_j)\) の相関関係を定義する関数 \(k(x_i, x_j)\)。例えば、Squared Exponential (RBF) カーネルは以下のように表されます。 $\(k(x_i, x_j) = \sigma_f^2 \exp \left( -\frac{||x_i - x_j||^2}{2l^2} \right)\)\( ここで、\)\sigma_f^2\( は信号の分散、\)l$ は特性長スケールです。

観測データ \(\{ (x_i, y_i) \}_{i=1}^N\) が与えられたとき、ガウス過程は未観測の点 \(x^*\) における関数値 \(f(x^*)\) の事後分布を予測します。この事後分布もガウス分布であり、その平均と分散は解析的に求めることができます。予測される平均値は、目的関数の推定値として、分散は推定値の不確実性として解釈されます。

3. 獲得関数 (Acquisition Function)

獲得関数は、探索 (Exploration: 未知の領域を探索する) と活用 (Exploitation: 既に良い値が見つかっている領域を詳しく探索する) のトレードオフをバランスさせ、次に最も情報量の多い点を提案します。本プログラムでサポートされる主要な獲得関数は以下の通りです。

  • Expected Improvement (EI): 現在の最良値 \(f_{best}\) をどれだけ改善できるかの期待値です。 $\(EI(x) = E[\max(f(x) - f_{best}, 0)]\)\( \)f(x)$ の事後分布がガウス分布である場合、この期待値は解析的に計算できます。

  • Probability of Improvement (PI): 現在の最良値 \(f_{best}\) を改善できる確率です。 $\(PI(x) = P(f(x) > f_{best})\)$

  • Thompson Sampling (TS): 代理モデルの事後分布から関数をサンプリングし、そのサンプリングされた関数の最大値を与える点を次の評価点とする手法です。

4. 標準化

記述子の標準化は、各記述子のスケールを均一にすることで、ガウス過程モデルの学習を安定させ、カーネル関数が特定の次元に過度に影響されるのを防ぐ効果があります。本プログラムでは、以下の式で標準化を行います。

\[X_{std} = \frac{X - \mu}{\sigma}\]

ここで、\(X\) は元の記述子の値、\(\mu\) は訓練データの平均、\(\sigma\) は訓練データの標準偏差です。

5. 逐次最適化 (Sequential Optimization)

本プログラムの特筆すべき機能の一つは、ユーザーの入力に基づいてベイズ最適化のサイクルを繰り返す「逐次最適化」です。各サイクルでは、現在のガウス過程モデルが提案する最適な候補点が、擬似的に新たな訓練データとして扱われます。これにより、モデルは新しい情報を学習し、次のサイクルでより洗練された探索戦略を立てることができます。これは、限られた予算内で効率的に最適解に収束させる上で非常に強力なアプローチとなります。

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

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

  • numpy: 数値計算

  • pandas: データ構造とデータ解析

  • openpyxl: Excelファイル（.xlsx）の読み書き

  • matplotlib: グラフ描画

  • physbo: ベイズ最適化のコアライブラリ

physbo は特に重要なライブラリであり、互換性のある Python 3.6 環境でのインストールが推奨されています。

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

pip install numpy pandas openpyxl matplotlib physbo

もし、physbo のインストールで問題が発生する場合は、プログラムの冒頭にあるメッセージを参考に、Python 3.6 の仮想環境をアクティブにしてから再試行してください。

# 例: Anaconda環境でpy36というPython 3.6仮想環境がある場合
conda activate py36
pip install physbo

必要な入力ファイル

プログラムは、実験データまたはシミュレーション結果を含む Excel (.xlsx) または CSV (.csv) ファイルを読み込みます。

期待されるファイル形式とデータ構造:

  • ヘッダー行: 記述子（説明変数）と目的関数（応答変数）のラベルを含みます。

  • 制御コード: ヘッダーのラベルには、その列が目的関数であるか、またその最適化の方向（最大化、最小化、特定値への収束）を示すための制御コードを含めることができます。

    • =VALUE:ラベル: その列の値を VALUE に近づけるように最適化します（-(f(x) - VALUE)^2 を最大化する）。例えば =10.0:Target1 は Target1 を 10.0 に近づけます。

    • maxN:ラベル または toN:ラベル: その列の値を最大化します。N はパレート図にプロットする際の軸のインデックス（0から始まる、例えば max0:Target1 は最初のパレート軸）。N を省略すると自動的に割り当てられます。

    • minN:ラベル: その列の値を最小化します（内部的には -f(x) を最大化します）。N はパレート図にプロットする際の軸のインデックス。

    • -: その列を記述子からも目的関数からも除外します。

    • 上記以外のラベルは、デフォルトで記述子として扱われます。

  • 欠損値: 目的関数の列に欠損値（例: Excelの空白セルや NaN）がある行は、テストデータとして扱われ、ガウス過程モデルによる予測の対象となります。記述子の列に欠損値がある行は、解析から除外されます。

入力ファイル例 (data_simple.xlsx):

X1

X2

max:Y1

1

1

10

1

2

12

2

1

15

2

2

3

1

3

2

この例では、X1 と X2 が記述子、max:Y1 が最大化すべき目的関数です。最後の3行は Y1 の値が欠損しているため、ベイズ最適化によって Y1 の値が予測されます。

生成される出力ファイル

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

  1. search_result{i}.npz:

    • ファイル名: 入力ファイルと同じディレクトリに生成されます。{i} は目的関数のインデックス（1から始まる）に置き換えられます。例えば、data_simple.xlsx を入力とした場合、search_result1.npz のようになります。

    • 内容: PHYSBO ライブラリが生成したベイズ最適化の探索履歴が NumPy の圧縮形式で保存されます。これには、各ステップで選択されたアクションのインデックス、それに対応する目的関数値、獲得関数値などが含まれます。

  2. {filebody}-predict{i}.xlsx:

    • ファイル名: 入力ファイルと同じディレクトリに生成されます。{filebody} は入力ファイル名から拡張子を除いたものに、{i} は目的関数のインデックス（1から始まる）に置き換えられます。例えば、data_simple.xlsx を入力とした場合、data_simple-predict1.xlsx のようになります。

    • 内容: 全てのデータ点（訓練データとテストデータ）について、以下の情報を含む Excel ファイルです。

      • 元の記述子（入力ファイルの X1, X2 など）

      • 元の目的関数値（訓練データについては実測値、テストデータについては空白または NaN）

      • mean: ガウス過程モデルによって予測された目的関数の平均値

      • mean-std: 予測平均値から標準偏差を引いた値（予測信頼区間の下限）

      • mean+std: 予測平均値に標準偏差を加えた値（予測信頼区間の上限）

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

プログラムはコマンドライン引数を受け取ります。

python bayes_gp_plain_shimizu.py infile max_num_probes num_rand_basis score_mode interval standardize

  • infile: 入力データファイル名（例: data.xlsx）。

  • max_num_probes: ベイズ最適化の探索ステップ数。収束するまで繰り返すための回数を指定します。

  • num_rand_basis: ランダムな基底関数（Random Fourier Features）の数。-1を指定すると非近似ガウス過程（厳密解）になりますが、大規模データでは計算コストが高くなります。通常は、訓練データを再現できる程度の大きな数（例: 200）を指定します。

  • score_mode: 使用する獲得関数のモード。

    • EI: Expected Improvement (期待改善量)

    • PI: Probability of Improvement (改善確率)

    • TS: Thompson Sampling (トンプソンサンプリング)

  • interval: ハイパーパラメータを更新するサイクル数。0を指定すると初回にのみ更新します。

  • standardize: 記述子を標準化するかどうかのフラグ。

    • 0: 標準化しない

    • 1: 標準化する

例:

python bayes_gp_plain_shimizu.py data_simple.xlsx 10 200 EI 0 1

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

ここでは、上記の「必要な入力ファイル」で示した data_simple.xlsx を使用する例を考えます。

1. 入力ファイル data_simple.xlsx の作成

テキストエディタで以下の内容を記述し、data_simple.xlsx として保存してください。

X1,X2,max:Y1
1,1,10
1,2,12
2,1,15
2,2,
3,1,
3,2,

注意: 上記はCSV形式ですが、Excelで開いて保存し直せば .xlsx ファイルになります。また、Excelで直接上記の内容を入力しても構いません。Y1 の欠損値は空白セルとして入力します。

2. プログラムの実行

以下のコマンドでプログラムを実行します。

python bayes_gp_plain_shimizu.py data_simple.xlsx 5 200 EI 0 1

  • data_simple.xlsx: 入力ファイル名。

  • 5: max_num_probes (探索ステップ数)。

  • 200: num_rand_basis (ランダム基底関数の数)。

  • EI: score_mode (獲得関数としてExpected Improvementを使用)。

  • 0: interval (ハイパーパラメータの更新間隔。0は初回のみ更新)。

  • 1: standardize (記述子を標準化する)。

3. 実行結果の説明

プログラムはまず、コンソールにデータの読み込み状況、記述子、目的関数の情報、そして標準化の有無などを表示します。

次に、ベイズ探索が開始され、以下のようなプロンプトが表示されます。

繰り返し回数を入力>>>

ここで、最適化のサイクルを繰り返したい回数を入力します。例えば 3 と入力すると、3回の逐次最適化サイクルが実行されます。各サイクルの後には、その時点でのベストな候補点に関する情報がコンソールに表示され、その候補点が次のサイクルの訓練データに追加されます。

その後、ガウス過程モデルの予測結果と獲得関数値に基づいて、以下のグラフウィンドウが順次開いて表示されます。

  • contour: mean (予測平均値の等高線図): 記述子 X1, X2 を軸にとり、目的関数 Y1 の予測平均値の等高線が表示されます。既存の訓練データ点は黄色い星印で示され、その大きさは目的関数値の大小に対応します。

  • contour: std (予測標準偏差の等高線図): 同様に、目的関数 Y1 の予測標準偏差の等高線が表示されます。不確実性の高い領域が視覚的に確認できます。

  • prediction and aquisition functions (予測値と獲得関数のプロット): 各データ点（インデックス）を横軸にとり、Y1 の予測平均値（黒線）、その信頼区間（青い帯）、訓練データ点の実測値（青い丸）、そして獲得関数 EI の値（赤い破線）がプロットされます。EI が最も高い点が次の探索候補点として星印で示されます。

これらのグラフは、プログラムが終了するまで表示され続けます。グラフをクリックすると、その点に最も近いデータの詳細情報がコンソールに表示されます。

コンソール出力例 (一部):

Read data from [data_simple.xlsx]
  descriptors: ['X1', 'X2']
  objective functions:
    i=0 Pareto plot index=0 label=Y1 mode=max

# of all data: 6
# of training data: 3
# of descriptors        : 2
# of objective variables: 1
  Objective variables: ['Y1']
# of descriptors: 2
  Descriptors: ['X1', 'X2']
  Descriptors for plot: {'x': 'X1', 'y': 'X2'}

Make policy:
  Training data indexes:
[0 1 2]
    descriptors from X_train:
[[1. 1.]
 [1. 2.]
 [2. 1.]]
    descriptors from X_all:
[[1. 1.]
 [1. 2.]
 [2. 1.]]
    objective values:
[[10.]
 [12.]
 [15.]]

Standardize descriptors:
       X1: X(standardized)[i] = (X[i] -       1.333) /      0.4714
       X2: X(standardized)[i] = (X[i] -       1.333) /      0.4714

Objective #1

Start Bayes search:

Bayse search:
show_search_results
   num_probe best_action best_fx action_X
           0           2   15.00  [2.0, 1.0]
           1           2   15.00  [2.0, 1.0]
           2           2   15.00  [2.0, 1.0]
           3           2   15.00  [2.0, 1.0]
           4           2   15.00  [2.0, 1.0]

  Best candidate    : 2 [2. 1.] 15.0
Save seach result to [search_result1.npz]

Save predictions to [data_simple-predict1.xlsx]
X_all.T=
[[1. 1. 2. 2. 3. 3.]
 [1. 2. 1. 2. 1. 2.]]

Plot
  ntargets=1 in graphs 1 x 1
  ndescriptors=2
  indexes for 2D contour plot: idx=0 for x-axis, 1 for y-axis

繰り返し回数を入力>>> 3
[15.0]
[[1. 1.]
 [1. 2.]
 [2. 1.]
 [2. 1.]]
[[10.]
 [12.]
 [15.]
 [15.]]
[0 1 2 2]
[15.0]
[[1. 1.]
 [1. 2.]
 [2. 1.]
 [2. 1.]
 [2. 1.]]
[[10.]
 [12.]
 [15.]
 [15.]
 [15.]]
[0 1 2 2 2]
[15.0]
[[1. 1.]
 [1. 2.]
 [2. 1.]
 [2. 1.]
 [2. 1.]
 [2. 1.]]
[[10.]
 [12.]
 [15.]
 [15.]
 [15.]
 [15.]]
[0 1 2 2 2 2]

====================================================================================
Please cite designated referneces for PHYSBO. See
              Citation: https://issp-center-dev.github.io/PHYSBO/manual/master/en/introduction.html
====================================================================================

Press ENTER to terminate

この例では、max_num_probes = 5 ですが、show_search_results ではベストアクションは常にインデックス2 ([2.0, 1.0]) であり、best_fx = 15.0 のようです。これは、初期の訓練データで既に最適な点が学習されているか、あるいは探索空間が単純で速やかに収束したことを示唆します。繰り返し回数を 3 と入力したことで、同じ [2.0, 1.0] のデータが3回訓練データに追加され、モデルが更新されています。

生成されるファイル:

  • search_result1.npz

  • data_simple-predict1.xlsx

data_simple-predict1.xlsx を開くと、X1, X2 の列に加えて Y1、mean、mean-std、mean+std の列が追加されており、欠損していた Y1 の箇所に予測値が記入されていることを確認できます。