\( \def\RR {\mathbb R } \def\BS(#1) {{\boldsymbol {#1} } } \def\ND(#1/#2/#3) {{\cal N}({#1}\, |\,{#2} ,{#3}) } \def\PD(#1/#2) {\frac{\partial #1}{\partial #2} } \def\CD(#1/#2) {p\,({#1}\,|\,{#2}) } \def\IP(#1/#2) {{\bf {#1}}\TP \bf {#2} } \def\EV(#1) {{\mathbb E}\,\lbrack{#1}\rbrack } \def\AG(#1/#2/#3) {\underset{\small {#2}} {{\rm arg\,#1}\ } {\displaystyle {#3}} } \def\TP{^\intercal} \)

8 ガウス過程回帰のハイパーパラメータ推定

　前回、ガウス過程回帰によって観測値から自由度が高く柔軟な関数の確率モデルが導かれ、複数のテスト入力に対する回帰値とバラツキを表す分散が推定できることを紹介した。
　ガウス過程回帰は、カーネル行列を共分散行列とするガウス分布が主役を務めていることも紹介した。そして、カーネル行列を構成するカーネル関数は、ハイパーパラメータをもつ関数でもあった。前回までは、固定したハイパーパラメータによって推定したが、より良い推定を行うためのハイパーパラメータをチューニングすることが求められる。このチューニングをハイパーパラメータ推定という。ここでは、動径基底関数を例に、最適なハイパーパラメータを推定する方法を紹介する。

8.1 尤度関数の最大化

　次式で表すように、３つのハイパーパラメータをまとめて\(\ \BS(\vartheta) =(\theta_1,\theta_2,\theta_3)\ \)とおくと、カーネル関数は\(\ \BS(\vartheta) \ \)に依存するので、明示的に 　 \[k({\bf x},{\bf x}’\,|\,\BS(\vartheta) ) = \theta_1 \exp \left( – \frac{|{\bf x}-{\bf x}’|^2}{\theta_2}\right) + \theta_3\,\delta({\bf x},{\bf x}’) \tag{8.1}\] と表すことができる。さらに、このカーネルから計算されるカーネル行列も\(\ \BS( \vartheta) \ \)に依存することになるので\(\ K_\BS(\vartheta) \ \)と表すこととする。 観測データ\(\ {\cal D}=(X,{\bf y})\ \)が与えられたときの出力ベクトルの尤度関数は、 \begin{eqnarray} p\,({\bf y}\,|\,X,\BS(\vartheta) )&=& {\cal N}({\bf y}\,|\,{\bf 0},K_\BS(\vartheta) ) \\ &=&\frac{1}{(2\pi)^{N/2}} \frac{1}{|K_\BS(\vartheta) |^{1/2}}\exp \left( -\frac12 {\bf y}\TP K_\BS(\vartheta) ^{-1}{\bf y} \right)\tag{8.2} \end{eqnarray} となる。最適なハイパーパラメータ推定は、この尤度関数を最大化する\(\ \BS(\vartheta) \ \)を求めることと等価であるので、式 (8.2) の対数を取った対数尤度関数\(\ L\ \)の最大化問題を考えることとする。　 \[ L=\log p\,({\bf y}\,|\,X,\BS(\vartheta) )= -\frac N2 \log(2\pi) -\frac 12 \log |K_\BS(\vartheta) | -\frac 12 {\bf y}\TP K_\BS(\vartheta) ^{-1}{\bf y} \tag{8.3} \] 　\(\ {\cal D} =\{\, (-3.5, 0.7), (-1.5, 1.8), (0, 1.6), (0.5, 2.3), (3.0, 1.0)\,\}\ \)の5個の観測値が与えられたとき、式 (8.1) のカーネル関数で定義されるガウス過程回帰モデルを考えることとする。ハイパーパラメータを固定した場合と、３つのハイパーパラメータの対数を等間隔\(\ 0.025\ \)のグリッド上に配置し、式 (8.3) を直接計算して\(\ L\ \)が最大になるグリッドを推定した場合を比較する。ハイパーパラメータの推定結果は、\(\ \log \BS(\vartheta) _{opt}= (0.400, 3.225, -1.950)\ \)で、\(\ L_{max}= -5.9189\ \)であった。図 8.1（左上）に\(\ \log \theta_1 = 0.4\ \)としたときの対数尤度を２次元平面上に等高線図としてプロットした。また、図 8.1（右下）には\(\ \BS(\vartheta) =(1,1,0.01)\ \)と固定した場合、図 8.1（左下）には\(\ \BS(\vartheta) _{opt}\ \)の場合に推定したガウス回帰過程を示した。対数尤度で\(\ 3.236\ \)の差（尤度ではおよそ25倍に相当）があることが示された。参考までに、固定したハイパーパラメータ近傍の対数尤度を図 8.1（右上）に示す。

　ハイパーパラメータ空間内のグリッドごとに式 (8.3) によって直接対数尤度を求めるのは計算負荷が大きく、この最大化問題を解く方法として様々な提案がなされているが、代表的な手法として

• MCMC法（Markov chain Monte Carlo method, マルコフ連鎖モンテカルロ法）
• SCG法やL-BFGS法などの勾配法（gradient method）

が挙げられる。ここでは、勾配法についてその概要を紹介する。なお、最適なハイパーパラメータの推定のために、MATLABやPythonのプラットフォームで利用できるライブラリが用意されているので改めて紹介の機会を設けることにする。

8.2 勾配法

　勾配法は、図 8.1(右上）の青色+を出発点とし、最良推定値である図 8.1（左上）の赤色+である頂上へ辿り着く過程に例えることができる。定性的で直感的ではあるが、出発点における勾配の方向へ進み、峠に到達したら改めて勾配を評価することを繰り返すことで頂上に到達することができる。
　そこで、式 (8.3) で定義する対数尤度の勾配をハイパーパラメータによる微分で評価することを考える。式 (8.3) の係数と定数項を省略した次式を微分する。 　 \[ {\cal L}= 2L=-\log |K_\BS(\vartheta) | -{\bf y}\TP K_\BS(\vartheta) ^{-1}{\bf y} \tag{8.4} \] 　式 (8.4) に以下の公式を適用することで、\(\ {\cal L}\ \)の\(\ \BS(\vartheta) \ \)の成分\(\ \theta \in \BS(\vartheta) \ \)による微分は、

\[\PD({\cal L}/\theta) = -{\rm tr} \left( K_{\BS(\vartheta) }^{-1} \PD(K_{\BS(\vartheta) }/\theta) \right)+(K_{\BS(\vartheta) }^{-1}{\bf y} )\TP \PD(K_{\BS(\vartheta) }/\theta) (K_{\BS(\vartheta) }^{-1} {\bf y}) \tag{8.5}\] と表すことができる。式 (8.5) は任意のカーネル関数について成立し、カーネル関数のハイパーパラメータの数にも制約がない重要な表現である。
　式 (8.1) をカーネル関数とするカーネル行列\(\ K_{\BS(\vartheta) }\ \)の微分\(\ {\partial K_{\BS(\vartheta) }}/\partial{\theta} \ \)を式 (8.5) に適用するとき、\(\ \theta_1 \gt 0\,,\theta_2 \gt0\,,\theta_3\gt 0\ \)であるので、 \[\theta_1 = \exp(\tau)\,,\, \theta_2=\exp(\sigma)\,,\,\theta_3=\exp(\eta)\ \tag{8.6}\] と変換すると、カーネル関数は \[k({\bf x},{\bf x}’\,|\,\BS(\vartheta) ) = {\rm e}^\tau \exp \left( – \frac{|{\bf x}-{\bf x}’|^2}{{\rm e}^\sigma} \right) + {\rm e}^\eta\,\delta({\bf x},{\bf x}’) \tag{8.7}\] さらに、カーネル関数の微分は \begin{equation} \left\{ \begin{array}{l} {\displaystyle \PD(k({\bf x},{\bf x}’\,|\,\BS(\vartheta) )/\tau) } &=& {\rm e}^\tau \exp \left( -{\displaystyle \frac{|{\bf x}-{\bf x}’|^2}{{\rm e}^\sigma} }\right) = k({\bf x},{\bf x}’\,|\,\BS(\vartheta) ) – {\rm e}^\eta\,\delta({\bf x},{\bf x}’) \\ {\displaystyle \PD(k({\bf x},{\bf x}’\,|\,\BS(\vartheta) )/\sigma) } &=& {\rm e}^\tau \exp \left( -{\displaystyle \frac{|{\bf x}-{\bf x}’|^2}{{\rm e}^\sigma} }\right) \cdot {\displaystyle\PD(/\sigma) }\left( -{\displaystyle \frac{|{\bf x}-{\bf x}’|^2}{{\rm e}^\sigma} }\right) \\ &=&\left( k({\bf x},{\bf x}’)-{\rm e}^\eta\,\delta({\bf x},{\bf x}’) \right) \cdot {\rm e}^{-\sigma}|{\bf x}-{\bf x}’|^2\\ {\displaystyle \PD(k({\bf x},{\bf x}’\,|\,\BS(\vartheta) )/\eta) } &=& {\rm e}^\eta\,\delta({\bf x},{\bf x}’)\tag{8.8} \end{array} \right. \end{equation} となる。１次微分で求めた式 (8.5) による勾配は最適解ではなく、２次微分を含めたより効率的なアルゴリズムである共役勾配法（CG法、conjugate gradient method）が提案されている。

8.3 ガウス過程回帰のまとめ

　ハイパーパラメータ推定については概要のみの紹介となったが、カーネル関数を異なるカーネルの線形結合として定義することもできる。たとえば、線形カーネルとガウスカーネルによって、長期間にわたる気温変動をモデル化することもできる。 \[ k({\bf x},{\bf x}’\,|\ \BS(\vartheta) )= \theta_1+ \theta_2\,{\bf x}\TP {\bf x}’+\theta_3\exp(-\frac{{|\bf x}-{\bf x}’|^2}{\theta_4})+\theta_5 \,\delta({\bf x},{\bf x}’)\] 　これまで線形回帰モデルから出発してガウス過程回帰を紹介してきたが、ガウス過程回帰は「観測値集合から柔軟な回帰モデルの関数\(\ y=f({\bf x})\ \)を推定する確率モデル」と位置付けることができる。
　モノづくりのライフサイクルにおいて、設計プロセスは「性能・品質・コスト」を決定する重要なプロセスである。図 8.2 に示すようにCAEツールと連携することで、部品・アセンブリ・システムそれぞれの階層において試作モデルを作成することなく、「性能・品質・コスト」のバラツキを予測することができる強力なツールを入手したことになる。
　ガウス過程回帰を実装した代理モデル（Emulator）をより効果的に運用するためには、以下の課題が残されている。

• 設計空間が３次元を超える場合のバラツキを可視化するツール
• 効率的な設計条件を探索する実験計画法（design of experiments, DOE）

　前者については名案はなく、頭の中で想像するか、主成分分析（principal component analysis, PCA）など次元削減手法の適用が現実的な対応である。
　次回は後者について紹介することとする。