チュートリアル09:ガウス関数を用いた自由エネルギー面のフィッティング

任意次元の自由エネルギー面の離散的な計算値を、なめらかな解析関数でフィッティングする。

計算の基礎

任意次元の離散データをガウス関数の重ね合わせでフィッティングする。

$f(x^{1},...,x^{n})=\sum _{j=1}^{N_{\text{Gauss}}}a(j)\exp \left[-{\frac {1}{2}}\sum _{l=1}^{n}\left({\frac {x^{l}-\mu (l,j)}{\sigma (l,j)}}\right)^{2}\right]$

• n=1次元 のとき、

$f(x)=\sum _{j=1}^{N_{\text{Gauss}}}a(j)\exp \left[-{\frac {1}{2}}\left({\frac {x-\mu (j)}{\sigma (j)}}\right)^{2}\right]$

• n=2次元 のとき、

$f(x,y)=\sum _{j=1}^{N_{\text{Gauss}}}a(j)\exp \left[-{\frac {1}{2}}\left({\frac {x-\mu (1,j)}{\sigma (1,j)}}\right)^{2}-{\frac {1}{2}}\left({\frac {y-\mu (2,j)}{\sigma (2,j)}}\right)^{2}\right]$

フィッティングされる離散データのセットを$\{f_{i}^{\text{Data}}\}$ (i=1, 2, ..., N) とするとき、フィッティングは残差の2乗

$L={\frac {1}{N}}\sum _{i=1}^{N}w_{i}\left|f(x_{i}^{1},...,x_{i}^{n})-f_{i}^{\text{Data}}\right|^{2}$

を最小化するように決められる。Nは離散データ点の数、$\{w_{i}\}$ (i=1, 2, ..., N)は重み因子である。

$N_{\text{Gauss}}$を与えたとき、最適化されるパラメータは $a(j),\mu (l,j),\sigma (l,j)$ $(j=1,2,...,N_{\text{Gauss}};l=1,2,...,n)$である。これらを$\mathbf {q} =\{q_{j}\}$ $(j=1,2,...,N_{\text{Fit}})$とまとめて表す。このとき、残差Lを最小化する（勾配${\partial L}/{\partial q_{j}}=0$にする）パラメータを求めるNewton-Raphson法の方程式は

${\begin{cases}\sum _{k=1}^{N_{\text{Fit}}}({\frac {\partial ^{2}L}{\partial q_{1}\partial q_{k}}})^{0}(q_{k}-q_{k}^{0})&+({\frac {\partial L}{\partial q_{1}}})^{0}&=0\\\vdots \\\sum _{k=1}^{N_{\text{Fit}}}({\frac {\partial ^{2}L}{\partial q_{N_{\text{Fit}}}\partial q_{k}}})^{0}(q_{k}-q_{k}^{0})&+({\frac {\partial L}{\partial q_{N_{\text{Fit}}}}})^{0}&=0\end{cases}}$

$\leftrightarrow {\begin{bmatrix}({\frac {\partial ^{2}L}{\partial q_{1}\partial q_{1}}})^{0}&\cdots &({\frac {\partial ^{2}L}{\partial q_{1}\partial q_{N_{\text{Fit}}}}})^{0}\\\vdots &\ddots &\vdots \\({\frac {\partial ^{2}L}{\partial q_{N_{\text{Fit}}}\partial q_{1}}})^{0}&\cdots &({\frac {\partial ^{2}L}{\partial q_{N_{\text{Fit}}}\partial q_{N_{\text{Fit}}}}})^{0}\end{bmatrix}}{\begin{bmatrix}q_{1}-q_{1}^{0}\\\vdots \\q_{N_{\text{Fit}}}-q_{N_{\text{Fit}}}^{0}\end{bmatrix}}+{\begin{bmatrix}({\frac {\partial L}{\partial q_{1}}})^{0}\\\vdots \\({\frac {\partial L}{\partial q_{N_{\text{Fit}}}}})^{0}\end{bmatrix}}={\begin{bmatrix}0\\\vdots \\0\end{bmatrix}}$

$\leftrightarrow \mathbf {A} (\mathbf {q} -\mathbf {q} ^{0})+\mathbf {b} =\mathbf {0}$

