Umbrella samping/ガウス過程回帰による自由エネルギー曲線の予測

分子動力学シミュレーションにおけるUmbrella sampingは，対象とした自由度に対してバイアスポテンシャルを設定し，その平衡点の周辺を重点的に計算することで効率的に自由エネルギー曲線（平均力ポテンシャル）を計算する手法です。
反応経路の目星がある程度ついている状況において特に威力を発揮します。
Umbrella sampingの結果から自由エネルギー曲線を解析する手法としては，Weighted Histogram Analysis Method（WHAM）が実質デファクトスタンダードとなっていますが，2014年にガウス過程回帰を用いた解析方法（参考文献１）が提案されました。
本記事では，Umbrella samping/ガウス過程回帰による自由エネルギー曲線の予測手法を解説します。基本的には参考文献１を抜粋した内容となっておりますが，数式展開については参考文献１よりも詳しく記載しています。

ガウス過程回帰について

そもそもガウス過程回帰とは何ぞやと思われる方もいらっしゃいますかと思いますが，参考文献３によりますと以下のようになります。

ガウス過程回帰は，入力変数$${\mathbf{x}}$$から出力変数である実数値$${y}$$への関数$${y=f(\mathbf{x})}$$を推定するモデルの一つである．その特徴の一つはその非線形性であり，線形回帰ではうまくフィッテイングできない場合にも有効である．もう一つの重要な特徴はベイズ推定を用いる点である．推定される関数は一つの関数ではなく，関数の分布として得られるので，推定の不確実性を表現することができる．

参考文献３

ガウス過程回帰を解説することが目的では在りませんので，本記事ではこれ以上ガウス過程回帰そのものの説明はしません。
興味ある方は参考文献３等をご参照ください。

Umbrella Sampingについて

系内の配置に関する全自由度を$${\mathbf{q}}$$，ポテンシャルエネルギーを$${V(\mathbf{q})}$$，及び$${\mathbf{q}}$$で定義される興味ある自由度を$${\hat{s}(\mathbf{q})}$$とします。
配置の出現確率がボルツマン分布に従うとき，$${\hat{s}(\mathbf{q})=x}$$となる確率$${P(x)}$$は以下の式で表されます。

$$
\begin{align*}
P(x)&=\frac{1}{Z}\int{\rm d}\mathbf{q}\exp\left(-\beta V(\mathbf{q})\right)\delta(\hat{s}(\mathbf{q})-x)\\
\beta&:=\frac{1}{k_{\rm B}T}\\
Z&:=\int{\rm d}\mathbf{q}\exp\left(-\beta V(\mathbf{q})\right)
\end{align*}\tag{1}
$$

ここで，$${k_{\rm B}, T}$$はそれぞれボルツマン定数，絶対温度です。
Umbrella Sampingでは，ポテンシャルエネルギー以外に$${\hat{s}(\mathbf{q})}$$を変数とするバイアスポテンシャル$${w_{\theta}(\hat{s}(\mathbf{q}))}$$を課すことにより，特定領域を重点的にサンプリングすることを考えます。
ある状態間の反応経路を重点的にサンプリングするため，平衡点の異なるバイアスポテンシャルを複数用意し，それぞれのバイアスポテンシャルを課したMD計算を個別実行，もしくはレプリカ交換サンプリングすることが一般的です。ここでは複数あるバイアスポテンシャルのインデックスを$${\theta}$$とします。
また，バイアスポテンシャルの関数形としてはハーモニックポテンシャルが用いることが多いため，$${\theta}$$番目のバイアスポテンシャルは

$$
\begin{align*}
w_{\theta}(\hat{s}(\mathbf{q}))&=\frac{\kappa_{\theta}}{2}(\hat{s}(\mathbf{q})-x_{\theta})^2
\end{align*}\tag{2}
$$

のようになります。$${V(\mathbf{q})+w_{\theta}(\hat{s}(\mathbf{q}))}$$の系はumbrella windowと称されます。
$${w_{\theta}(\hat{s}(\mathbf{q}))}$$が課せられた系では$${\hat{s}(\mathbf{q})=x}$$となる確率$${P_{\theta}^b(x)}$$は以下の式で表されます。

$$
\begin{align*}
P_{\theta}^b(x)&=\frac{1}{Z_{\theta}^b}\int{\rm d}\mathbf{q}\exp\left(-\beta( V(\mathbf{q})+w_{\theta}(\hat{s}(\mathbf{q})))\right)\delta(\hat{s}(\mathbf{q})-x)\\
&=\frac{\exp\left(-\beta w_{\theta}(x)\right)}{Z_{\theta}^b}\int{\rm d}\mathbf{q}\exp\left(-\beta V(\mathbf{q})\right)\delta(\hat{s}(\mathbf{q})-x)\\
&=\exp\left(-\beta w_{\theta}(x)\right)\frac{Z}{Z_{\theta}^b}P(x)\\
Z_{\theta}^b&:=\int{\rm d}\mathbf{q}\exp\left(-\beta( V(\mathbf{q})+w_{\theta}(\hat{s}(\mathbf{q})))\right)
\end{align*}\tag{3}
$$

式(3)より，$${\theta}$$番目のumbrella windowのサンプリング結果から目的の自由エネルギー曲線$${f(x)}$$は

$$
\begin{align*}
f(x)&=-\frac{1}{\beta}\ln P(x)\\
&=-\frac{1}{\beta}\ln P_{\theta}^b(x)-w_{\theta}(x)+\frac{1}{\beta}\ln\left(\frac{Z}{Z_{\theta}^b}\right)\\
&=y_{\theta}(x)-f_{0}^{\theta}\\
y_{\theta}(x)&:=-\frac{1}{\beta}\ln P_{\theta}^b(x)-w_{\theta}(x)\\
f_{0}^{\theta}&=-\frac{1}{\beta}\ln\left(\frac{Z}{Z_{\theta}^b}\right)
\end{align*}\tag{4}
$$

$${\theta}$$番目のumbrella windowで$${N_{\theta}}$$回サンプリングした際，$${x}$$となる状態が$${n_{\theta}(x)}$$回カウントされたとすると，$${P_{\theta}^b(x)}$$は

$$
\begin{align*}
P_{\theta}^b(x)&\simeq\frac{n_{\theta}(x)}{N_{\theta}}
\end{align*}\tag{5}
$$

と推定されます。
異なるumbrella window間でサンプリングされた領域が重なる場合，辻褄が合うように$${f(x)}$$を決める必要があります。

$$
\begin{align*}
y_{\theta}(x)-f_{0}^{\theta} &=y_{\theta'}(x)-f_{0}^{\theta'}
\end{align*}\tag{6}
$$

自由エネルギー曲線の予測分布

