ガウス・ニュートン法

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

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

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

## 概要
データフィッティングにおいて、与えられたモデル関数 *y*&nbsp;=&nbsp;*f* (*x* , ***&beta;***) が*m* 個のデータ点 {(*x<sub>i</sub>* , *y<sub>i</sub>* ); *i* = 1, ... , *m* } に最もよくフィットするような*n* （&le; *m* ）個のパラメータ***&beta;*** = (&beta;1 , ... , &beta;*n* )を見つけることが目的である。

このとき、残差を
: r_i(\boldsymbol{\beta})= y_i - f(x_i, \boldsymbol{\beta})
とする。

このとき、ガウス・ニュートン法は残差の平方和
: S(\boldsymbol \beta)= \sum_{i=1}^m (r_i(\boldsymbol \beta))^2
の最小値を反復計算で求める。初期推測値***&beta;***(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}
は***&beta;***(*s* ) における***r*** のヤコビアン、*J<sub>*'r**</sub>''<sup>T</sup> は行列''J<sub>**r**</sub>'' の[[転置行列|転置]]を表す。

''m''&nbsp;=&nbsp;''n'' ならば、この反復計算は
:<math> \boldsymbol{\beta}^{(s+1)} = \boldsymbol{\beta}^{(s)} - J_r^{-1} \boldsymbol{r}(\boldsymbol \beta^{(s)}) </math>
のように簡略化される。これは1次元ニュートン法の直接的な一般化である。

ガウス・ニュートン法は関数''f'' の[[ヤコビアン]]''J<sub>f</sub>'' を用いて次のように表すこともできる：
:<math> \boldsymbol{\beta}^{(s+1)} = \boldsymbol{\beta}^{(s)} + ({J_f}^\mathrm{T} {J_f} )^{-1} {J_f}^\mathrm{T} \boldsymbol{r}(\boldsymbol \beta^{(s)}). </math>

== 注釈 ==
ガウス・ニュートン法は関数''r<sub>i</sub>'' を並べたベクトルの[[線形近似]]で与えられる。[[テイラーの定理]]を用いれば、各反復において次式が成り立つ：
: <math>\boldsymbol{r}(\boldsymbol{\beta})\approx \boldsymbol{r}(\boldsymbol{\beta}^s)+J_r(\boldsymbol{\beta}^s)\boldsymbol{\Delta}.</math>
ここで **&Delta;** = **&beta;** - **&beta;**<sup>''s''</sup> である。右辺の残差平方和を最小化する**&Delta;**を見つけること、すなわち
: <math>\operatorname{min}\left[{\|\boldsymbol{r}(\boldsymbol{\beta}^s) + J_r(\boldsymbol{\beta}^s)\boldsymbol{\Delta}\|_2}^2\right]</math>
は線形[[最小二乗法]]の問題であるため、[[陽関数|陽的]]に解くことができ、正規方程式を与える。ここで ||*||<sub>2</sub> は 2-ノルム（[[ユークリッドノルム]]）である。

正規方程式は未知の増分**&Delta;**についての''m'' 本の線形同次方程式である。これは[[コレスキー分解]]を用いることで、またはより良い方法としては''J<sub>**r**</sub>'' の[[QR分解]]を用いることで、1ステップで解ける。大きな系に対しては、[[共役勾配法]]のような[[反復解法]]が有効である。''J<sub>**r**</sub>'' の列ベクトルが[[線形従属]]である場合、''J<sub>**r**</sub>''<sup>T</sup>''J<sub>**r**</sub>'' が非正則になるため反復解法は失敗する。

==例==

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

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

これらのデータに対し、次の形のモデル曲線（[[ミカエリス・メンテン式]]）のパラメータ''V''<sub>max</sub> と''K''<sub>M</sub> を、最小二乗の意味で最もよくフィットするように決定したい：
:<math>v=\frac{V_\mathrm{max}[\mathrm{S}]}{K_\mathrm{M}+[\mathrm{S}]}</math>

''x<sub>i</sub>'' と''y<sub>i</sub>'' (''i'' = 1, ... , 7) で [S] と''v'' のデータを表す。また、''&beta;''<sub>1</sub> = ''V''<sub>max</sub> 、''&beta;''<sub>2</sub> = ''K<sub>M</sub>'' とする。残差
: <math>r_i = y_i - \frac{\beta_1 x_i}{\beta_2+x_i} \quad (i=1,\dots, 7)</math>
の平方和を最小化する''&beta;''<sub>1</sub> と''&beta;''<sub>2</sub> を見つけることが目的となる。

