ガウス・ニュートン法

出典: フリー百科事典『ウィキペディア(Wikipedia)』
移動: 案内検索

ガウス・ニュートン法(ガウス・ニュートンほう、: Gauss-Newton method)とは、非線形最小二乗法を解く手法の一つである。これは関数の最大・最小値を見出すニュートン法の修正とみなすことができる。ニュートン法とは違い、ガウス・ニュートン法は二乗和の最小化にしか用いることができないが、計算するのが困難な2階微分が不要と言う長所がある。

非線形最小二乗法は非線形回帰英語版などで、観測データを良く表すようにモデルのパラメータを調整するために必要となる。

この手法の名称はカール・フリードリヒ・ガウスアイザック・ニュートンにちなむ。

概要[編集]

データフィッティングにおいて、与えられたモデル関数 y = f (x , β) がm 個のデータ点 {(xi , yi ); i = 1, ... , m } に最もよくフィットするようなn (≤ m )個[1]のパラメータβ = (β1 , ... , βn )を見つけることが目的である。

このとき、残差

r_i(\boldsymbol{\beta})= y_i - f(x_i, \boldsymbol{\beta})

とする。

このとき、ガウス・ニュートン法は残差の平方和

 S(\boldsymbol \beta)= \sum_{i=1}^m (r_i(\boldsymbol \beta))^2

の最小値を反復計算で求める[2]。初期推測値β(0) から初めて、この方法は以下の計算を繰り返す。

 \boldsymbol \beta^{(s+1)} = \boldsymbol \beta^{(s)} - ({J_r}^\mathrm{T} {J_r} )^{-1} {J_r}^\mathrm{T} \boldsymbol{r}(\boldsymbol \beta^{(s)}),\quad (s=0,1,2,\dots).

ここで

 J_r = \frac{\partial r_i (\boldsymbol{\beta}^{(s)})}{\partial \beta_j}

β(s ) におけるrヤコビアンJrT は行列Jr転置を表す。

m = n ならば、この反復計算は

 \boldsymbol{\beta}^{(s+1)} = \boldsymbol{\beta}^{(s)} - J_r^{-1} \boldsymbol{r}(\boldsymbol \beta^{(s)})

のように簡略化される。これは1次元ニュートン法の直接的な一般化である。

ガウス・ニュートン法は関数fヤコビアンJf を用いて次のように表すこともできる:

 \boldsymbol{\beta}^{(s+1)} = \boldsymbol{\beta}^{(s)} + ({J_f}^\mathrm{T} {J_f} )^{-1} {J_f}^\mathrm{T} \boldsymbol{r}(\boldsymbol \beta^{(s)}).

注釈[編集]

ガウス・ニュートン法は関数ri を並べたベクトルの線形近似で与えられる。テイラーの定理を用いれば、各反復において次式が成り立つ:

\boldsymbol{r}(\boldsymbol{\beta})\approx \boldsymbol{r}(\boldsymbol{\beta}^s)+J_r(\boldsymbol{\beta}^s)\boldsymbol{\Delta}.

ここで Δ = β - βs である。右辺の残差平方和を最小化するΔを見つけること、すわなち

\operatorname{min}\left[{\|\boldsymbol{r}(\boldsymbol{\beta}^s) + J_r(\boldsymbol{\beta}^s)\boldsymbol{\Delta}\|_2}^2\right]

は線形最小二乗法の問題であるため、陽的に解くことができ、正規方程式を与える。ここで ||*||2 は 2-ノルム(ユークリッドノルム)である。

正規方程式は未知の増分Δについてのm 本の線形同次方程式である。これはコレスキー分解を用いることで、またはより良い方法としてはJrQR分解を用いることで、1ステップで解ける。大きな系に対しては、共役勾配法のような反復解法が有効である。Jr の列ベクトルが線形従属である場合、JrTJr が非正則になるため反復解法は失敗する。

[編集]

この例で得られているデータ(赤点)と、\hat\beta_1=0.362, \hat\beta_2=0.556から計算されたモデル曲線(青線)

ここでは例として、ガウス・ニュートン法を使ってデータとモデルによる予測値の間の残差平方和を最小化し、データにモデルをフィットさせる。

生物実験において、酵素媒介反応における基質濃度[S] と反応率v の関係として次表のようなデータが得られたとする(右図の赤点)。

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

これらのデータに対し、次の形のモデル曲線(ミカエリス・メンテン式)のパラメータVmaxKM を、最小二乗の意味で最もよくフィットするように決定したい[3]

v=\frac{V_\mathrm{max}[\mathrm{S}]}{K_\mathrm{M}+[\mathrm{S}]}

xiyi (i = 1, ... , 7) で [S] とv のデータを表す。また、β1 = Vmax 、β2 = KM とする。残差