Umbrella samplingによって各umbrella windowで$${\{x_i,y_{\theta}(x_i)\}, i=1,2,\cdots, N_{\theta}}$$が得られたとします。
全データを以下のようにまとめて表現することにします。

$$
\begin{align*}
\mathbf{x}&=\begin{bmatrix}\mathbf{x}_{1}^{\rm T}&\cdots&\mathbf{x}_{\theta}^{\rm T}&\cdots&\mathbf{x}_{w}^{\rm T}\end{bmatrix}^{\rm T}\\
\mathbf{x}_{\theta}^{\rm T}&=\begin{bmatrix}x_1&\cdots&x_{N_{\theta}}\end{bmatrix}^{\rm T}\\
\mathbf{y}&=\begin{bmatrix}\mathbf{y}_{1}^{\rm T}&\cdots&\mathbf{y}_{\theta}^{\rm T}&\cdots&\mathbf{y}_{w}^{\rm T}\end{bmatrix}^{\rm T}\\
\mathbf{y}_{\theta}^{\rm T}&=\begin{bmatrix}y_{\theta}(x_1)&\cdots&y_{\theta}(x_{N_{\theta}})\end{bmatrix}^{\rm T}\\
\end{align*}\tag{7}
$$

ここで，$${w}$$はumbrella window数を表します。
目的は，得られたデータ$${\mathbf{y}}$$から尤もらしい$${f(x)}$$を予測することですが，ガウス過程回帰では観測データ$${\mathbf{y}}$$の事後分布として$${f(x)}$$が満たす確率分布を考えます。
予測したいデータ点が$${M}$$個あるとし，

$$
\begin{align*}
\mathbf{x}^*&:=\begin{bmatrix}x_1^*&\cdots&x_M^*\end{bmatrix}^{\rm T}\\
\mathbf{f}(\mathbf{x}^*)&:=\begin{bmatrix}f(x_1^*)&\cdots&f(x_M^*)\end{bmatrix}^{\rm T}\\
\end{align*}\tag{8}
$$

と定義すると，$${p(\mathbf{f}(\mathbf{x}^*)|\mathbf{y})}$$が対象となる確率分布です。
ベイズの定理を利用すると，$${p(\mathbf{f}(\mathbf{x}^*)|\mathbf{y})}$$は

$$
\begin{align*}
p(\mathbf{f}(\mathbf{x}^*)|\mathbf{y})&=\frac{p(f(x^*), \mathbf{y})}{p(\mathbf{y})}\\
&=\frac{\int{\rm d}\mathbf{f}(\mathbf{x})p(\mathbf{f}(\mathbf{x}), \mathbf{f}(\mathbf{x}^*), \mathbf{y})}{p(\mathbf{y})}\\
&=\frac{\int{\rm d}\mathbf{f}(\mathbf{x})p(\mathbf{y}|\mathbf{f}(\mathbf{x}), f(x^*))p(\mathbf{f}(\mathbf{x}), \mathbf{f}(\mathbf{x}^*))}{p(\mathbf{y})}\\
&=\frac{\int{\rm d}\mathbf{f}(\mathbf{x})p(\mathbf{y}|\mathbf{f}(\mathbf{x}))p(\mathbf{f}(\mathbf{x}), \mathbf{f}(\mathbf{x}^*))}{p(\mathbf{y})}\\
&\propto \int{\rm d}\mathbf{f}(\mathbf{x})p(\mathbf{f}(\mathbf{x}), \mathbf{f}(\mathbf{x}^*))p(\mathbf{y}|\mathbf{f}(\mathbf{x}))
\end{align*}\tag{8}
$$

と式変形できます。
式(8)の積分の中には２つの確率分布がありますが，ガウス過程回帰ではその２つがガウス分布であると仮定します。
$${p(\mathbf{f}(\mathbf{x}))}$$は，

$$
\begin{align*}
p(\mathbf{f}(\mathbf{x}))&=\mathcal{N}(\mathbf{f}(\mathbf{x})|\mathbf{0},\mathbf{K}(\mathbf{x}))
\end{align*}\tag{9}
$$

と定義されます。
$${\mathbf{K}(\mathbf{x})}$$は各要素がカーネル関数で構成される行列になります。

$$
\begin{align*}
\mathbf{K}(\mathbf{x})&=\begin{bmatrix}\mathbf{k}(x_1,\mathbf{x})&\cdots&\mathbf{k}(x_N,\mathbf{x})\end{bmatrix}\\
\mathbf{k}(x_n, \mathbf{x})&=\begin{bmatrix}k(x_1,x_n)&\cdots& k(x_N,x_n)\end{bmatrix}^{\rm T}
\end{align*}\tag{10}
$$

一方，$${p(\mathbf{y}|\mathbf{f}(\mathbf{x}))}$$は，

$$
\begin{align*}
p(\mathbf{y}|\mathbf{f}(\mathbf{x}))&=\mathcal{N}(\mathbf{y}|\mathbf{f}(\mathbf{x}),\boldsymbol\Sigma_y)
\end{align*}\tag{11}
$$

と定義されます。
ガウス分布の性質を利用することにより，式(8)の積分を実際に計算すると

$$
\begin{align*}
p(\mathbf{f}(\mathbf{x}^*)|\mathbf{y})&=\mathcal{N}\left(\mathbf{f}(\mathbf{x}^*)\left|\boldsymbol\mu,\boldsymbol\Sigma\right.\right)\\
\boldsymbol\mu&=\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\left(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y\right)^{-1}\mathbf{y}\\
\boldsymbol\Sigma&=\mathbf{K}(\mathbf{x}^*)-\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\left(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y\right)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})
\end{align*}\tag{12}
$$

となります。計算手順の詳細はAppendix Aをご参照ください。
ただし，式(12)では式(4)の$${f_{0}^{\theta}}$$の存在が考慮されていないことに注意してください。
$${\mathbf{f}_0=\begin{bmatrix}f_{0}^{1}&\cdots&f_{0}^{w}\end{bmatrix}^{\rm T}}$$が一様分布に従うと仮定した場合，式(12)は以下のように拡張されます。

$$
\begin{align*}\boldsymbol\mu&=\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\widetilde{\mathbf{K}}(\mathbf{x})\mathbf{y}\\
\boldsymbol\Sigma&=\mathbf{K}(\mathbf{x}^*)-\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\widetilde{\mathbf{K}}(\mathbf{x})\mathbf{K}(\mathbf{x}^*,\mathbf{x})
\end{align*}\tag{13}
$$

ここで，$${\mathbf{H}}$$は$${\mathbf{y}}$$の$${i}$$成分が$${\theta}$$番目のumbrella windowに含まれているときに$${H_{\theta i}=1}$$，そうでない場合に $${H_{\theta i}=0}$$の要素を持つ行列です。
また，$${\widetilde{\mathbf{K}}(\mathbf{x})}$$は

