ガウス・ニュートン法のソースを表示

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

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

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

== 概要 ==
[[曲線あてはめ|データフィッティング]]において、与えられたモデル関数 ''y''&nbsp;=&nbsp;''f'' (''x'' , '''''&beta;''''') が''m'' 個のデータ点 {(''x<sub>i</sub>'' , ''y<sub>i</sub>'' ); ''i'' = 1, ... , ''m'' } に最もよくフィットするような''n'' （&le; ''m'' ）個<ref>アルゴリズム内の''m''&nbsp;&ge;&nbsp;''n'' という仮定は必要である。そうでなければ、行列''J<sub>r</sub>''<sup>T</sup>''J<sub>r</sub>'' の逆行列を計算できず、正規方程式の解（少なくとも唯一解）を求めることができない。</ref>のパラメータ'''''&beta;''''' = (&beta;<sub>1</sub> , ... , &beta;<sub>''n''</sub> )を見つけることが目的である。

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

このとき、ガウス・ニュートン法は残差の平方和
:<math> S(\boldsymbol \beta)= \sum_{i=1}^m (r_i(\boldsymbol \beta))^2</math>
の最小値を[[反復計算]]で求める<ref name=ab>Björck (1996)</ref>。初期推測値'''''&beta;'''''<sup>(0)</sup> から初めて、この方法は以下の計算を繰り返す。
:<math> \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). </math>
ここで
:<math> J_r = \frac{\partial r_i (\boldsymbol{\beta}^{(s)})}{\partial \beta_j}</math>
は'''''&beta;'''''<sup>(''s'' )</sup> における'''''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>'' が非正則になるため反復解法は失敗する。

==例==
[[File:Gauss Newton illustration.png|thumb|right|280px|この例で得られているデータ（赤点）と、<math>\hat\beta_1=0.362, \hat\beta_2=0.556</math>から計算されたモデル曲線（青線）]]
ここでは例として、ガウス・ニュートン法を使ってデータとモデルによる予測値の間の残差平方和を最小化し、データにモデルをフィットさせる。

生物実験において、[[酵素媒介反応]]<!--enzyme-mediated reaction-->における[[基質]]濃度[S] と[[反応率]]<!--reaction rate-->''v'' の関係として次表のようなデータが得られたとする（右図の赤点）。
:{|class="wikitable" style="text-align: center;"
|''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''<sub>max</sub> と''K''<sub>M</sub> を、最小二乗の意味で最もよくフィットするように決定したい<ref>[[ミカエリス・メンテン式#定数の決定]]で説明するように、実際は変数に[S]<sup>-1</sup>と''v''<sup>-1</sup> を選ぶことで、この問題は線形最小二乗法として解ける。</ref>：
:<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'' の減少方向を向いていることは証明されている<ref>Björck (1996) p260</ref>。もしこのアルゴリズムが収束すれば、その極限は''S'' の[[極値|停留点]]<!--stationary point-->である。しかし収束については、ニュートン法では保証されている局所収束さえも保証されていない。

ガウス・ニュートン法の{{仮リンク|収束の速さ|en|rate of convergence}}は2次である<ref>Björck (1996) p341, 342</ref>。もし初期推測値が最小値から遠いか、または行列''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;|で反復ごとに漸近的に減少する<!--If |λ| < 1 then the method converges linearly and the error decreases asymptotically with a factor |λ| at every iteration.-->。しかし|&lambda;| > 1 なら、この方法はもはや局所的にも収束しない<ref>Fletcher (1987) p.113</ref>。

== ニュートン法からの導出 ==
後で示すように、ガウス・ニュートン法は近似関数の最適化にニュートン法を適用して導出される<!--In what follows, the Gauss–Newton algorithm will be derived from [[Newton's method in optimization|Newton's method]] for function optimization via an approximation.-->。そのことから、ガウス・ニュートン法の収束の速さはほとんど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つの条件の下であり、これらが成り立つ場合には収束が期待される<ref>Nocedal (1997) {{Page needed|date=December 2010}}</ref>：
# ''r<sub>i</sub>'' は十分小さい。少なくとも最小値付近。
# 関数の非線形性は穏やかであり、<math>{\partial^2 r_i}/{\partial \beta_j \partial \beta_k}</math> が比較的小さくなる。

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

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

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

最適な縮小因子''&alpha;''が 0 に近いような場合、発散を回避する別の方法は[[レーベンバーグ・マルカート法]]（[[信頼領域]]法）を使うことである<ref name="ab"/>。増分ベクトルが[[最急降下法|最急降下]]方向に向くように正規方程式は修正される。
:<math>(J^\mathrm{T} J+\lambda D)\Delta = J^\mathrm{T} \boldsymbol{r},</math>
ここで''D'' は[[正定値]]対角行列である。&lambda;が 0 から増大するにつれて増分ベクトル'''&Delta;'''は長さが単調に減少し、かつ方向は最急降下方向に近づくため、&lambda;を十分大きくすれば必ずより小さい''S'' の値を見出せることが保証されている<ref>中川、小柳 (1982) p.102</ref>。

いわゆるマルカートパラメータ&lambda;は直線探索により最適化されるが、&lambda;が変わるたびに毎回シフトベクトルの再計算をしなければならないため非効率的である。より効率的な方法は、発散が起きた時に&lambda;を''S'' が減少するまで増加させる。そしてその値を1回の反復から次まで維持する、しかし&lambda;を 0 に設定することができるときにもしカットオフ値に届くまで可能なら減少させる。このとき''S'' の最小値は標準のガウス・ニュートン法の最小化になる<!--but decrease it if possible until a cut-off value is reached when the Marquardt parameter can be set to zero; the minimization of ''S'' then becomes a standard Gauss–Newton minimization.-->。

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

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

== 脚注 ==
{{脚注ヘルプ}}
{{Reflist}}

==参考文献==
*{{cite book | last = Björck | first = A. | title = Numerical methods for least squares problems | publisher = SIAM, Philadelphia | year = 1996 | isbn = 0-89871-360-9 }}
* {{Cite book | last1=Fletcher | first1=Roger | title=Practical methods of optimization | publisher=John Wiley & Sons | location=New York | edition=2nd | isbn=978-0-471-91547-8 | year=1987 }}.
*{{cite book | last = Nocedal | first = Jorge | coauthors = Wright, Stephen | title = Numerical optimization | publisher = New York: Springer | year = 1999 | isbn = 0-387-98793-2 }}
*{{cite book |和書 |author=中川徹|author2=[[小柳義夫]] |title=最小二乗法による実験データ解析 |publisher=東京大学出版会 |year=1982 |isbn=4-13-064067-4 }}

{{最適化アルゴリズム}}

{{DEFAULTSORT:かうすにゆうとんほう}}
[[Category:最適化アルゴリズム]]
[[Category:統計アルゴリズム]]
[[Category:アイザック・ニュートン]]
[[Category:カール・フリードリヒ・ガウス]]
[[Category:数学のエポニム]]
[[Category:最小二乗法]]
[[Category:数学に関する記事]]