$\mathbf {A}$は正定値の対称行列である。パラメータの初期値$\mathbf {q} ^{0}$を代入して上式を解き、得られた$\mathbf {q}$を新しい$\mathbf {q} ^{0}$とする、という手続きを繰り返してパラメータを収束させる。ガウス関数の個数$N_{\text{Gauss}}$は、入力で与えられたフィッティングの許容誤差を満たすように決める。

なお、実際には残差Lのかわりにpenalty functionを加算した$L^{\text{tot}}$を用いている。

$L^{\text{tot}}=L+L^{\text{penalty}}=L+\lambda _{1}L_{1}^{\text{p}}+\lambda _{2}L_{2}^{\text{p}}+\lambda _{3}L_{3}^{\text{p}}$

ここで、

• $L_{1}^{\text{p}}$ : 極端に小さなパラメータ$\sigma$を禁止する。
• $L_{2}^{\text{p}}$ : 離散データの存在する座標領域から極端に離れた$\mu$を禁止する。
• $L_{3}^{\text{p}}$ : ガウス関数同士の極端に大きな重なりを禁止する。

$\lambda _{1},\lambda _{2},\lambda _{3}(>0)$はペナルティの強さを調整するための係数である。

gauss_fit.exeの内部

ガウス関数のパラメータに加え、使用するガウス関数の数$N_{\text{Gauss}}$も最適化しなければならない。本プログラムでは以下の手順を繰り返してフィッティングを遂行する。

1. 使用するガウス関数をひとつ増やし、初期パラメータを与える。

2. いま加えたガウス関数のパラメータのみを変数とし、ほかは全て固定してフィッティングする。

3. いま加えたガウス関数とひとつ古いガウス関数の2つだけを動かしてフィッティングする。

4. 同様に動かすガウス関数の数をひとつずつ増やしてゆき、最後には全てのパラメータをフィッティングする。

5. 終了判定。ユーザー指定の許容誤差を満たすか、あるいは残差Lがほとんど変化しなくなったら終了。

6. 1. に戻り、結果が確定するまで同じ手順を繰り返す。

終了判定について、まず許容誤差を満たすかどうかは次式により判定される。$L^{\text{init}}$はガウス関数が0個のときの残差、$\tau ^{\text{finish}}$はユーザー指定のthresholdである。

$L<L^{\text{init}}\times \tau ^{\text{finish}}$

また残差不変の判定は次式による。$L^{\text{old}}$はガウス関数が1つ少ないときの残差、$\tau ^{\text{unchange}}$はユーザー指定のthresholdである。

$|L-L^{\text{old}}|<L^{\text{old}}\times \tau ^{\text{unchange}}$

個々のフィッティングの収束判定も同様である。$L^{\text{oldstep}}$は繰り返し計算の1 step前の残差、$\tau ^{\text{iter}}$はユーザー指定のthresholdである。

$|L-L^{\text{oldstep}}|<L^{\text{oldstep}}\times \tau ^{\text{iter}}$

入力された離散データは、フィッティング前に以下の規則で変換される。当然、最終出力はもとのスケールに依拠する。

• 座標$\{x_{i}\}(i=1,2,...,n)$: すべての次元で最小値-1.0、最大値1.0となるように平行移動・スケーリング。
• 数値$y$: 絶対値の最大値が1.0となるようにスケーリング。

Newton-Raphson法では、初期値によってはパラメータが収束しない危険がある。これを回避するため、全パラメータ・全stepの変位の符号を必ず収束に向かうよう修正し、かつ絶対変位に上限を設定している。前述のスケーリングにより全パラメータに同じ上限値を適用可能。ユーザーは変換後のスケールで上限値を指定する。

プログラム gauss_fit の使用（1次元）

gauss_fit.exeの入力ファイルはネームリスト&gauss_fitを含む。最低限、入出力ファイル名を指定すれば適当にフィッティングしてくれる。

> cat input_1.nml 
&gauss_fit
  input_file_name = "data_1.txt"
  param_file_name = "param_1.output"
  bias_file_name = "bias_1.output"
  plot_file_name = "plot_1.output"
  log_file_name = "log_1"
  spcoord_ID = ""
  ngauss_max = 50
  niter_max = 100
  ncolumn_data = 1
  ncolumn_weight = 0
  omp_number = 1
  plot_result = .TRUE.