$$
\begin{align*}
\widetilde{\mathbf{K}}(\mathbf{x})&:=(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y)^{-1}\left[\mathbf{I}-\mathbf{H}^{\rm T}\left\{\mathbf{H}(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y)^{-1}\mathbf{H}^{\rm T}\right\}^{-1}\mathbf{H}(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y)^{-1}\right]
\end{align*}
$$

です。計算手順の詳細はAppendix Bをご参照ください。
式(13)には複雑な行列計算が含まれますので，可能であればより簡単に計算できる定式化が望まれます。
各umbrella windowでリファレンス値となる自由エネルギー値を１点選び，リファレンス値との相対的な自由エネルギー差をベースに定式化することにより，式(13)と等価であるにも関わらずより簡単に計算できる式に帰着させることができます。証明についてはAppendix Cをご参照ください。

カーネル関数について

ガウス過程回帰のカーネル関数として様々な関数形がこれまで提案されてきており，理想的には適用する系の性質に併せて選択するべきです。
参考文献１では，少なくとも文献内で扱った系ではガウスカーネルで十分にパフォーマンスが得られたと結論付けられています。

$$
\begin{align*}k(x_i,x_j)&=\sigma_f^2\exp\left(-\frac{(x_i-x_j)^2}{2l^2}\right)
\end{align*}\tag{14}
$$

$${\sigma_f, l}$$はハイパーパラメータであり，系の性質に併せて適切な値に設定することが求められます。

∑_yの推定

$${p(\mathbf{f}(\mathbf{x}^*)|\mathbf{y})}$$の分散共分散行列である$${\boldsymbol\Sigma_y}$$をどう用意するかも解析に大きく影響しえます。
Umbrella samplingの結果から$${\boldsymbol\Sigma_y}$$を推定することにします。
Umbrella window $${\theta}$$で得られたサンプル数を$${N_{\theta}}$$，ビン$${\theta_i}$$がサンプリングされる確率を$${p_{\theta_i}}$$とします。
このとき，各ビンのサンプル数$${\{n_{\theta_i}\}}$$は多項分布に従うため，$${N_{\theta}}$$が全て独立試行で得られたと仮定すると，

$$
\begin{align*}
{\rm cov}(n_{\theta_i},n_{\theta_j})&=N_{\theta}(\delta_{ij}p_{\theta_i}-p_{\theta_i}p_{\theta_j})
\end{align*}\tag{15}
$$

となります。
しかしながら，現実的には実際のサンプリングデータ間には相関が生じてしまう場合が懸念されます。
そこで，実際には独立試行で得られたサンプリングデータ数が$${N_{\theta,{\rm eff}}=\sum_{\theta_i}n_{\theta_i,{\rm eff}}}$$であったとすると，

$$
\begin{align*}
{\rm cov}(n_{\theta_i},n_{\theta_j})&={\rm cov}\left(\left(\frac{N_{\theta}}{N_{\theta,{\rm eff}}}\right)n_{\theta_i,{\rm eff}},\left(\frac{N_{\theta}}{N_{\theta,{\rm eff}}}\right)n_{\theta_j,{\rm eff}}\right)\\
&=\left(\frac{N_{\theta}}{N_{\theta,{\rm eff}}}\right)^2{\rm cov}\left(n_{\theta_i,{\rm eff}},n_{\theta_j,{\rm eff}}\right)\\
&=\left(\frac{N_{\theta}}{N_{\theta,{\rm eff}}}\right)^2N_{\theta,{\rm eff}}(\delta_{ij}p_{\theta_i}-p_{\theta_i}p_{\theta_j})
\end{align*}\tag{16}
$$

となります。
$${y_{\theta_i}=-\beta^{-1}\ln(n_{\theta_i}/N_{\theta})-w_{\theta}(x_{\theta_{i}})}$$及び，1次の誤差伝搬を用いることにより，

$$
\begin{align*}
{\rm cov}(y_{\theta_i},y_{\theta_j})&\simeq\frac{\partial y_{\theta_i}}{\partial n_{\theta_i}}\frac{\partial y_{\theta_j}}{\partial n_{\theta_j}}{\rm cov}(n_{\theta_i},n_{\theta_j})\\
&=\frac{1}{\beta^2n_{\theta_i}n_{\theta_j}}\left(\frac{N_{\theta}}{N_{\theta,{\rm eff}}}\right)^2N_{\theta,{\rm eff}}(\delta_{ij}p_{\theta_i}-p_{\theta_i}p_{\theta_j})\\
&\simeq\frac{1}{\beta^2n_{\theta_i}n_{\theta_j}}\left(\frac{N_{\theta}}{N_{\theta,{\rm eff}}}\right)^2N_{\theta,{\rm eff}}\left(\delta_{ij}\frac{n_{\theta_i}}{N_{\theta}}-\frac{n_{\theta_i}}{N_{\theta}}\frac{n_{\theta_j}}{N_{\theta}}\right)\\
&=\frac{1}{\beta^2N_{\theta,{\rm eff}}}\left(\delta_{ij}\frac{N_{\theta}}{n_{\theta_i}}-1\right)\\
\end{align*}\tag{17}
$$

が得られます。
式(17)はumbrella samplingによるサンプリングデータから$${\boldsymbol\Sigma_y}$$を計算できる形となっています。

実行手順

以上の内容から実行する手順を以下に記します。

1. $${\sigma_f,l}$$を事前に用意

2. Umbrella samplingを実行し，サンプリングデータから$${\mathbf{y}(\mathbf{x})}$$と$${N_{\theta,{\rm eff}}}$$を解析

3. $${\boldsymbol\Sigma_y}$$を計算

4. $${\mathbf{K}(\mathbf{x}),\mathbf{K}(\mathbf{x}^*),\mathbf{K}(\mathbf{x}^*,\mathbf{x})}$$を計算

5. 式(C3)の$${\boldsymbol\mu=\left(\mathbf{B}_{\Delta}\mathbf{K}(\mathbf{x}^*,\mathbf{x})\right)^{\rm T}\left(\mathbf{B}_{\Delta}\mathbf{K}_y\mathbf{B}_{\Delta}^{\rm T}\right)^{-1}\mathbf{B}_{\Delta}\mathbf{y}}$$を計算し，$${\boldsymbol\mu}$$を$${\mathbf{f}(\mathbf{x}^*)}$$と見なす

