備忘録。非線形最小二乗法をニュートン法で解くだけの記事。


・概要

 \(n\)個のデータの組\((x_1,y_1)...(x_n,y_n)\)を、\(m\)個のパラメーター\(\beta_1...\beta_m\)を持つ関数\(f(x;\boldsymbol{\beta})\)にフィッティングさせたい。ここで$$\boldsymbol{\beta}=\left[\begin{array}{cc}\beta_1\\\vdots\\\beta_m\end{array}\right]$$である。
 パラメータを持つ関数というのは、例えば正規分布における\(\mu\) や \(\sigma\)、ミカエリス・メンテン式における\(V_{\mathrm{max}}\) や \(K_{\mathrm m}\) のようなものである。別に二次関数 \(ax^2+bx+c\) の \(a,b,c\) のようなものでもいい。とにかくグラフの形を特徴づける値である。

・方針

 目標は、残差$$ S(\boldsymbol{\beta})=\sum_i (y_i-f(x_i;\boldsymbol{\beta}))^2$$を最小化することである。ここでデータの(フィッティング曲線からの)ばらつき具合はすべて等しいものだと仮定している。
 残差が最小ということは、必要条件として$$ \frac{\partial}{\partial\boldsymbol{\beta}}S(\boldsymbol{\beta})=0\quad\ldots(1)$$であることが求められる。
 ここで、その条件を満たす真の値 \(\boldsymbol\beta\) の近傍に初期値 \(\boldsymbol\beta_0\) を取る。すなわち、ある小さい誤差 \(\boldsymbol\delta\) があって、\(\boldsymbol\beta=\boldsymbol\beta_0+\boldsymbol\delta\) であるものとする。そこで \(f(x;\boldsymbol{\beta})\) を \(\boldsymbol\beta_0\) まわりでテイラー展開し、高次項を無視すると$$f(x;\boldsymbol{\beta})\approx f(x;\boldsymbol{\beta_0})+\boldsymbol{g}\boldsymbol{\delta}\quad\ldots(2)$$
というように \(\boldsymbol{\delta}\) に関する1次の式にすることができる。ただし、めんどくさいので \(f\) の勾配(grad) \(\displaystyle \frac{\partial}{\partial\boldsymbol{\beta}}f(x;\boldsymbol{\beta})\) は簡略化して \(\boldsymbol{g}\) とおいた。以降も同様。
 なお(2)式の右辺第2項は \(\boldsymbol{g}\) が行ベクトル、\(\boldsymbol{\delta}\) が列ベクトルであり、項全体としてはスカラーになると思うと都合がいい(内積を取っていると思えばいい)。成分表記をすれば \(\displaystyle \sum_j g_j\delta_j\)である。

・計算

 (2)の近似を(1)式にぶちこんでみると、$$ \begin{eqnarray}\frac{\partial}{\partial\boldsymbol{\beta}}S(\boldsymbol{\beta})&=&-2\sum_i (y_i-f(x_i;\boldsymbol{\beta}))\boldsymbol{g}\\
&=&-2\sum_i (y_i-f(x_i;\boldsymbol{\beta_0})-\boldsymbol{g}\boldsymbol{\delta})\boldsymbol{g}\\
&=&-2\sum_i ( (y_i-f(x_i;\boldsymbol{\beta_0}))\boldsymbol{g}-\boldsymbol{g}\boldsymbol{\delta}\boldsymbol{g})\quad\ldots(3)
\end{eqnarray}
$$となってごちゃごちゃなのでまとめてみよう。
 まず第1項はスカラー×ベクトルなのでこのまま適当に$$ \sum_i (y_i-f(x_i;\boldsymbol{\beta_0}))^t\!\boldsymbol{g}=\boldsymbol{\gamma}$$とでも置いてしまう。ただし\(^t\)は転置である。いま転置を付けたのは後々のためである。
 続いて第2項の \(\boldsymbol{g}\boldsymbol{\delta}\boldsymbol{g}\) だが、行ベクトル×列ベクトル×行ベクトルの形であり、前半の行ベクトル×列ベクトルは内積だと思ってよいので順序を入れ替えることができて、\(^t\!\boldsymbol{\delta}\ ^t\!\boldsymbol{g}\boldsymbol{g}\) と書き換えることができる。また行列の積の結合順序は任意であるので \(^t\!\boldsymbol{g}\boldsymbol{g}\) を先に計算してしまってよい。従って、\(\displaystyle \sum_i {}^t\!\boldsymbol{g}\boldsymbol{g}=\boldsymbol{\alpha}\) (これは\(m\times m\)対称行列になる)とおく。成分表記すれば \(\displaystyle \alpha_{jk}=\sum_i g_jg_k\) になる。

 そうすると(3)式は$$\frac{\partial}{\partial\boldsymbol{\beta}}S(\boldsymbol{\beta})=-2(^t\!\boldsymbol{\gamma}-{}^t\!\boldsymbol{\delta}\boldsymbol{\alpha})$$となる。これが\(0\)に等しいので、$$\begin{eqnarray}&&\quad^t\!\boldsymbol{\gamma}-{}^t\!\boldsymbol{\delta}\boldsymbol{\alpha}&=&0\\
&\Leftrightarrow&\quad\quad\quad ^t\!\boldsymbol{\delta}\boldsymbol{\alpha}&=&^t\!\boldsymbol{\gamma}\\
&\Leftrightarrow&\quad\quad\quad \boldsymbol{\alpha}\boldsymbol{\delta}&=&\boldsymbol{\gamma}\\
\end{eqnarray}$$である。最後は両辺を転置しているが、\(\boldsymbol{\alpha}\)は対称行列であるので \(^t\!\boldsymbol{\alpha}=\boldsymbol{\alpha}\)である。

 あとは両辺に左から \(\boldsymbol{\alpha}^{-1}\) を掛けてやれば晴れて \(\boldsymbol{\delta}=\boldsymbol{\alpha}^{-1}\boldsymbol{\gamma}\) と求まる。こうして求まった \(\boldsymbol{\beta}\approx\boldsymbol{\beta_0}+\boldsymbol{\delta}\) を再度初期値に設定して同じことをやることで、限りなく真の値に近づいていくという寸法である。ただしニュートン法はあまり収束が早くないうえに、発散することもあるので万能ではない。


 ヤコビアン?そんな高尚なものなど知らん。

非線形最小二乗法

投稿ナビゲーション


コメントを残す

メールアドレスが公開されることはありません。