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

このドキュメントは、Pythonスクリプト bayes_gp_plain.py の技術的な側面について解説します。

プログラムの動作

bayes_gp_plain.py は、PHYSBOライブラリを使用して、ガウス過程に基づくベイズ最適化を実行し、最適な材料設計パラメータやプロセス条件を効率的に探索することを目的としたプログラムです。

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

データ読み込み: Excel (.xlsx) または CSV (.csv) 形式の入力ファイルから、記述子 (特徴量) と目的関数 (ターゲット) のデータを読み込みます。
目的関数設定: 入力ファイルのヘッダー行に記述された特殊なプレフィックス (例: max:, min:, =value:) に基づいて、目的関数の最適化モード（最大化、最小化、特定値への適合）を自動的に設定します。
データ前処理: 欠損値 (NaN) を含む行を処理し、学習データと予測対象のテストデータを識別します。オプションで記述子の標準化を行います。
ベイズ探索の実行: PHYSBOライブラリを利用し、指定された獲得関数 (EI, PI, TS) と試行回数でベイズ最適化を実行します。ランダム基底近似を用いた高速化も可能です。
結果の保存: ベイズ探索の履歴 (.npz ファイル) と、各データ点における目的関数の予測平均値、標準偏差、信頼区間を含む詳細な予測結果 (.xlsx ファイル) を出力します。
結果の可視化: Matplotlib を用いて、予測された目的関数の平均値と標準偏差の等高線図、目的関数が複数ある場合のパレート図、およびデータ点ごとの予測値と獲得関数の推移をグラフ表示します。
インタラクティブな情報表示: 生成されたグラフ上の任意の点をクリックすることで、その点に最も近いデータ点の記述子、予測値、獲得関数値などの詳細情報をコンソールに出力します。

このプログラムは、多次元の設計空間において、実験やシミュレーションの評価コストが高い場合に、効率的に最適解を見つけ出すための強力なツールを提供します。

原理

bayes_gp_plain.py は、ベイズ最適化の主要な手法であるガウス過程 (Gaussian Process, GP) を利用しています。

ベイズ最適化の概要: ベイズ最適化は、未知の目的関数 $f(x)$ を最小化（または最大化）する $x$ を探索する際に、評価コストが高い場合に特に有効な最適化手法です。これは、以下の2つの主要なコンポーネントで構成されます。

代理モデル (Surrogate Model): これまでに観測されたデータ点 $\{(x_1, f(x_1)), \dots, (x_t, f(x_t))\}$ を用いて、目的関数 $f(x)$ の挙動を近似するモデル。本プログラムではガウス過程が用いられます。
獲得関数 (Acquisition Function): 代理モデルの予測（平均値と不確実性）を利用して、次に評価すべき最も有望な $x$ を決定するための基準。

ガウス過程 (Gaussian Process, GP): GPは、関数に対する確率分布を定義するノンパラメトリックなモデルです。これは、関数の任意の有限個の点における値が、同時ガウス分布に従うと仮定します。GPは「平均関数」と「共分散関数（カーネル関数）」によって完全に定義されます。