6. 式(C3)の$${\boldsymbol\Sigma=\mathbf{K}(\mathbf{x}^*)-\mathbf{K}_{\Delta}(\mathbf{x}^*,\mathbf{x})^{\rm T}\mathbf{K}_{\Delta y}^{-1}\mathbf{K}_{\Delta}(\mathbf{x}^*,\mathbf{x})}$$を計算し，$${\boldsymbol\Sigma}$$の対角成分を$${\mathbf{f}(\mathbf{x}^*)}$$の推定誤差（の2乗）と見なす

実際に試した際の所感

手元にあったumbrella samplingのデータ（公開できないデータですみません。。。）を用いて挙動を調べてみたのですが，想像していた以上に$${\sigma_f,l}$$の値に対する結果への影響が大きく，慎重に決める必要があるように見受けられました。
（確かに論文が主張するように殆ど結果が変わらない範囲があることも確認できましたが，それらを調査なしに見出すことは困難だと思います）
また，解析結果は$${N_{\theta,{\rm eff}}}$$にもそれなりに依存するので，$${N_{\theta,{\rm eff}}}$$の値を慎重に決める必要があることも分かりました。
当初は，論文の内容をよく理解せずに結果の画だけを見て画期的な手法なのではと過剰に期待してしまっていたのですが，現在はむしろWHAMが如何に優れた手法であるかを思い知らされている感じがしています。
理論的に面白いとは今でも感じていますが，WHAMと比較して，

• 理論内に導入されている仮定が多い

• 結果に（それなりに大きく）影響するパラメータが複数ある

ことを考慮すると実践に適用するには相当に難しいですし，WHAMにそれほど不満がないため，それらの困難を乗り越えてでも使用するモチベーションがありません。
ガウス過程回帰による自由エネルギー曲線の予測論文が公開されて８年が経過していますが，いまだにWHAMがデファクトスタンダードとなっていることに納得した次第です。

参考文献

1. J. Chem. Theory Comput. 2014, 10, 9, 4079–4097

2. Christopher Bishop, Pattern Recognition and Machine Learning

3. システム/制御/情報, Vol. 62, No. 10, pp. 390-395, 2018

Appendix 

A. 式(12)の導出

$${\hat{\mathbf{x}}:=\begin{bmatrix}\mathbf{x}^{\rm T}&(\mathbf{x}^*)^{\rm T}\end{bmatrix}^{\rm T}}$$とおくと，

$$
\begin{align*}
p(\mathbf{f}(\mathbf{x}),\mathbf{f}(\mathbf{x}^*))&=p(\mathbf{f}(\hat{\mathbf{x}}))\\
&=\mathcal{N}(\mathbf{f}(\hat{\mathbf{x}})|\mathbf{0},\mathbf{K}(\hat{\mathbf{x}}))\\
&\propto \exp\left(-\frac{1}{2}\mathbf{f}(\hat{\mathbf{x}})^{\rm T}\mathbf{K}(\hat{\mathbf{x}})^{-1}\mathbf{f}(\hat{\mathbf{x}})\right)
\end{align*}\tag{A1}
$$

と表すことができます。
ここで，$${\mathbf{K}(\hat{\mathbf{x}})}$$は，

$$
\begin{align*}
\mathbf{K}(\hat{\mathbf{x}})&=\begin{bmatrix}\mathbf{K}(\mathbf{x})&\mathbf{K}(\mathbf{x}^*,\mathbf{x})\\ \mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}&\mathbf{K}(\mathbf{x}^*)\end{bmatrix}\\
\mathbf{K}(\mathbf{x}^*,\mathbf{x})&=\begin{bmatrix}\mathbf{k}(x_1^*,\mathbf{x}) &\cdots & \mathbf{k}(x_M^*,\mathbf{x})\end{bmatrix}
\end{align*}\tag{A2}
$$

です。
区分行列に関する逆行列の公式より，

$$
\begin{align*}
\mathbf{K}(\hat{\mathbf{x}})^{-1}&=\begin{bmatrix}\mathbf{M}&-\mathbf{M}\mathbf{K}(\mathbf{x}^*,\mathbf{x})\mathbf{K}(\mathbf{x}^*)^{-1}\\ -\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\mathbf{M}&\mathbf{K}(\mathbf{x}^*)^{-1}+\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\mathbf{M}\mathbf{K}(\mathbf{x}^*,\mathbf{x})\mathbf{K}(\mathbf{x}^*)^{-1}\end{bmatrix}\\
\mathbf{M}&:=\left(\mathbf{K}(\mathbf{x})-\mathbf{K}(\mathbf{x}^*,\mathbf{x})\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\right)^{-1}
\end{align*}
$$

となることを利用して$${p(\mathbf{f}(\mathbf{x}),\mathbf{f}(\mathbf{x}^*))}$$の指数を展開すると，

$$
\begin{align*}
&\mathbf{f}(\hat{\mathbf{x}})^{\rm T}\mathbf{K}(\hat{\mathbf{x}})^{-1}\mathbf{f}(\hat{\mathbf{x}})\\
&=\mathbf{f}(\mathbf{x})^{\rm T}\mathbf{M}\left(\mathbf{f}(\mathbf{x})-\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})\mathbf{f}(\mathbf{x}^*)\right)\\
&\ \ \ \ -\mathbf{f}(\mathbf{x}^*)^{\rm T}\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})\mathbf{M}\mathbf{f}(\mathbf{x})\\
&\ \ \ \ +\mathbf{f}(\mathbf{x}^*)^{\rm T}\left\{\mathbf{K}(\mathbf{x}^*)^{-1}+\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\mathbf{M}\mathbf{K}(\mathbf{x}^*,\mathbf{x})\mathbf{K}(\mathbf{x}^*)^{-1}\right\}\mathbf{f}(\mathbf{x}^*)\\
&=\mathbf{f}(\mathbf{x}^*)^{\rm T}\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{f}(\mathbf{x}^*)\\
&\ \ \ \ +\left\{\mathbf{f}(\mathbf{x})-\mathbf{K}(\mathbf{x}^*,\mathbf{x})\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{f}(\mathbf{x}^*)\right\}^{\rm T}\mathbf{M}\left\{\mathbf{f}(\mathbf{x})-\mathbf{K}(\mathbf{x}^*,\mathbf{x})\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{f}(\mathbf{x}^*)\right\}
\end{align*}\tag{A3}
$$

が得られます。
式(A3)の結果をから式(A1)は，

$$
\begin{align*}
p(\mathbf{f}(\mathbf{x}),\mathbf{f}(\mathbf{x}^*))&\propto \exp\left(-\frac{1}{2}\mathbf{f}(\mathbf{x}^*)^{\rm T}\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{f}(\mathbf{x}^*)\right)\\
&\ \ \ \ \ \ \ \ \times\mathcal{N}\left(\mathbf{f}(\mathbf{x})\left|\mathbf{K}(\mathbf{x}^*,\mathbf{x})\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{f}(\mathbf{x}^*), \mathbf{M}^{-1}\right.\right)
\end{align*}\tag{A4}
$$

