0. この論文を一言で言うと

ベイズ最適化を用いて、結晶界面の原子構造の決定において有望な結果を提示するオープンソースライブラリCOMBO(COMmon Bayesian Optimization library)を開発した.
COMBO github

##1. Introduction
実験の設計は, 実験結果から次の手順を適用する繰り返し手順であり,多くの場合は実験はシミュレーションで代用される.
数学的には, ブラックボックス関数の最適化問題として表される.
ベイズ最適化は, ブラックボックス最適化に対する確立された技術である.

材料科学では, 少なくとも４つのベイズ最適化のアプリケーション（弾性特性, 融点, 格子熱伝導率, 粒界界面の設計）が存在するが, これらのアプリケーションで扱うことのできる学習データの数が少ない（数百程度）という問題がある.
今後大規模な問題を扱うことを考えると, 学習データの数は数万程度まで増えるが, 既存のパッケージ（scikit-learnなど）では, 計算時間が遅いことが知られている.

我々は, 最適に効率化した最新の機械学習アルゴリズムを開発し、オープンソースライブラリCOMBOを作った.
COMBOは候補数の増加に対し, 計算量が線形でしか増えないため、大規模問題の扱いに向いている.

##2. Technical Features and Details
###2.0 はじめに
候補点のセットをd次元のベクトルとして表す.
$\boldsymbol{x_1}, . . . , \boldsymbol{x_m}\in \mathbb{R}^d$
実験により, 候補点でのブラックボックス関数の値を得ることができる.
我々は, できるだけ少ない実験で, その値を最大にする候補点を同定したい.

連続した実験設計において, 候補点 $\boldsymbol{x_1}, ... , \boldsymbol{x_n}$

から値 $y_1, . . . , y_n$ が既に測定されていることを想定する.

学習データ $D=(\boldsymbol{x_i}, y_i)_{i=1}^n$から学んだ予測モデルに基づいて,

アルゴリズムは, 残りの$m-n$のデータからどの1つを適用すべきかを決定する.
この過程は, 予め決めた実験回数が終わるまで繰り返される.

###2.1 Thompson Sampling
トンプソンサンプリングは, 連続的な実験設計に対する人気のアルゴリズムである. トンプソンサンプリングの基本的な考え方は, 確率マッチングであり, 最適である確率に従って候補点を選択するアルゴリズムである.

トンプソンサンプリングを容易にするためには,予測モデルが必要であり, ここでは, ベイジアン線形回帰モデルを採用する.
$$ y = \boldsymbol{w}^{{\rm T}}\phi(\boldsymbol{x}) + \varepsilon$$

$・\boldsymbol{x}\in\mathbb{R}^d　　　: 入力ベクトル(候補点)$
$・\phi:\mathbb{R}^d→\mathbb{R}^l :特徴ベクトル$
$・\boldsymbol{w} \in \mathbb{R}^l　　　:重み行列$
$・\varepsilon　　　　　　:N(0, \sigma^2)に従うノイズ$

$\Phi$ を $i$ 列目が $\phi(\boldsymbol{x_i})$ に相当する $l×n$ 行列とすると,事後分布 $\boldsymbol{w}{\hspace 1mm}|{\hspace 1mm} D$ は以下となる.
$$\boldsymbol{w}{\hspace 1mm}|{\hspace 1mm}D \sim N(\mu,\sigma^2)$$なお, $\mu=(ΦΦ +\sigma^2I)^{-1} Φ_\boldsymbol{y}, Σ=\sigma^2(ΦΦ^T +\sigma^2I)^{-1}$である.