/

フィッティングされる離散データのセットはinput_file_nameに指定する。 ただし、ファイルの一番上の行に次の#plot文を追加する必要がある。

#plot a b c
  a - 座標の次元
  b - 座標上のデータの数
  c - 座標点の数（＝#plot文と空行を除いた行数）
ex) #plot 2 1 40000

input_file_nameには重み因子wを指定可能である。例えばファイルの各行が次の形式で書かれていたとする。

座標x_1　座標x_2　...　座標x_n　数値y_1　数値y_2　数値y_3

上の入力ネームリスト（およびデフォルトの設定）ではy_1を離散データの値として読み込み、重み因子として全離散点に1.0を適用する。これをy_2を入力データ、y_3を重み因子として読み込みたいときは、ネームリストの該当項目を次のように変更すればよい。

ncolumn_data = 2
ncolumn_weight = 3

さて、上のinput_1.nmlを入力としてgauss_fit.exeを起動する。チュートリアル09のディレクトリに入り、以下のコマンドを実行する。

> ../../../gauss_fit.exe input_1.nml

フィッティング終了後、自動的にグラフが表示される。不要ならばネームリストでplot_result=.FALSE.を指定する。

凡例	説明
data	入力離散データ
gauss	フィッティング結果
diff	残差（入力 - 結果）

plot_file_nameにはグラフと同じデータが出力される。左から座標x, data, gauss, diff。

> head -n 5 plot_1.output 
#plot          1         3       128
   -0.6350000000E-09    0.4216994350E-07    0.8711169105E-04   -0.8706952111E-04
   -0.6250000000E-09   -0.1550848046E-06    0.2878606196E-04   -0.2894114676E-04
   -0.6150000000E-09   -0.2710527917E-06   -0.1433599580E-04    0.1406494301E-04
   -0.6050000000E-09   -0.6318615346E-06   -0.4313367571E-04    0.4250181418E-04

log_file_nameにはフィッティングの経過が逐次出力される。ファイル名を省略した場合は標準出力される。これを参照すると、今回は14個目のガウス関数を追加した時点で残差Lが収束したためフィッティングを終了したことが分かる。パラメータは各行a, mu, sigの順に並んでいる。

> tail -n 25 log_1

# === Finished ===
# Residual does not change.
# Linit        =     0.1942166169E-02
# Lfin         =     0.3966710300E-06
# Lfin / Linit =     0.2042415506E-03

# Final parameters:
0.2749113636E-001  0.3413507899E-010  0.3669218306E-010
-0.1048779438E+001 -0.1608140032E-009  0.5772141359E-010
-0.4126744080E-002  0.5758140452E-009  0.3253049596E-010
0.1083004078E-001 -0.3860694684E-010  0.5897058793E-011
0.7875520328E-001 -0.1244001532E-009  0.2280624675E-010
-0.2463903832E+000  0.5059173022E-009  0.6521035761E-009
-0.9371230311E-002  0.5332822313E-010  0.1162332110E-010
0.8092789523E+000 -0.1676087679E-009  0.5315597783E-010
0.2251484120E+000 -0.9264616861E-010  0.9706211079E-010
0.2890017656E-001 -0.3187295130E-009  0.1998884124E-009
-0.2726601873E-002  0.1584640210E-009  0.2369364899E-010
0.1578408916E+000  0.5677966197E-009  0.4493305049E-009
0.7063114353E-001 -0.2606451312E-010  0.7996983269E-009
-0.1391012356E-001 -0.1492487262E-010  0.1173747077E-008

##### GAUSS_FIT END #####

param_file_nameにはフィッティングされたガウスパラメータが出力される。例によって左からa, mu, sigである。

> head -n 5 param_1.output 
  0.2749113636E-001  0.3413507899E-010  0.3669218306E-010
 -0.1048779438E+001 -0.1608140032E-009  0.5772141359E-010
 -0.4126744080E-002  0.5758140452E-009  0.3253049596E-010
  0.1083004078E-001 -0.3860694684E-010  0.5897058793E-011
  0.7875520328E-001 -0.1244001532E-009  0.2280624675E-010