となることが分かります。
式(A4)を式(1)に代入すると，

$$
\begin{align*}
p(\mathbf{f}(\mathbf{x}^*)|\mathbf{y})&\propto \exp\left(-\frac{1}{2}\mathbf{f}(\mathbf{x}^*)^{\rm T}\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{f}(\mathbf{x}^*)\right)\\
&\ \ \ \ \ \ \ \ \times\int{\rm d}\mathbf{f}(\mathbf{x})\mathcal{N}\left(\mathbf{f}(\mathbf{x})\left|\mathbf{K}(\mathbf{x}^*,\mathbf{x})\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{f}(\mathbf{x}^*), \mathbf{M}^{-1}\right.\right)\mathcal{N}(\mathbf{y}|\mathbf{f}(\mathbf{x}),\boldsymbol\Sigma_y)
\end{align*}
$$

となります。
参考文献２の式(2.115)を利用して積分を計算すると，

$$
\begin{align*}
&\int{\rm d}\mathbf{f}(\mathbf{x})\mathcal{N}\left(\mathbf{f}(\mathbf{x})\left|\mathbf{K}(\mathbf{x}^*,\mathbf{x})\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{f}(\mathbf{x}^*), \mathbf{M}^{-1}\right.\right)\mathcal{N}(\mathbf{y}|\mathbf{f}(\mathbf{x}),\boldsymbol\Sigma_y)\\
&=\mathcal{N}\left(\mathbf{y}\left|\mathbf{K}(\mathbf{x}^*,\mathbf{x})\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{f}(\mathbf{x}^*),\boldsymbol\Sigma_y+\mathbf{M}^{-1}\right.\right)\\
\end{align*}\tag{A5}
$$

となります。
式(A5)を用いて式(A4)の指数を式変形すると，

$$
\begin{align*}
&-\frac{1}{2}\mathbf{f}(\mathbf{x}^*)^{\rm T}\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{f}(\mathbf{x}^*)-\frac{1}{2}\left(\mathbf{y}-\mathbf{K}(\mathbf{x}^*,\mathbf{x})\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{f}(\mathbf{x}^*)\right)^{\rm T}\left(\boldsymbol\Sigma_y+\mathbf{M}^{-1}\right)^{-1}\left(\mathbf{y}-\mathbf{K}(\mathbf{x}^*,\mathbf{x})\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{f}(\mathbf{x}^*)\right)\\
&=-\frac{1}{2}\mathbf{f}(\mathbf{x}^*)^{\rm T}\left\{\mathbf{K}(\mathbf{x}^*)^{-1}+\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\left(\boldsymbol\Sigma_y+\mathbf{M}^{-1}\right)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})\mathbf{K}(\mathbf{x}^*)^{-1}\right\}\mathbf{f}(\mathbf{x}^*)+\mathbf{f}(\mathbf{x}^*)^{\rm T}\left\{\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\left(\boldsymbol\Sigma_y+\mathbf{M}^{-1}\right)^{-1}\mathbf{y}\right\}+{\rm (others)}\\
&=-\frac{1}{2}\mathbf{f}(\mathbf{x}^*)^{\rm T}\left\{\mathbf{K}(\mathbf{x}^*)^{-1}+\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\left(\boldsymbol\Sigma_y+\mathbf{K}(\mathbf{x})-\mathbf{K}(\mathbf{x}^*,\mathbf{x})\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\right)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})\mathbf{K}(\mathbf{x}^*)^{-1}\right\}\mathbf{f}(\mathbf{x}^*)+\mathbf{f}(\mathbf{x}^*)^{\rm T}\left\{\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\left(\boldsymbol\Sigma_y+\mathbf{K}(\mathbf{x})-\mathbf{K}(\mathbf{x}^*,\mathbf{x})\mathbf{K}(\mathbf{x}^*)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\right)^{-1}\mathbf{y}\right\}+{\rm (others)}\\
&=-\frac{1}{2}\mathbf{f}(\mathbf{x}^*)^{\rm T}\left\{\mathbf{K}(\mathbf{x}^*)-\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\left(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y\right)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})\right\}^{-1}\mathbf{f}(\mathbf{x}^*)+\mathbf{f}(\mathbf{x}^*)^{\rm T}\left\{\mathbf{K}(\mathbf{x}^*)-\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\left(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y\right)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})\right\}^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\left(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y\right)^{-1}\mathbf{y}+{\rm (others)}\\
\end{align*}
$$

となります。式変形の途中で参考文献２の式(C.5)，(C.6)を利用しました。
以上より，

$$
\begin{align*}
p(\mathbf{f}(\mathbf{x}^*)|\mathbf{y})&=\mathcal{N}\left(\mathbf{f}(\mathbf{x}^*)\left|\boldsymbol\mu,\boldsymbol\Sigma\right.\right)\\
\boldsymbol\mu&=\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\left(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y\right)^{-1}\mathbf{y}\\
\boldsymbol\Sigma&=\mathbf{K}(\mathbf{x}^*)-\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\left(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y\right)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})
\end{align*}\tag{A6}
$$

となることが分かります。

B. 式(13)の導出

$${p(\mathbf{y}|\mathbf{f}(\mathbf{x}))}$$を$${p(\mathbf{y}-\mathbf{H}^{\rm T}\mathbf{f}_0|\mathbf{f}(\mathbf{x}))}$$に拡張した際の$${p(\mathbf{f}(\mathbf{x}^*)|\mathbf{y})}$$を導出します。
$${p(\mathbf{f}_0)}$$が一様分布である場合，

$$
\begin{align*}
p(\mathbf{f}(\mathbf{x}^*)|\mathbf{y})&=\int{\rm d}\mathbf{f}_0p(\mathbf{f}(\mathbf{x}^*),\mathbf{f}_0|\mathbf{y})\\
&=\int{\rm d}\mathbf{f}_0p(\mathbf{f}(\mathbf{x}^*)|\mathbf{f}_0,\mathbf{y})p(\mathbf{f}_0|\mathbf{y})\\
&=\int{\rm d}\mathbf{f}_0p(\mathbf{f}(\mathbf{x}^*)|\mathbf{f}_0,\mathbf{y})\left\{\frac{p(\mathbf{y}|\mathbf{f}_0)p(\mathbf{f}_0)}{p(\mathbf{y})}\right\}\\
&=\frac{p(\mathbf{f}_0)}{p(\mathbf{y})}\int{\rm d}\mathbf{f}_0p(\mathbf{f}(\mathbf{x}^*)|\mathbf{f}_0,\mathbf{y})p(\mathbf{y}|\mathbf{f}_0)\\
&\propto\int{\rm d}\mathbf{f}_0p(\mathbf{f}(\mathbf{x}^*)|\mathbf{f}_0,\mathbf{y})p(\mathbf{y}|\mathbf{f}_0)
\end{align*}\tag{B1}
$$

