Fortranで非線形最小二乗法(Gauss-Newton法)

こんにちは。Qiita初投稿です。
FortranとLapackを使って非線形最小二乗法を作ってみました。コード全文はGitHubに載せました。文末にURLを貼っております。

#記事を読んでいただく上で
この記事は著者が独学で分かった気になっていることを書いています。説明や解釈などが間違っていることがありますので、そこはそうじゃないだろ、というところがあれば指摘していただけると著者の勉強になります。よろしくお願いします。

#非線形最小二乗法
最小二乗法とは、測定したデータに合いそうな関数を当てはめる手法で、モデル関数(元になる形の近い関数)とデータの差の絶対値の二乗が最小になるように当てはめます。
モデル関数を$f(x)$、各点を$(x_i, y_i)$とすると、$f(x)$のパラメータを変えつつ
$$S = \sum(y_i – f(x_i))^2$$
を計算し、$S$が最小になるパラメータを探せばいいわけです。このときのパラメータとは、例えば$f(x) = ax + b$ならば、$(a, b)$のことです。
モデル関数が線形($y = ax + b$のような)であれば、連立一次方程式を解くことでパラメータが求まりますが、モデル関数が非線形のときは反復計算によって求めなければなりません。モデル関数が非線形のときの最小二乗法を非線形最小二乗法と呼びます。

#ガウス・ニュートン法
非線形最小二乗法の解法はいくつかあり、最急降下法やニュートン法、Marquardt
(たぶんマルカールって読むんだと思う)法、共役勾配法などがあります。ガウス・ニュートン法もその一つです。
ガウス・ニュートン法はニュートン法の派生法と言えます。ニュートン法は方程式を反復計算で解く方法で、微分可能であれば非線形な方程式でも解くことができるので、非線形最小二乗法だけでなくいろいろな問題に使えます。ニュートン法で非線形最小二乗法を解くには、関数の二回微分を含むヘッセ行列(ヘッシアン)が必要になります。二回微分を求めるのを回避する方法がガウス・ニュートン法です。
ガウス・ニュートン法では、ヘッシアン$H$をヤコビ行列(ヤコビアン)$J$の組み合わせで近似します。
$$H \approx J^TJ $$
ヤコビアンは一回微分までしか含まないため、計算がより簡単になります。

#式
ある程度解に近い初期値をパラメータ$\beta$として

$$r = y_i – f(x_i, \beta)$$
$$S = \sum (r^2)$$
が最小になるように
$$J^TJ \delta x = -J^T r $$
を解き、その解からパラメータを更新します。上記の式は線形なので、QR分解などで解けます。今回はLapackに解いてもらいます。
$S$が最小かどうか(収束したか)は、$S$がほぼ変化しなくなる($S$の変化が規定値以下になる)
まで、で判断することが多いようです。
今回は収束具合が見たいので、10回繰り返して$S$の変化を見ます。

#プログラム
例としてWikipediaのガウス・ニュートン法の記事¹

i	1	2	3	4	5	6	7
[S]	0.038	0.194	0.425	0.626	1.253	2.500	3.740
v	0.050	0.127	0.094	0.2122	0.2729	0.2665	0.3317
これを、ミカエリス・メンテン式をモデル関数としてフィッティングします。

$$v = \frac{V_{max}[S]}{K_m + [S]}$$
生化学を学んだことのある人には馴染み深い式ですね。$V_{max}$と$K_m$が最適化すべきパラメータです。
ちなみに、この関数は非線形ですが、変数を$\frac{1}{[S]}$と$\frac{1}{v}$と取ることで
$$\frac{1}{v} = \frac{K_m}{V_{max}}\frac{1}{[S]} + \frac{1}{V_{max}}$$
となり、線形の式になります。実際の実験ではこの式をプロットして、切片と傾きから$V_{max}$と$K_m$を決めます。これをラインウィーバー・バーク プロットと言います。

まず、ヤコビアンとその転置行列を計算します。
ヤコビアン$J$は(データ数$\times$パラメータ数)の行列です。例では、$7\times2$行列です。
Wikipediaの例のようにモデル関数を手計算で微分してもいいですが、ゆくゆくは様々なモデル関数に対応できるようにしたいので、$\Delta \beta = 0.1$として
$$\frac{\partial f(\beta_1)}{\partial \beta_1} = \frac{f(\beta_1 + \Delta \beta) – f(\beta_1)}{\Delta \beta}$$
で数値計算します。

do i=1, m
      func = -(beta1*x0(i))/(beta2 + x0(i))
      deltafunc1 = -((beta1+dbeta)*x0(i))/(beta2 + x0(i))
      deltafunc2 = -(beta1*x0(i))/(beta2 + dbeta + x0(i))
      jacobian(1, i) = (deltafunc1 - func)/dbeta
      jacobian(2, i) = (deltafunc2 - func)/dbeta
      e(i) = y0(i) + func
 end do