r_i = y_i - \frac{\beta_1 x_i}{\beta_2+x_i} \quad (i=1,\dots, 7)

の平方和を最小化するβ1 とβ2 を見つけることが目的となる。

未知パラメータβに関する残差ベクトルr のヤコビアンJr は7×2の行列で、i 番目の行は次の要素を持つ:

\begin{pmatrix}\frac{\partial r_i}{\partial \beta_1},& \frac{\partial r_i}{\partial \beta_2}\end{pmatrix}
= \begin{pmatrix}-\frac{x_i}{\beta_2+x_i},& \frac{\beta_1x_i}{\left(\beta_2+x_i\right)^2}\end{pmatrix}.

初期推定値としてβ1 = 0.9、β2 = 0.2から始め、ガウス・ニュートン法による5回の反復計算を行うと、最適値 \hat\beta_1=0.362, \hat\beta_2=0.556 が得られる。残差平方和は5回の反復計算で初期の1.445から0.00784まで減少する。右図はこれらの最適パラメータを用いたモデルで決まる曲線と、実験データとの比較を示す。

収束性[編集]

増分ΔS の減少方向を向いていることは証明されている[4]。もしこのアルゴリズムが収束すれば、その極限はS停留点である。しかし収束については、ニュートン法では保証されている局所収束さえも保証されていない。

ガウス・ニュートン法の収束の速さ英語版は2次である[5]。もし初期推測値が最小値から遠いか、または行列JrT Jr悪条件であれば収束は遅いか、全くしなくなる。例えば、m = 2本の方程式とn = 1個の変数のある次の問題を考える:

 \begin{align}
r_1(\beta) &= \beta + 1 \\
r_2(\beta) &= \lambda \beta^2 + \beta - 1.
\end{align}

この問題の最適値はβ = 0 である。もしλ = 0 なら実質的に線形問題であり、最適値は一回の計算で見つかる。もし|λ| < 1 なら、この手法は線形に収束し残差は係数|λ|で反復ごとに漸近的に減少する。しかし|λ| > 1 なら、この方法はもはや局所的にも収束しない[6]

ニュートン法からの導出[編集]

後に示すように、ガウス・ニュートン法は近似関数の最適化に用いられるニュートン法から与えられる。その結果、ガウス・ニュートン法の収束の速さはほとんど2次である。

パラメータβを持つ関数S の最小化をするとき、ニュートン法による漸化式

 \boldsymbol\beta^{(s+1)} = \boldsymbol\beta^{(s)} - H^{-1} \boldsymbol{g} \,

である。ここでgS勾配ベクトルHSヘッシアンである。 S = \sum_{i=1}^m r_i^2であるから、勾配g は次で与えられる:

\boldsymbol{g}=2{J_r}^\mathrm{T} \boldsymbol{r}, \quad\text{or,}\quad g_j=2\sum_{i=1}^m r_i\frac{\partial r_i}{\partial \beta_j}.

ヘッシアンH は勾配gβ で微分することで計算される:

H_{jk}=2\sum_{i=1}^m \left(\frac{\partial r_i}{\partial \beta_j}\frac{\partial r_i}{\partial \beta_k}+r_i\frac{\partial^2 r_i}{\partial \beta_j \partial \beta_k} \right).

2階微分項(右辺第2項)を無視することでガウス・ニュートン法を得る。つまり、ヘッシアンは

H\approx 2{J_r}^\mathrm{T} J_r, \quad\text{or,}\quad H_{jk}\approx 2\sum_{i=1}^m \frac{\partial r_i}{\partial \beta_j}\frac{\partial r_i}{\partial \beta_k}
 =2\sum_{i=1}^m J_{ij}J_{ik}

と近似される。ここで

J_{ij}=\frac{\partial r_i}{\partial \beta_j}

はヤコビアンJr の成分である。

これらの表現を上述の漸化式に代入して、次式を得る:

 \boldsymbol{\beta}^{(s+1)} = \boldsymbol\beta^{(s)}+\boldsymbol{\Delta};\quad
\boldsymbol{\Delta} = -({J_r}^\mathrm{T} J_r)^{-1} {J_r}^\mathrm{T} \boldsymbol{r}.

ガウス・ニュートン法の収束は常に保証されているわけではない。2階微分項を無視するという近似、すなわち

\left|r_i\frac{\partial^2 r_i}{\partial \beta_j \partial \beta_k}\right| \ll \left|\frac{\partial r_i}{\partial \beta_j}\frac{\partial r_i}{\partial \beta_k}\right|