が成立します。
$${p(\mathbf{f}(\mathbf{x}^*)|\mathbf{f}_0,\mathbf{y})}$$については，式(A6)の$${\mathbf{y}}$$を$${\mathbf{y}-\mathbf{H}^{\rm T}\mathbf{f}_0}$$に入れ替えればよいので，

$$
\begin{align*}
p(\mathbf{f}(\mathbf{x}^*)|\mathbf{f}_0,\mathbf{y})&=\mathcal{N}\left(\mathbf{f}(\mathbf{x}^*)\left|\boldsymbol\mu^{(1)},\boldsymbol\Sigma^{(1)}\right.\right)\\
\boldsymbol\mu^{(1)}&=\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\left(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y\right)^{-1}(\mathbf{y}-\mathbf{H}^{\rm T}\mathbf{f}_0)\\
\boldsymbol\Sigma^{(1)}&=\mathbf{K}(\mathbf{x}^*)-\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\left(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y\right)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})
\end{align*}\tag{B2}
$$

一方，$${p(\mathbf{y}|\mathbf{f}_0)}$$については，$${p(\mathbf{f}(\mathbf{x})|\mathbf{f}_0)=p(\mathbf{f}(\mathbf{x}))}$$であることを考慮すると，

$$
\begin{align*}
p(\mathbf{y}|\mathbf{f}_0)&=\int{\rm d}\mathbf{f}p(\mathbf{y},\mathbf{f}(\mathbf{x})|\mathbf{f}_0)\\
&=\int{\rm d}\mathbf{f}p(\mathbf{y}|\mathbf{f}(\mathbf{x}),\mathbf{f}_0)p(\mathbf{f}(\mathbf{x})|\mathbf{f}_0)\\
&=\int{\rm d}\mathbf{f}p(\mathbf{y}|\mathbf{f}(\mathbf{x}),\mathbf{f}_0)p(\mathbf{f}(\mathbf{x}))\\
&=\int{\rm d}\mathbf{f}\mathcal{N}\left(\mathbf{y}\left|\mathbf{f}(\mathbf{x})+\mathbf{H}^{\rm T}\mathbf{f}_0,\boldsymbol\Sigma_y\right.\right)\mathcal{N}\left(\mathbf{f}(\mathbf{x})\left|\mathbf{0},\mathbf{K}(\mathbf{x})\right.\right)\\
&=\mathcal{N}\left(\mathbf{y}\left|\mathbf{H}^{\rm T}\mathbf{f}_0,\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y\right.\right)\\
&\propto \mathcal{N}\left(\mathbf{f}_0\left|\boldsymbol\mu^{(2)},\boldsymbol\Sigma^{(2)}\right.\right)\\
&\ \ \ \ \ \ \ \ \times\exp\left(-\frac{1}{2}\mathbf{y}^{\rm T}(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y)^{-1}\mathbf{y}+\frac{1}{2}(\boldsymbol\mu^{(2)})^{\rm T}(\boldsymbol\Sigma^{(2)})^{-1}\boldsymbol\mu^{(2)}\right)\\
&\propto \mathcal{N}\left(\mathbf{f}_0\left|\boldsymbol\mu^{(2)},\boldsymbol\Sigma^{(2)}\right.\right)\\
\boldsymbol\mu^{(2)}&:=\left\{\mathbf{H}(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y)^{-1}\mathbf{H}^{\rm T}\right\}^{-1}\mathbf{H}(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y)^{-1}\mathbf{y}\\
\boldsymbol\Sigma^{(2)}&:=\left\{\mathbf{H}(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y)^{-1}\mathbf{H}^{\rm T}\right\}^{-1}
\end{align*}\tag{B3}
$$

となります。
式(B2)，式(B3)を式(B1)に代入すると，

$$
\begin{align*}
p(\mathbf{f}(\mathbf{x}^*)|\mathbf{y})&\propto\int{\rm d}\mathbf{f}_0p(\mathbf{f}(\mathbf{x}^*)|\mathbf{f}_0,\mathbf{y})p(\mathbf{y}|\mathbf{f}_0)\\
&=\int{\rm d}\mathbf{f}_0\mathcal{N}\left(\mathbf{f}(\mathbf{x}^*)\left|\boldsymbol\mu^{(1)},\boldsymbol\Sigma^{(1)}\right.\right)\mathcal{N}\left(\mathbf{f}_0\left|\boldsymbol\mu^{(2)},\boldsymbol\Sigma^{(2)}\right.\right)\\
&=\mathcal{N}\left(\mathbf{f}(\mathbf{x}^*)\left|\boldsymbol\mu^{(3)},\boldsymbol\Sigma^{(3)}\right.\right)\\
\end{align*}\tag{B4}
$$

と計算できます。
ここで，

$$
\begin{align*}\boldsymbol\mu^{(3)}&:=\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y)^{-1}\left(\mathbf{y}-\mathbf{H}^{\rm T}\boldsymbol\mu^{(2)}\right)\\
&=:\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\widetilde{\mathbf{K}}(\mathbf{x})\mathbf{y}\\
\widetilde{\mathbf{K}}(\mathbf{x})&:=(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y)^{-1}\left[\mathbf{I}-\mathbf{H}^{\rm T}\left\{\mathbf{H}(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y)^{-1}\mathbf{H}^{\rm T}\right\}^{-1}\mathbf{H}(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y)^{-1}\right]\\
\boldsymbol\Sigma^{(3)}&:=\boldsymbol\Sigma^{(1)}+\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y)^{-1}\mathbf{H}\boldsymbol\Sigma^{(2)}\mathbf{H}^{\rm T}(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})\\
&=\mathbf{K}(\mathbf{x}^*)-\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\left(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y\right)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})+\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y)^{-1}\mathbf{H}\left\{\mathbf{H}(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y)^{-1}\mathbf{H}^{\rm T}\right\}^{-1}\mathbf{H}^{\rm T}(\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y)^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})\\
&=\mathbf{K}(\mathbf{x}^*)-\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\widetilde{\mathbf{K}}(\mathbf{x})\mathbf{K}(\mathbf{x}^*,\mathbf{x})
\end{align*}
$$

です。

C. 自由エネルギー差をベースにした定式化