bias_file_nameには同じパラメータが、FreeFlex入力ファイルのセクション&biasdataの形式で出力される。ただしパラメータaはkcal/molからJへと単位換算され、かつ符号が反転している。自由エネルギーを相殺するバイアスU0としての利用を想定。spcoord番号の部分はネームリストのspcoord_IDで指定できる。

> cat bias_1.output
&biasdata
GAUSS_1D   1 1   3 -0.1910000197E-21  0.3413507899E-10  0.3669218306E-10  0
GAUSS_1D   1 1   3  0.7286599243E-20 -0.1608140032E-09  0.5772141359E-10  0
GAUSS_1D   1 1   3  0.2867135757E-22  0.5758140452E-09  0.3253049596E-10  0
GAUSS_1D   1 1   3 -0.7524381591E-22 -0.3860694684E-10  0.5897058793E-11  0
GAUSS_1D   1 1   3 -0.5471671008E-21 -0.1244001532E-09  0.2280624675E-10  0
...
&end

プログラム gauss_fit の使用（2次元）

1次元の場合と同様である。

> cat input_2.nml 
&gauss_fit
  input_file_name = "data_2.txt"
  param_file_name = "param_2.output"
  bias_file_name = "bias_2.output"
  plot_file_name = "plot_2.output"
  log_file_name = "log_2"
  spcoord_ID = ""
  ngauss_max = 50
  niter_max = 100
  ncolumn_data = 1
  ncolumn_weight = 0
  omp_number = 1
  plot_result = .TRUE.
/

> ../../../gauss_fit.exe input_2.nml

> head -n 5 plot_2.output 
#plot          2         3       400
    0.0000000000E+00    0.0000000000E+00    0.1043823086E+01    0.1042287943E+01    0.1535142352E-02
    0.0000000000E+00    0.1000000000E-09    0.1057056936E+01    0.1055504872E+01    0.1552063903E-02
    0.0000000000E+00    0.2000000000E-09    0.1072252054E+01    0.1070630300E+01    0.1621753392E-02
    0.0000000000E+00    0.3000000000E-09    0.1088987338E+01    0.1087241260E+01    0.1746077385E-02

ログによれば、今回は3個目のガウス関数を追加した時点でフィッティングの目標精度を達成し終了した。2次元のためパラメータはa, mu_1, mu_2, sig_1, sig_2の順に並んでいる。

> tail -n 14 log_2

# === Finished ===
# Success.
# Linit        =     0.6647634059E+01
# Lfin         =     0.3145653153E-05
# Lfin / Linit =     0.4731989044E-06

# Final parameters:
0.4023808365E+001  0.9997095094E-009  0.1000379804E-008  0.3993443851E-009  0.6029095616E-009
0.8948314372E+000 -0.2419818741E-009  0.7666091501E-011  0.1106353303E-008  0.2450188450E-007
0.8897099946E+000  0.2177034840E-008  0.1085057848E-010  0.1097390843E-008  0.2133223238E-007

##### GAUSS_FIT END #####

> head param_2.output 
  0.4023808365E+001  0.9997095094E-009  0.1000379804E-008  0.3993443851E-009  0.6029095616E-009
  0.8948314372E+000 -0.2419818741E-009  0.7666091501E-011  0.1106353303E-008  0.2450188450E-007
  0.8897099946E+000  0.2177034840E-008  0.1085057848E-010  0.1097390843E-008  0.2133223238E-007

> cat bias_2.output 
&biasdata
GAUSS_2D   2 1 2   5 -0.2795619166E-19  0.9997095094E-09  0.1000379804E-08  0.3993443851E-09  0.6029095616E-09  0
GAUSS_2D   2 1 2   5 -0.6217015547E-20 -0.2419818741E-09  0.7666091501E-11  0.1106353303E-08  0.2450188450E-07  0
GAUSS_2D   2 1 2   5 -0.6181433328E-20  0.2177034840E-08  0.1085057848E-10  0.1097390843E-08  0.2133223238E-07  0
&end

ネームリスト変数一覧

こちらを参照。