に正当性があるのは次の2つの条件のもとでであり、これらが成り立つ場合には収束が期待される[7]

  1. ri は十分小さい。少なくとも最小値付近。
  2. 関数の非線形性は穏やかであり、{\partial^2 r_i}/{\partial \beta_j \partial \beta_k} が比較的小さくなる。

改善バージョン[編集]

ガウス・ニュートン法は、初期推定値が真の解から大きく離れていたり、モデル関数の非線形性が大きい場合には安定性が悪い。また、残差平方和S は反復ごとに必ずしも減少するわけではない。そのため、実用上は安定化が必要である[8]

S は反復ごとに必ずしも減少するわけではないが、増分ベクトルΔS が減少する方向を向いているから、S (βs) が停留点にない限り、任意の十分に小さなα> 0 に対して S (βs + αΔ) < S (βs) が成り立つ。したがって、発散したときに、更新方程式に縮小因子[9]αを導入して

 \boldsymbol \beta^{s+1} = \boldsymbol \beta^s+\alpha\boldsymbol{\Delta}

とすることが解決法の一つとなる。

言い換えれば、増分ベクトルΔは目的関数S の下り方向を指してはいるが長すぎるので、その道のほんの一部を行くことで、S を減少させようというアイディアである。縮小因子αの最適値は直線探索で見つけることができる。つまり、直接探索法を(通常 0 < α ≤ 1 の区間で)用いて、S を最小化する値を探すことでαの大きさは決められる。

最適な縮小因子αが 0 に近いような場合、発散を回避する別の方法はLevenberg-Marquardtアルゴリズム英語版信頼範囲英語版法)を使うことである[2]。増分ベクトルが最急降下方向に向くように正規方程式は修正される。

(J^\mathrm{T} J+\lambda D)\Delta = J^\mathrm{T} \boldsymbol{r},

ここでD正定値対角行列である。λが 0 から増大するにつれて増分ベクトルΔは長さが単調に減少し、かつ方向は最急降下方向に近づくため、λを十分大きくすれば必ずより小さいS の値を見出せることが保証されている[10]

いわゆるMarquardtパラメータλは直線探索により最適化されるが、λが変わるたびに毎回シフトベクトルの再計算をしなければならないため非効率的である。より効率的な方法は、発散が起きた時にλをS が減少するまで増加させる。そしてその値を1回の反復から次まで維持する、しかしλを 0 に設定することができるときにもしカットオフ値に届くまで可能なら減少させる。このときS の最小値は標準のガウス・ニュートン法の最小化になる。

関連するアルゴリズム[編集]

DFP法英語版BFGS法英語版のような擬似ニュートン法英語版では、ヘッシアン{\partial^2 S}/{\partial \beta_j \partial\beta_k}の推定は1階微分{\partial r_i}/{\partial\beta_j}のみを用いて数値的になされる。したがってn 回の反復計算による修正の後、この方法はパフォーマンスにおいてニュートン法を近似する。ガウス・ニュートン法やLevenberg-Marquardt法などは非線形最小二乗問題にのみ適用できるのに対して、擬似ニュートン法は一般的な実数値関数を最小化できることに注意する。

1階微分を使って最小化問題を解く他の方法は、最急降下法のみである。しかし、この方法は2階微分を近似的にも考慮しておらず、その結果、この方法は多くの関数に対して、特にパラメータが強い相互作用を持っている場合は非常に非効率的である。

脚注[編集]

  1. ^ アルゴリズム内のm ≥ n という仮定は必要である。そうでなければ、行列JrTJr の逆行列を計算できず、正規方程式の解(少なくとも唯一解)を求めることができない。
  2. ^ a b Björck (1996)
  3. ^ ミカエリス・メンテン式#実験によるパラメータの決定で説明するように、実際は変数に[S]-1v-1 を選ぶことで、この問題は線形最小二乗法として解ける。
  4. ^ Björck (1996) p260
  5. ^ Björck (1996) p341, 342
  6. ^ Fletcher (1987) p.113
  7. ^ Nocedal (1997) [要ページ番号]
  8. ^ 中川、小柳 (1982) p.98
  9. ^ 中川、小柳 (1982) p.98
  10. ^ 中川、小柳 (1982) p.102

参考文献[編集]

  • Björck, A. (1996). Numerical methods for least squares problems. SIAM, Philadelphia. ISBN 0-89871-360-9. 
  • Fletcher, Roger (1987). Practical methods of optimization (2nd ed.). New York: John Wiley & Sons. ISBN 978-0-471-91547-8. .
  • Nocedal, Jorge; Wright, Stephen (1999). Numerical optimization. New York: Springer. ISBN 0-387-98793-2. 
  • 中川徹; 小柳義夫 『最小二乗法による実験データ解析』 東京大学出版会、1982年ISBN 4-13-064067-4