ガウス過程

PHYSBOではガウス過程回帰を実行しながらベイズ最適化を行なっています。

そのため、学習データが与えられた際にガウス過程回帰を実行することもでき、学習済みモデルを利用したテストデータの予測も行うことができます。

ここでは、その手順について紹介します。

探索候補データの準備

本チュートリアルでは例として、Cuの安定した界面構造の探索問題を扱います。目的関数の評価にあたる構造緩和計算には、実際には1回あたり数時間といったオーダーの時間を要しますが、本チュートリアルでは既に評価済みの値を使用します。問題設定については、以下の文献を参照してください。 S. Kiyohara, H. Oda, K. Tsuda and T. Mizoguchi, “Acceleration of stable interface structure searching using a kriging approach”, Jpn. J. Appl. Phys. 55, 045502 (2016).

データセットファイル s5-210.csv を data ディレクトリ以下に保存し、次のように読み出します。

[1]:

import numpy as np

import matplotlib.pyplot as plt

import physbo


def load_data():
    A =  np.asarray(np.loadtxt('data/s5-210.csv',skiprows=1, delimiter=',') )
    X = A[:,0:3]
    t  = -A[:,3]
    return X, t

X, t = load_data()
X = physbo.misc.centering( X )

学習データの定義

対象データのうち、ランダムに選んだ1割をトレーニングデータとして利用し、別のランダムに選んだ1割をテストデータとして利用します。

[2]:

N = len(t)
Ntrain = int(N*0.1)
Ntest = min(int(N*0.1), N-Ntrain)

id_all   = np.random.choice(N, N, replace=False)
id_train  = id_all[0:Ntrain]
id_test = id_all[Ntrain:Ntrain+Ntest]

X_train = X[id_train]
X_test = X[id_test]

t_train = t[id_train]
t_test = t[id_test]

print("Ntrain =", Ntrain)
print("Ntest =", Ntest)

Ntrain = 1798
Ntest = 1798

ガウス過程の学習と予測

以下のプロセスでガウス過程を学習し、テストデータの予測を行います。

ガウス過程のモデルを生成します。
X_train（学習データのパラメータ）, t_train（学習データの目的関数値）を用いてモデルを学習します。
学習されたモデルを用いてテストデータ(X_test)に対する予測を実行します。

共分散の定義(ガウシアン)

[3]:

cov = physbo.gp.cov.Gauss( X_train.shape[1],ard = False )

平均の定義

[4]:

mean = physbo.gp.mean.Const()

尤度関数の定義(ガウシアン)

[5]:

lik = physbo.gp.lik.Gauss()

ガウス過程モデルの生成

[6]:

gp = physbo.gp.Model(lik=lik,mean=mean,cov=cov)
config = physbo.misc.set_config()

ガウス過程モデルを学習

[7]:

gp.fit(X_train, t_train, config)

Start the initial hyper parameter searching ...
Done

Start the hyper parameter learning ...
0 -th epoch marginal likelihood 13760.853923583709
50 -th epoch marginal likelihood 6031.859927457567
100 -th epoch marginal likelihood 3418.8688786331622
150 -th epoch marginal likelihood 1713.752090407802
200 -th epoch marginal likelihood 329.20339957667966
250 -th epoch marginal likelihood -482.55368524374444
300 -th epoch marginal likelihood -982.9079141420298
350 -th epoch marginal likelihood -1140.7031705370855
400 -th epoch marginal likelihood -1367.476428948909
450 -th epoch marginal likelihood -1508.1346001789386
500 -th epoch marginal likelihood -1566.2691297825002
Done

学習されたガウス過程におけるパラメタを出力

[8]:

gp.print_params()



likelihood parameter =   [-3.23999709]
mean parameter in GP prior:  [-1.0613511]
covariance parameter in GP prior:  [-1.38290386 -2.25189235]

テストデータの平均値（予測値）および分散を計算

[9]:

gp.prepare(X_train, t_train)
fmean = gp.get_post_fmean(X_train, X_test)
fcov = gp.get_post_fcov(X_train, X_test)

予測の結果

[10]:

fmean

[10]:

array([-1.03065227, -1.0563664 , -1.04071129, ..., -0.97478745,
       -0.98638275, -1.09557148], shape=(1798,))