$候補点\boldsymbol{x_i}$に対する予測値は, $\boldsymbol{w}^{{\rm T}}\phi(\boldsymbol{x})$である.
$領域\boldsymbol{w}$は, 以下の最大化で定義される.
$$W_i = \lbrace w \in R^l {\hspace 1mm}|{\hspace 1mm} w^{\rm T}\phi(\boldsymbol{x_i}) = \underset{j}{max}{\hspace 1mm} \boldsymbol{w}^{\rm T}\phi(\boldsymbol{x_j}) \rbrace$$トンプソンサンプリングは、最適である確率$p_i=P(\boldsymbol{w} \in \boldsymbol{W_i}{\hspace 1mm}|{\hspace 1mm}D)$に従って候補点を選択する.
次のように$p_i$を計算せずに実行できる.
・$P(\boldsymbol{w}{\hspace 1mm}|{\hspace 1mm}D)$からベクトル $\boldsymbol{s}$ をサンプリングする.
・$\boldsymbol{s}$に関して最大スコアを持つ候補点を以下の式で決定する, $$i^∗ = \underset{i}{argmin}{\hspace 1mm} \boldsymbol{s}^{\rm T}\boldsymbol{x_i}.$$ベイジアン線形モデルのトンプソンサンプリングは, 候補点を評価する計算時間が$O(l)$でかかるため, 特に効率的である.