Appendix Bで示したように各umbrella windowで得られた平均力ポテンシャルの不定定数$${\mathbf{f}_0}$$の存在を露わに考慮した場合の計算式は複雑な行列計算を必要とします。
ここでは，各umbrella windowで得られた自由エネルギー値の一つをreference値として設定し，それ以外の値に対してreference値との差分を取ることにより，式(A7)の表式が式(B5)と等価であることを示します。
reference値との差分を取った自由エネルギー差のデータをベクトルとしてまとめたものを

$$
\begin{align*}
\Delta\mathbf{y}&=\begin{bmatrix}\Delta\mathbf{y}_1^{\rm T}&\cdots&\Delta\mathbf{y}_i^{\rm T}&\cdots&\Delta\mathbf{y}_w^{\rm T}\end{bmatrix}^{\rm T}
\end{align*}\tag{C1}
$$

と表すことにします。
ここで，$${\Delta\mathbf{y}_i}$$は$${i}$$番目のumbrella windowから得られた自由エネルギーの相対値であり，元のbin数$${b_i}$$から1を引いた$${b_i-1}$$次元のベクトルになります。また，$${w}$$はumbrella window数を表します。
各umbrella windowのbinを合わせた総数を$${N}$$すると，$${\Delta\mathbf{y}}$$は$${N-w}$$次元のベクトルになります。
$${\mathbf{y}_i}$$の最後の成分をreference値に選ぶことにすると，$${\Delta\mathbf{y}_i}$$と$${\mathbf{y}_i}$$には以下の関係が成立します。

$$
\begin{align*}
\Delta\mathbf{y}&=\mathbf{B}_{\Delta}\mathbf{y}\\
\mathbf{B}_{\Delta}&:=\begin{bmatrix}\mathbf{B}_{\Delta}^1&\mathbf{0}&\cdots&\mathbf{0}\\
\mathbf{0}&\mathbf{B}_{\Delta}^2&\mathbf{0}&\vdots\\
\vdots & \mathbf{0} &\ddots & \vdots\\
\mathbf{0} & \cdots & \cdots & \mathbf{B}_{\Delta}^w
\end{bmatrix}\\
\mathbf{B}_{\Delta}^i&=\begin{bmatrix}\mathbf{I}_{b_i-1}&-\mathbf{1}_{b_i-1}\end{bmatrix}
\end{align*}\tag{C2}
$$

ここで，$${\mathbf{I}_{b_i-1}}$$は$${(b_i-1)\times (b_i-1)}$$の単位行列，$${\mathbf{1}_{b_i-1}}$$は全成分が1の$${b_i-1}$$次元ベクトルです。$${\mathbf{B}_{\Delta}^i}$$は$${(b_i-1)\times b_i}$$行列となります。
$${\mathbf{K}_y:=\mathbf{K}(\mathbf{x})+\boldsymbol\Sigma_y}$$とし，$${\mathbf{B}_{\Delta}}$$を用いて式(A6)を以下のように書き換えます。

$$
\begin{align*}
\boldsymbol\mu&=\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\mathbf{K}_y^{-1}\mathbf{y}\\
&\rightarrow\left(\mathbf{B}_{\Delta}\mathbf{K}(\mathbf{x}^*,\mathbf{x})\right)^{\rm T}\left(\mathbf{B}_{\Delta}\mathbf{K}_y\mathbf{B}_{\Delta}^{\rm T}\right)^{-1}\mathbf{B}_{\Delta}\mathbf{y}\\
&=\mathbf{K}_{\Delta}(\mathbf{x}^*,\mathbf{x})^{\rm T}\mathbf{K}_{\Delta y}^{-1}\Delta\mathbf{y}\\
\boldsymbol\Sigma&=\mathbf{K}(\mathbf{x}^*)-\mathbf{K}(\mathbf{x}^*,\mathbf{x})^{\rm T}\mathbf{K}_y^{-1}\mathbf{K}(\mathbf{x}^*,\mathbf{x})\\
&\rightarrow\mathbf{K}(\mathbf{x}^*)-\mathbf{K}_{\Delta}(\mathbf{x}^*,\mathbf{x})^{\rm T}\mathbf{K}_{\Delta y}^{-1}\mathbf{K}_{\Delta}(\mathbf{x}^*,\mathbf{x})\\
\mathbf{K}_{\Delta}(\mathbf{x}^*,\mathbf{x})&:=\mathbf{B}_{\Delta}\mathbf{K}(\mathbf{x}^*,\mathbf{x})\\
\mathbf{K}_{\Delta y}&:=\mathbf{B}_{\Delta}\mathbf{K}_y\mathbf{B}_{\Delta}^{\rm T}
\end{align*}\tag{C3}
$$

次に

$$
\begin{align*}
\mathbf{B}_{\Delta}&=\mathbf{I}_{-w}\mathbf{B}\\
\mathbf{H}_w\mathbf{B}&=\mathbf{H}
\end{align*}\tag{C4}
$$

を満たす$${N\times N}$$の行列$${\mathbf{B}}$$を定義します。
ここで，$${\mathbf{I}_{-w}}$$と$${\mathbf{H}_w}$$は，$${N\times N}$$の単位行列$${\mathbf{I}_N}$$を構成する$${(N-w)\times N}$$，$${w\times N}$$の区分行列です。

$$
\begin{align*}
\mathbf{I}_N&=\begin{bmatrix}\mathbf{I}_{-w}\\\mathbf{H}_w\end{bmatrix}
\end{align*}\tag{C5}
$$

$${\mathbf{B}_{\Delta}}$$の成分ベクトルを$${\{\mathbf{b}_{\Delta i}\}, i=1,\cdots,N-w}$$，$${\mathbf{H}}$$の成分ベクトルを$${\{\mathbf{h}_i\}, i=1,\cdots,w}$$とおくと，

$$
\begin{align*}
\mathbf{b}_{\Delta i}^{\rm T}\mathbf{h}_j&=0\\
\mathbf{h}_i^{\rm T}\mathbf{h}_j&=b_i\delta_{ij}
\end{align*}\tag{C6}
$$

が成立します。この性質を利用すると，