dbetaが$\Delta \beta$、deltafunc1が$f(\beta_1 + \Delta \beta)$、deltafunc2が$f(\beta_2 + \Delta \beta)$、e(i)が残差$r$のことです。

次に、$J^TJ$を$A$、$\delta x$を$x$、$J^Tr$を$B$として$Ax = B$の連立方程式を解きます。
行列積、転置はFortran90の組み込み関数であるmatmul()、transpose()で行います。

a = matmul(jacobian, transpose(jacobian))
b = -matmul(e, transpose(jacobian))

右辺$B$がベクトルの実係数をもつ連立一次方程式はLapackのDGESVサブルーチンで計算できます。

call dgesv(n, nrhs, a, lda, ipiv, b, ldb, info)

• 引数

• n : 行列$A$の次数(パラメータ数)
• nrhs : 行列$B$の列数(Bはベクトルなので1)
• a : $A(lda, n)$
• lda : 行列$A$の行数($J^TJ$は$2\times 2$行列なので2)
• ipiv : 次元(n)の整数配列(なんかたぶんルーチン内の計算で使うやつ)
• b : $B(ldb, nrhs)$ 入力には右辺ベクトル$B$、出力時には解行列$x$
• ldb : 行列$B$の行数(データ数)
• info : 0のとき正常終了、それ以外では異常

次に$\beta$を更新します。新しいbetaは
$\beta ← \beta + \delta x$
とします。

beta1 = beta1 + b(1)
beta2 = beta2 + b(2)

これを収束するまで繰り返します。
今回は10回繰り返して収束しているかどうか確認しました。

#実行結果
プログラムの全文はぎっはぶ²

10回ループののち、パラメータは以下のようになりました。
beta1 = 0.3624
beta2 = 0.5595
それをもとにモデル関数を計算し、グラフに当てはめます。

それっぽくデータの点の間を通ってますね。

次は収束具合を見てみましょう。
各繰り返しでのSをプロットしたのが以下です。

2回目の繰り返しで大きく下がっていますね。5回目以降では小数点5桁以下の変化しかありません。収束していると言えます。

ガウス・ニュートン法では初期値を近しい値にしておかなければうまく収束しません。
例えば、beta1とbeta2をそれぞれ2にしてみると10回の繰り返しでは収束しなくなります。
初期値がどれほど離れると収束しなくなるのか、見てみましょう。
beta1 = 2で固定して、beta2を0.01から2まで0.01刻みで変化させた結果が以下のグラフです。

beta2 = 1を超えたあたりで10回では収束しなくなっています。10回以上の繰り返しで収束する可能性もありますが、10回の時点で$\beta$が発散している点もあるため、期待薄です。

#まとめ
ガウス・ニュートン法では、適切な初期値を与えることで、パラメータを最適化できることがわかりました。上記で試したように、ガウス・ニュートン法は初期値が大きく離れていたり、モデル関数の非線形性が強かったりすると収束しません。より広い範囲で収束させる方法として、Marquardt法や共役勾配法があります。

非線形最小二乗法はデータのフィッティングだけでなく、物理や工学分野などで様々な応用がなされています。
例えば、粉末X線回折ピークを分析する手法であるリートベルト解析は、数学的には非線形最小二乗法と同義です。

最後に、参考にした記事、書籍を勝手ながら紹介させていただきます。

• MATLABで非線形最小二乗フィッティングする手順メモ
  : MATLABで同じWikipediaの例のパラメータを最適化していらっしゃいます。

• 行列計算パッケージLAPACK利用の手引　(小国力 著　丸善)
  : LAPACKの公式リファレンスはぜんぶ英語なので、日本語の解説はとても助かりました。特に、今回初めてLapackを使ったので、Lapackはなんぞや、というところから参考にしました。

以上です。ご一読ありがとうございます。

#参考文献

• https://ja.wikipedia.org/wiki/ガウス・ニュートン法
• https://ossyaritoori.hatenablog.com/entry/2017/01/31/MATLABで非線形最小二乗フィッティングする手順メモ_
• 行列計算パッケージLAPACK利用の手引　小国力 著　丸善

#コード

• https://github.com/PicricAcid/gauss_newton

1. https://ja.wikipedia.org/wiki/ガウス・ニュートン法
   の例を借ります。
   酵素の基質濃度と反応速度の関係のデータです。 ↩

2. https://github.com/PicricAcid/gauss_newton
   を見ていただくとして、プログラムを実行してみましょう。
   とりあえず、初期値はWikipediaの値
   beta1 = 0.9
   beta2 = 0.2
   とします。 ↩