###2.2 Quick Sampling
多次元ガウス分布からのサンプリングは，共分散行列の三角分解を必要とする．
我々の場合, 事後分布$P(\boldsymbol{w}{\hspace 1mm}|{\hspace 1mm}D)$からベクトル$\boldsymbol{s}$をサンプリングする必要がある.
$\boldsymbol{A}$を以下の式で定義する:
$$ A= \frac{1}{\sigma^2}\Phi\Phi^{\rm T}+I$$すると, 事後分布は以下の式で表される.
$$\boldsymbol{w}{\hspace 1mm}|{\hspace 1mm}D \sim N(\frac{1}{\sigma^2}\boldsymbol{A}^{-1})$$したがって，トンプソンサンプリングの各ステップで$A^{-1}$の三角分解が必要となる．実験が行われた後，新しい学習例$(\boldsymbol{x}',y')$が得られ，行列 $\boldsymbol{A}$ は次のように更新される.
$$ \boldsymbol{A'} = \boldsymbol{A} + \frac{1}{\sigma^2}\phi(\boldsymbol{x'})\phi(\boldsymbol{x'})^{\rm T}$$1ランク変更された行列のコレスキー分解は、元の行列から $O(l^2)$ 時間で計算される(すなわち, fast-update)のに対し, ゼロから計算すると $O(l^3)$ 時間がかかることはよく知られている. そこで, 我々は $\boldsymbol{A}$ のコレスキー分解Lをプロセス中ずっと維持し(すなわち $A = L^{\rm T}L$), 新しい例が追加されるたびにLに高速更新を適用する.
サンプル$\boldsymbol{s} $は $\boldsymbol{\mu}＋s_0$として得られ, ここで$s_0$は$N(0,\boldsymbol{A}^{-1})$からの試料である. (1)に $\boldsymbol{A} = L^{\rm T}L$ を代入すると,平均ベクトル $\boldsymbol{\mu}$ は次式の解として記述される.
$$L^{\rm T}L{\hspace 0.5mm}\boldsymbol{\mu}=\frac{1}{\sigma^2}\Phi{\hspace 0.5mm}\boldsymbol{y}$$
$\boldsymbol{\mu}$ は,1回目は$L^{\rm T}$に対して、2回目は$L$に対して逆置換をかけることで解が得られる. また, $z \sim N(0,I)$をサンプリングし, 1次方程式 $\boldsymbol{z} = L{\hspace 0.5mm}\boldsymbol{s_0}$を逆置換で解くことで$s_0$が得られる. サンプリング処理は$O(l^2)$時間で行われる.

###2.3 Random Feature Map
特徴マップ$\phi$の選択はベイズ最適化の性能に大きく影響する. COMBOでは, ガウスカーネルによる写像に近似したランダムな特徴量マップを採用している．ガウスカーネルを $k(∆)=exp(-∥∆∥^2/2)$と表現しましょう。ボヒナーの定理により, 次のように書き換えらる.
$$ k(\boldsymbol{x}−\boldsymbol{x′}) = \int\exp(j\boldsymbol{ω}^{\rm T}(\boldsymbol{x}−\boldsymbol{x′}))p(\boldsymbol{ω})d\boldsymbol{ω} $$$$ p(ω) = (2π)^{−d/2} exp(−∥\boldsymbol{\omega}∥^2/2) $$ $z_{\boldsymbol{ω},b}(\boldsymbol{x}) = \sqrt{2}cos(ω^{\rm T}\boldsymbol{x}+b)$と定義する. $\boldsymbol{\omega}$を$p(\omega)$から引き, $b$を$[0, 2\pi]$から一様に引くと, $E[z_{\boldsymbol{ω},b}(\boldsymbol{x})z_{ω,b}(\boldsymbol{x′})] = k(\boldsymbol{x}-\boldsymbol{x′})$となる. COMBOでは, 特徴量マップは$l$個のランダムサンプル$(w_i, b_i)^{l}_{i=1}$を用いて以下で定義される.

$$ \phi(x)=(z_{\omega_1,b_1}(\boldsymbol{x}),...,z_{\omega_l,b_l}(\boldsymbol{x}))^{\rm T} $$

すると, 特徴空間の内積はガウスカーネルとほぼ等しくなる. ガウスカーネルの幅は, ベクトル$\boldsymbol{x}$を再スケーリングすることで変更することができる.限界$l → ∞$では、ランダムな特徴写像を持つベイズ線形モデルはガウス過程に収束する.

##3. Benchmarking
COMBOは, 面心立方銅のΣ5[001](210)簡略化一致格子結晶粒界の剛体変換を最適化するために適用される.
並進のグリッドポイント17,983から, 最小の粒界エネルギーを生成するために最適な並進を見つける必要がある.$候補点x_iは, x,y,z$軸に沿った平行移動を表す3次元ベクトルである.$対応する値y_i$は, 経験的ポテンシャル法を使用した静的格子計算によって得られた粒界エネルギーである.

性能を測定するために,「成功」とは, 300回の実験の中でトップk点のうち少なくとも1点を見つけることと定義される. 最初の20個の実験は実行的に選択され, 残りの実験はCOMBOによって設計される.ランダム特徴量マップの特徴量lは2,000と5,000に設定されている.それぞれの場合, 異なる初期実験で30回のCOMBOを適用する.ｋの値を変えたときのCOMBOの成功確率を図２に示す.点線は,超幾何学的分布を用いて解析的に計算されたランダムデザインの成功確率を示している.ランダム設計と比較して, COMBOはトップポイントを見つけるのにはるかに効果的でした。例えば, k = 30の場合, ランダム設計の成功確率は40%であった. それがCOMBOによって83%(l = 2000)と90%(l = 5000)まで向上する.

図3は, scikit-learnとCOMBOの計算時間を比較したものである.scikit-learnの計算時間は2次関数的に伸びるのに対し, COMBOの計算時間は直線的に伸びている.非常に小さな問題ではscikit-learnの方が高速だが, 大規模な問題では予想通りCOMBOの方がはるかに効率が良いことがわかる.

##4. Discussion and Conclusion
COMBOにはユーザーの便宜のためにいくつかの機能がある.
第１に, COMBOには、ユーザーが独自のシミュレーターを使用するために簡単に変更できるインターフェース機能がある.
第2に, タイプIIの尤度を最大化することにより, カーネル幅やノイズ分散などの予測モデルのハイパーパラメーターをデータから自動的に学習できる.

学習データの数が増えると, ハイパーパラメーターをより正確に制御できる. COMBOは, ユーザーが指定した間隔でハイパーパラメーターを繰り返し更新するように設計されている.

COMBOは, より多くのモデルと機能を含むようにさらに開発され, 誰もがgithubの開発に参加するよう招待される.pythonパッケージとして, pymatgenなどの他の材料情報学パッケージと簡単に組み合わせることができる.