分散の結果

[11]:

fcov

[11]:

array([0.00047079, 0.00072114, 0.00054597, ..., 0.00067217, 0.00105654,
       0.00133083], shape=(1798,))

予測値の平均二乗誤差の出力

[12]:

np.mean((fmean-t_test)**2)

[12]:

np.float64(0.007645971280859177)

得られた回帰モデルについて、どの特徴量がどのくらい重要なのかを調べるための簡単な指標として permutation importance (PI)があります。 PHYSBOでは get_permutation_importance 関数で計算できます。

[13]:

pi_mean, pi_std = gp.get_permutation_importance(X_train, t_train, n_perm=100)

[14]:

features = list(range(len(pi_mean)))

plt.figure(figsize=(8, 5))
plt.barh(
    features,
    pi_mean,
    xerr=pi_std,
)
plt.gca().invert_yaxis()
plt.yticks(features)
plt.xlabel("Permutation Importance")
plt.ylabel("Features")
plt.tight_layout()

save permutation_importance.pdf

../_images/notebook_tutorial_Gaussian_process_26_1.png

このグラフにおいて、棒グラフはPIの平均値を、線分は標準偏差を示しており、0番の特徴量と比べて1,2番の特徴量のほうが重要そうだ、ということが読み取れます。

gp.model だけではなく, policy も get_permutation_importance 関数を持っています。使い方は policy.get_permutation_importance(n_perm) です。 policy では保存してある学習済みデータに対してPIを計算するため、 model とは違い、改めて渡す必要はありません。

訓練済みモデルによる予測

学習済みモデルのパラメタをgp_paramsとして読み出し、これを用いた予測を行います。

gp_paramsおよび学習データ(X_train, t_train)を記憶しておくことで、訓練済みモデルによる予測が可能となります。

学習されたパラメタを準備(学習の直後に行う必要あり)

[15]:

#学習したパラメタを１次元配列として準備
gp_params =  np.append(np.append(gp.lik.params, gp.prior.mean.params), gp.prior.cov.params)

gp_params

[15]:

array([-3.23999709, -1.0613511 , -1.38290386, -2.25189235])

学習に利用したモデルと同様のモデルをgpとして準備

[16]:

#共分散の定義 (ガウシアン)
cov = physbo.gp.cov.Gauss( X_train.shape[1],ard = False )

#平均の定義
mean = physbo.gp.mean.Const()

#尤度関数の定義 (ガウシアン)
lik = physbo.gp.lik.Gauss()

#ガウス過程モデルの生成
gp = physbo.gp.Model(lik=lik,mean=mean,cov=cov)

学習済みのパラメタをモデルに入力し予測を実行

[17]:

#学習済みのパラメタをガウス過程に入力
gp.set_params(gp_params)


#テストデータの平均値（予測値）および分散を計算
gp.prepare(X_train, t_train)
fmean = gp.get_post_fmean(X_train, X_test)
fcov = gp.get_post_fcov(X_train, X_test)

予測の結果

[18]:

fmean

[18]:

array([-1.03065227, -1.0563664 , -1.04071129, ..., -0.97478745,
       -0.98638275, -1.09557148], shape=(1798,))

分散の結果

[19]:

fcov

[19]:

array([0.00047079, 0.00072114, 0.00054597, ..., 0.00067217, 0.00105654,
       0.00133083], shape=(1798,))

予測値の平均二乗誤差の出力

[20]:

np.mean((fmean-t_test)**2)

[20]:

np.float64(0.007645971280859177)

(注) 上の例では事前に登録されているXと同じものを利用して予測を行いました。学習済みのモデルを利用してXに含まれていないパラメータ X_new に対して予測をしたい場合には、学習モデルで使用したデータXの平均(X_{mean})と標準偏差(X_{std})を求めていただいたうえで、 X_{new} = (X_{new} - X_{mean}) / X_{std} の変形を行うことで予測を行うことができます。また、渡す際のデータ形式はndarray形式になっています。そのため、X_{new}が一つのデータの場合には事前に変換する必要があります。例えば、X_{new}が実数である場合には、 X_new = np.array(X_new).reshape(1) などとして変換する必要があります。