$$
\begin{align*}
\mathbf{H}_w\mathbf{B}\mathbf{B}^{\rm T}&=\begin{bmatrix}\mathbf{h}_1^{\rm T}\mathbf{b}_{\Delta 1}&\cdots&\mathbf{h}_1^{\rm T}\mathbf{b}_{\Delta N-w}&\mathbf{h}_1^{\rm T}\mathbf{h}_1&\cdots&\mathbf{h}_1^{\rm T}\mathbf{h}_w\\
\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\
\mathbf{h}_w^{\rm T}\mathbf{b}_{\Delta 1}&\cdots&\mathbf{h}_w^{\rm T}\mathbf{b}_{\Delta N-w}&\mathbf{h}_w^{\rm T}\mathbf{h}_1&\cdots&\mathbf{h}_w^{\rm T}\mathbf{h}_w\end{bmatrix}\\
&=\begin{bmatrix}0&\cdots&0&b_1&0&\cdots&0\\
0&\cdots&0&0 & b_2 & 0&\vdots\\
0&\cdots&0&\vdots&0&\ddots&\vdots\\
0&\cdots&0&0&\cdots&0&b_w\end{bmatrix}\\
&=:{\rm diag}\left(\mathbf{h}^{\rm T}\mathbf{h}\right)\mathbf{H}_w\\
\therefore \mathbf{H}_w\left(\mathbf{B}^{\rm T}\right)^{-1}&={\rm diag}\left(\mathbf{h}^{\rm T}\mathbf{h}\right)^{-1}\mathbf{H}_w\mathbf{B}\\
&={\rm diag}\left(\mathbf{h}^{\rm T}\mathbf{h}\right)^{-1}\mathbf{H}
\end{align*}\tag{C7}
$$

が成立することが分かります。
また，参考文献２の式(2.91)より，逆行列をもつ任意の$${N\times N}$$行列$${A}$$に対して

$$
\begin{align*}
\left(\mathbf{I}_{-w}\mathbf{A}\mathbf{I}_{-w}^{\rm T}\right)^{-1}&=\mathbf{I}_{-w}\mathbf{A}^{-1}\mathbf{I}_{-w}^{\rm T}-\\
&\ \ \ \ \ \ \ \ \ \mathbf{I}_{-w}\mathbf{A}^{-1}\mathbf{H}_{w}^{\rm T}\left(\mathbf{H}_{w}\mathbf{A}^{-1}\mathbf{H}_{w}^{\rm T}\right)^{-1}\mathbf{H}_{w}\mathbf{A}^{-1}\mathbf{I}_{-w}^{\rm T}\\
\mathbf{I}_{-w}^{\rm T}\left(\mathbf{I}_{-w}\mathbf{A}\mathbf{I}_{-w}^{\rm T}\right)^{-1}\mathbf{I}_{-w}&=\mathbf{I}_{-w}^{\rm T}\mathbf{I}_{-w}\mathbf{A}^{-1}\mathbf{I}_{-w}^{\rm T}\mathbf{I}_{-w}-\\
&\ \ \ \ \ \ \ \ \ \mathbf{I}_{-w}^{\rm T}\mathbf{I}_{-w}\mathbf{A}^{-1}\mathbf{H}_{w}^{\rm T}\left(\mathbf{H}_{w}\mathbf{A}^{-1}\mathbf{H}_{w}^{\rm T}\right)^{-1}\mathbf{H}_{w}\mathbf{A}^{-1}\mathbf{I}_{-w}^{\rm T}\mathbf{I}_{-w}\\
&=\mathbf{A}^{-1}-\mathbf{A}^{-1}\mathbf{H}_{w}^{\rm T}\left(\mathbf{H}_{w}\mathbf{A}^{-1}\mathbf{H}_{w}^{\rm T}\right)^{-1}\mathbf{H}_{w}\mathbf{A}^{-1}\\
\end{align*}\tag{C8}
$$

が成立します。
式(C7)，(C8)を利用して$${\mathbf{B}_{\Delta}^{\rm T}\left(\mathbf{B}_{\Delta}\mathbf{K}_{y}\mathbf{B}_{\Delta}^{\rm T}\right)^{-1}\mathbf{B}_{\Delta}}$$を式変形すると，

$$
\begin{align*}
\mathbf{B}_{\Delta}^{\rm T}\left(\mathbf{B}_{\Delta}\mathbf{K}_{y}\mathbf{B}_{\Delta}^{\rm T}\right)^{-1}\mathbf{B}_{\Delta}&=\mathbf{B}^{\rm T}\mathbf{I}_{-w}^{\rm T}\left(\mathbf{I}_{-w}\mathbf{B}\mathbf{K}_{y}\mathbf{B}^{\rm T}\mathbf{I}_{-w}^{\rm T}\right)^{-1}\mathbf{I}_{-w}\mathbf{B}\\
&=\mathbf{B}^{\rm T}\left(\mathbf{B}\mathbf{K}_{y}\mathbf{B}^{\rm T}\right)^{-1}\mathbf{B}-\mathbf{B}^{\rm T}\left(\mathbf{B}\mathbf{K}_{y}\mathbf{B}^{\rm T}\right)^{-1}\mathbf{H}_{w}^{\rm T}\left(\mathbf{H}_{w}\left(\mathbf{B}\mathbf{K}_{y}\mathbf{B}^{\rm T}\right)^{-1}\mathbf{H}_{w}^{\rm T}\right)^{-1}\mathbf{H}_{w}\left(\mathbf{B}\mathbf{K}_{y}\mathbf{B}^{\rm T}\right)^{-1}\mathbf{B}\\
&=\mathbf{K}_{y}^{-1}-\mathbf{K}_{y}^{-1}\mathbf{B}^{-1}\mathbf{H}_{w}^{\rm T}\left(\mathbf{H}_{w}\left(\mathbf{B}^{\rm T}\right)^{-1}\mathbf{K}_{y}^{-1}\mathbf{B}^{-1}\mathbf{H}_{w}^{\rm T}\right)^{-1}\mathbf{H}_{w}\left(\mathbf{B}^{\rm T}\right)^{-1}\mathbf{K}_{y}^{-1}\\
&=\mathbf{K}_{y}^{-1}-\mathbf{K}_{y}^{-1}\mathbf{H}^{\rm T}{\rm diag}\left(\mathbf{h}^{\rm T}\mathbf{h}\right)^{-1}\left({\rm diag}\left(\mathbf{h}^{\rm T}\mathbf{h}\right)^{-1}\mathbf{H}\mathbf{K}_{y}^{-1}\mathbf{H}^{\rm T}{\rm diag}\left(\mathbf{h}^{\rm T}\mathbf{h}\right)^{-1}\right)^{-1}{\rm diag}\left(\mathbf{h}^{\rm T}\mathbf{h}\right)^{-1}\mathbf{H}\mathbf{K}_{y}^{-1}\\
&=\mathbf{K}_{y}^{-1}-\mathbf{K}_{y}^{-1}\mathbf{H}^{\rm T}\left(\mathbf{H}\mathbf{K}_{y}^{-1}\mathbf{H}^{\rm T}\right)^{-1}\mathbf{H}\mathbf{K}_{y}^{-1}\\
&=\widetilde{\mathbf{K}}(\mathbf{x})
\end{align*}
$$

となります。
これは，Appendix Bで得られた表式に他なりません。
以上より，自由エネルギー差による解析は自由エネルギーの不定値を露わに考慮した解析と等価であることが示されました。