未知パラメータ***&beta;***に関する残差ベクトル***r*** のヤコビアン''J<sub>r</sub>'' は7&times;2の行列で、''i'' 番目の行は次の要素を持つ：
:<math>\begin{pmatrix}\dfrac{\partial r_i}{\partial \beta_1},& \dfrac{\partial r_i}{\partial \beta_2}\end{pmatrix}
= \begin{pmatrix}-\dfrac{x_i}{\beta_2+x_i},& \dfrac{\beta_1x_i}{\left(\beta_2+x_i\right)^2}\end{pmatrix}.</math>

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

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

ガウス・ニュートン法のは2次である。もし初期推測値が最小値から遠いか、または行列''J<sub>r</sub>''<sup>T</sup> ''J<sub>r</sub>'' が[[条件数|悪条件]]であれば収束は遅いか、全くしなくなる。例えば、''m'' = 2本の方程式と''n'' = 1個の変数のある次の問題を考える：
:<math> \begin{align}
r_1(\beta) &= \beta + 1 \\
r_2(\beta) &= \lambda \beta^2 + \beta - 1.
\end{align} </math>
この問題の最適値は&beta; = 0 である。もし&lambda; = 0 なら実質的に線形問題であり、最適値は一回の計算で見つかる。もし|&lambda;| < 1 なら、この手法は線形に収束し残差は係数|&lambda;|で反復ごとに漸近的に減少する。しかし|&lambda;| > 1 なら、この方法はもはや局所的にも収束しない。

== ニュートン法からの導出 ==
後で示すように、ガウス・ニュートン法は近似関数の最適化にニュートン法を適用して導出される。そのことから、ガウス・ニュートン法の収束の速さはほとんど2次である。

パラメータ***&beta;***を持つ関数''S'' の最小化をするとき、ニュートン法による[[漸化式]]は
:<math> \boldsymbol\beta^{(s+1)} = \boldsymbol\beta^{(s)} - H^{-1} \boldsymbol{g} \, </math>
である。ここで***g*** は''S'' の[[勾配ベクトル]]、''H'' は''S'' の[[ヘッシアン]]である。<math> S = \sum_{i=1}^m r_i^2</math>であるから、勾配***g*** は次で与えられる：
:<math>\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}.</math>
ヘッシアン''H'' は勾配***g*** を***&beta;*** で微分することで計算される：
:<math>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).</math>

2階微分項（右辺第2項）を無視することでガウス・ニュートン法を得る。つまり、ヘッシアンは
:<math>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}</math>
と近似される。ここで
:<math>J_{ij}=\frac{\partial r_i}{\partial \beta_j}</math>
はヤコビアン''J<sub>r</sub>'' の成分である。

これらの表現を上述の漸化式に代入して、次式を得る：
:<math> \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}. </math>

ただしガウス・ニュートン法の収束は常に保証されているわけではない。2階微分項を無視するという近似、すなわち
:<math>\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|</math>
に正当性があるのは次の2つの条件の下であり、これらが成り立つ場合には収束が期待される：
# ''r<sub>i</sub>'' は十分小さい。少なくとも最小値付近。
# 関数の非線形性は穏やかであり、<math>{\partial^2 r_i}/{\partial \beta_j \partial \beta_k}</math> が比較的小さくなる。

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

''S'' は反復ごとに必ずしも減少するわけではないが、増分ベクトル**&Delta;**は''S'' が減少する方向を向いているから、''S'' (***&beta;***<sup>s</sup>) が停留点にない限り、任意の十分に小さな''&alpha;''> 0 に対して ''S'' (***&beta;***<sup>s</sup> + ''&alpha;***&Delta;*') s) が成り立つ。したがって、発散したときに、更新方程式に**縮小因子***&alpha;*を導入して
: \boldsymbol \beta^{s+1} = \boldsymbol \beta^s+\alpha\boldsymbol{\Delta}
とすることが解決法の一つとなる。

言い換えれば、増分ベクトル**&Delta;**は目的関数*S* の下り方向を指してはいるが長すぎるので、その道のほんの一部を行くことで、*S* を減少させようというアイディアである。縮小因子*&alpha;*の最適値は直線探索で見つけることができる。つまり、直接探索法を（通常 0 (J^\mathrm{T} J+\lambda D)\Delta = J^\mathrm{T} \boldsymbol{r},
ここで*D* は正定値対角行列である。&lambda;が 0 から増大するにつれて増分ベクトル**&Delta;**は長さが単調に減少し、かつ方向は最急降下方向に近づくため、&lambda;を十分大きくすれば必ずより小さい*S* の値を見出せることが保証されている。

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

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

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

## 脚注

### 注釈

### 出典

## 参考文献
*
* .
*
*