平均関数 $m(x)$: 関数の期待値。通常は0または線形関数が使われます。
共分散関数 $k(x, x')$: 2つの入力 $x$ と $x'$ の関数値の類似度を定義します。これにより、近い入力の関数値は似ているという事前知識をモデルに組み込むことができます。

観測データ $\{(x_1, f(x_1)), \dots, (x_t, f(x_t))\}$ が与えられると、GPはベイズの定理に基づいて事後分布を更新します。これにより、未観測点 $x^*$ における目的関数の予測平均値 $\mu(x^*)$ と予測分散 $\sigma^2(x^*)$ を得ることができます。

\[ f(x^*) \sim \mathcal{N}(\mu(x^*), \sigma^2(x^*)) \]

獲得関数 (Acquisition Functions): 予測平均値 $\mu(x^*)$ と予測分散 $\sigma^2(x^*)$ を利用して、次の探索点を決定します。本プログラムでは、以下の獲得関数モードをサポートしています。

Expected Improvement (EI): 改善の期待値。これまでの最適値 $f_{best}$ よりも改善するであろう期待値が高い点を優先します。 $$ EI(x) = E[\max(0, f(x) - f_{best})] $$
Probability of Improvement (PI): 改善確率。これまでの最適値 $f_{best}$ よりも改善する確率が高い点を優先します。 $$ PI(x) = P(f(x) > f_{best}) $$
Thompson Sampling (TS): GPの事後分布から関数をサンプリングし、そのサンプリングされた関数の最大値（または最小値）を与える点を次の評価点とします。

目的関数の変換: プログラムでは、入力データのヘッダーに基づいて、目的関数をベイズ最適化が扱いやすい形に変換します。

最小化 (min モード): 目的関数 $f(x)$ を最大化問題に変換するため、$-f(x)$ を新しい目的関数として扱います。
特定値への適合 (value モード): 特定の値 $V$ に近づけることを目標とする場合、目的関数を $-(f(x) - V)^2$ に変換し、この新しい目的関数の最大化を目指します。これにより、$V$ に近づくほど値が大きくなります。

これらの原理に基づいて、bayes_gp_plain.py は、与えられたデータからガウス過程モデルを構築し、獲得関数を用いて次に評価すべき最適な点を繰り返し探索することで、目的関数の最適化を効率的に進めます。

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

bayes_gp_plain.py の実行には、以下の非標準ライブラリが必要です。

numpy: 高度な数値計算をサポートします。
pandas: データ構造とデータ解析ツールを提供します。ExcelやCSVファイルの読み書きに利用されます。
openpyxl: pandas がExcelファイルを読み込む際に内部的に使用されます。
matplotlib: グラフの描画に使用されます。
physbo: ガウス過程に基づくベイズ最適化のコア機能を提供します。

これらのライブラリは、Pythonのパッケージマネージャーである pip を使用してインストールできます。

インストール手順:

推奨環境: physbo が特定のPythonバージョン（例: Python 3.6）での使用を推奨している場合があります。仮想環境の使用をお勧めします。
```
# 例: Python 3.6 の仮想環境を作成してアクティブ化する場合 (Anaconda/Minicondaを使用)
conda create -n py36 python=3.6
conda activate py36
```
(一般的な環境であれば、仮想環境の作成は必須ではありませんが、推奨されます。)
ライブラリのインストール: アクティブ化した仮想環境、または通常のPython環境で以下のコマンドを実行します。
```
pip install numpy pandas openpyxl matplotlib physbo
```
Pythonの実行ファイルが python3 の場合は、pip3 を使用します。
```
pip3 install numpy pandas openpyxl matplotlib physbo
```
または、Pythonモジュールとして実行します。
```
python -m pip install numpy pandas openpyxl matplotlib physbo
```

必要な入力ファイル

プログラムは、構造化された実験データまたはシミュレーション結果を含むExcel (.xlsx) または CSV (.csv) ファイルを入力として期待します。

ファイル形式とデータ構造:

形式: Microsoft Excel Workbook (.xlsx) または Comma Separated Values (.csv)
ヘッダー行:
- 各列の最初の行（ヘッダー）は、その列が記述子であるか目的関数であるかを定義します。
- 目的関数: 以下のいずれかのプレフィックスで開始します。
  - maxN:目的関数名: その列を最大化する目的関数として扱います。N はオプションの整数で、パレート図のX軸/Y軸の優先順位を決定します (0 または 1 が指定された場合、それぞれ第一・第二目的関数としてプロットされます)。N が省略された場合も自動的に割り当てられます。
  - minN:目的関数名: その列を最小化する目的関数として扱います。内部的には値が反転されます。N の扱いは maxN: と同様です。
  - =VALUE:目的関数名: その列の値を VALUE に近づける目的関数として扱います。内部的には -(値 - VALUE)^2 が最大化されます。
  - 例: max:Yield, min1:Cost, =100:Purity
- 記述子: 目的関数のプレフィックスを持たない列は、記述子として扱われます。
  - 記述子名に x:, y:, z: といったプレフィックスをつけることで、グラフ描画時の軸を指定できます (例: x:Temperature, y:Pressure)。
- 除外: - で始まるヘッダーを持つ列は、処理から除外されます。
  - 例: -Comment
データ行:
- ヘッダー行に続く各行は、1つの実験またはシミュレーションの結果を表します。
- 学習データ: 目的関数の列に数値が入力されている行は、ガウス過程モデルの学習データとして使用されます。
- テストデータ (予測対象): 目的関数の列が空 (NaN または空白セル) の行は、ベイズ最適化の候補点 (テストデータ) として扱われ、予測の対象となります。記述子の値は入力されている必要があります。

例 (data_simple.xlsx の一部を想定):

x:Temp	y:Press	max:Yield	max1:Purity
300	5	85	92
320	6	88	95
310	5.5
290	4.8	82	90
330	6.2

この例では、Temp と Press が記述子、Yield と Purity が目的関数です。3行目と5行目は Yield と Purity の値が空であるため、ベイズ最適化の候補点（テストデータ）として扱われます。

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

プログラムの実行後、以下のファイルがカレントディレクトリ、または入力ファイルと同じディレクトリに生成されます。

PHYSBO履歴ファイル (.npz):
- ファイル名: {入力ファイル名}-save{i}.npz
  - {入力ファイル名}: 元の入力ファイルのファイル本体名 (拡張子を除く)。
  - {i}: 目的関数のインデックス (1から開始)。例えば、入力ファイルに2つの目的関数がある場合、data_simple-save1.npz と data_simple-save2.npz が生成されます。
- 内容: PHYSBOライブラリによって生成されたベイズ探索の履歴オブジェクトがNumPyの .npz 形式で保存されます。これには、探索の全シーケンス、各ステップで選択された行動、対応する目的関数値、およびモデルのハイパーパラメータに関する情報が含まれます。これは、後から探索プロセスを分析したり、再開したりするために使用できます。
予測結果ファイル (.xlsx):
- ファイル名: {入力ファイル名}-predict{i}.xlsx
  - {入力ファイル名}: 元の入力ファイルのファイル本体名。
  - {i}: 目的関数のインデックス (1から開始)。履歴ファイルと同様に、目的関数の数だけファイルが生成されます。
- 内容: ガウス過程モデルによって予測された結果を含むExcelファイルです。各行は入力データの1点に対応し、以下の情報を含みます。
  - 全記述子の値: 入力ファイルに存在するすべての記述子の値。
  - 目的関数の元の値: 学習データ点については、入力ファイルから読み込まれた目的関数の値。テストデータ点では空欄またはNaN。
  - mean: ガウス過程モデルによって予測された目的関数の平均値。
  - mean-std: 予測平均値から予測標準偏差を引いた値（信頼区間の下限）。
  - mean+std: 予測平均値に予測標準偏差を加えた値（信頼区間の上限）。

これらの出力ファイルは、ベイズ最適化の結果を詳細に分析し、新たな実験計画を立てるための基礎データとして利用できます。

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

bayes_gp_plain.py は、以下の書式でコマンドラインから実行します。

python bayes_gp_plain.py infile max_num_probes num_rand_basis score_mode interval standardize

各引数の説明は以下の通りです。

infile:
- 説明: 入力データファイル名。Excel (.xlsx) または CSV (.csv) 形式のファイルを指定します。
- デフォルト値: data_simple.xlsx (コード内)
max_num_probes:
- 説明: ベイズ探索を実行する最大試行回数。モデルが収束に達するまで、この数を増やしてください。
- デフォルト値: 1 (コード内)
num_rand_basis:
- 説明: ガウス過程の計算に使用するランダム基底近似の基底数。学習データを忠実に再現するには大きな数を指定します。-1 を指定するとランダム基底近似を使用せず、正確な（しかし計算コストの高い）GP計算を行います。
- デフォルト値: 200 (コード内)
score_mode:
- 説明: 獲得関数 (Acquisition Function) のモードを指定します。以下のいずれかの値を設定できます。
  - EI: Expected Improvement (改善の期待値)
  - PI: Probability of Improvement (改善確率)
  - TS: Thompson Sampling (トンプソンサンプリング)
- デフォルト値: EI (コード内)
interval:
- 説明: ガウス過程のハイパーパラメータを更新するサイクル数。0 を指定すると、ハイパーパラメータは初期値から更新されません。1 を指定すると各プローブで更新されます。
- デフォルト値: 0 (コード内)
standardize:
- 説明: 記述子データを標準化するかどうかのフラグです。
  - 0: 標準化を行いません。
  - 1: 標準化を行います（平均0、標準偏差1）。
- デフォルト値: 1 (コード内)

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

ここでは、data_simple.xlsx という入力ファイル（上記の「必要な入力ファイル」セクションの例を参考にしてください）を使用して、具体的な実行例とその結果を説明します。

実行コマンド:

入力ファイル data_simple.xlsx を使用し、ベイズ最適化を最大10回試行します。ランダム基底近似には200の基底を使用し、獲得関数には EI (Expected Improvement) を用います。ハイパーパラメータは更新せず (interval 0)、記述子は標準化します (standardize 1)。

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

実行結果の説明:

初期化とパラメータ表示:
- プログラムが起動し、コマンドライン引数で与えられた、またはデフォルトで設定されたパラメータ（入力ファイル名、試行回数、基底数、獲得関数モード、インターバル、標準化フラグなど）がコンソールに表示されます。
- data_simple.xlsx からのデータ読み込みが開始され、記述子のリストと目的関数の設定（ラベル、モード、パレートプロットインデックス）が示されます。
- 学習データの数、記述子の数、目的変数の数などが表示され、記述子が標準化される場合はその詳細も出力されます。
ベイズ探索の開始:
- 指定されたmax_num_probes (ここでは10) に基づいて、PHYSBOによるベイズ探索が開始されます。
- 各探索ステップにおいて、ガウス過程モデルが更新され、獲得関数 (EI) の値が計算されます。最も高い獲得関数値を持つ点が次の探索候補として選択されます。
- コンソールには、探索の進行状況と、各プローブで選択された行動やその結果（もし利用可能であれば）の概要が表示されます。
結果の保存:
- ベイズ探索が完了すると、以下のファイルが生成されます（入力ファイルに複数の目的関数がある場合、{i} の部分が目的関数の数だけ増えます）。
  - data_simple-save1.npz: 最初の目的関数に関するPHYSBOの探索履歴が保存されます。
  - data_simple-predict1.xlsx: 最初の目的関数に関する予測結果（全記述子、元の目的関数値、予測平均値、平均値±標準偏差）がExcel形式で保存されます。
グラフの表示:
- Matplotlibによって複数のグラフィカルウィンドウが表示されます。
  - contour: mean ウィンドウ: 選択された2つの記述子 (例: x:Temp, y:Press) に対する、目的関数（例: max:Yield）の予測平均値の等高線図が表示されます。学習データ点は星印で示されます。
  - contour: std ウィンドウ: 同様に、目的関数の予測標準偏差の等高線図が表示されます。不確実性の高い領域が視覚化されます。
  - Pareto ウィンドウ (目的関数が2つ以上の場合): 2つの目的関数（例: max:Yield と max1:Purity）の予測平均値と、既存の学習データ点の散布図が表示され、パレートフロントの探索状況を示します。
  - prediction and aquisition functions ウィンドウ: 各データ点（インデックス）に対する、学習データ点、予測平均値と信頼区間（平均値±標準偏差）、および獲得関数（EI）の値がプロットされます。これにより、どの領域で不確実性が高く、どの領域が探索に有望であるかを確認できます。
インタラクティブな操作:
- 表示されているグラフの任意の点をマウスでクリックすると、その点に最も近いデータ点の詳細情報がコンソールに出力されます。これには、記述子の値、予測された目的関数の平均値とその標準偏差、および獲得関数の値が含まれます。等高線図では、クリックされた位置に最も近い学習データ点とその目的関数値も表示されます。
プログラムの終了:
- すべての処理が完了し、グラフが表示された後、「Press ENTER to terminate」というプロンプトがコンソールに表示されます。Enterキーを押すと、プログラムは終了し、グラフウィンドウも閉じられます。

この一連の動作により、ユーザーは手元のデータに基づいて効率的に最適な条件を探索し、その結果を視覚的に理解することができます。