ガウス・ニュートン法

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

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

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

概要

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

このとき、残差を

${\displaystyle r_{i}({\boldsymbol {\beta }})=y_{i}-f(x_{i},{\boldsymbol {\beta }})}$

とする。

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

${\displaystyle S({\boldsymbol {\beta }})=\sum _{i=1}^{m}(r_{i}({\boldsymbol {\beta }}))^{2}}$

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

${\displaystyle {\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 ).}$

ここで

${\displaystyle J_{r}={\frac {\partial r_{i}({\boldsymbol {\beta }}^{(s)})}{\partial \beta _{j}}}}$

はβ^{(s )} におけるr のヤコビアン、J_r^T は行列J_r の転置を表す。

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

${\displaystyle {\boldsymbol {\beta }}^{(s+1)}={\boldsymbol {\beta }}^{(s)}-J_{r}^{-1}{\boldsymbol {r}}({\boldsymbol {\beta }}^{(s)})}$

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

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

${\displaystyle {\boldsymbol {\beta }}^{(s+1)}={\boldsymbol {\beta }}^{(s)}+({J_{f}}^{\mathrm {T} }{J_{f}})^{-1}{J_{f}}^{\mathrm {T} }{\boldsymbol {r}}({\boldsymbol {\beta }}^{(s)}).}$

注釈

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

${\displaystyle {\boldsymbol {r}}({\boldsymbol {\beta }})\approx {\boldsymbol {r}}({\boldsymbol {\beta }}^{s})+J_{r}({\boldsymbol {\beta }}^{s}){\boldsymbol {\Delta }}.}$

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

${\displaystyle \operatorname {min} \left[{\|{\boldsymbol {r}}({\boldsymbol {\beta }}^{s})+J_{r}({\boldsymbol {\beta }}^{s}){\boldsymbol {\Delta }}\|_{2}}^{2}\right]}$

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

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

例

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

生物実験において、酵素媒介反応における基質濃度[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

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

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

x_i とy_i (i = 1, ... , 7) で [S] とv のデータを表す。また、β₁ = V_max 、β₂ = K_M とする。残差

${\displaystyle r_{i}=y_{i}-{\frac {\beta _{1}x_{i}}{\beta _{2}+x_{i}}}\quad (i=1,\dots ,7)}$

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

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

${\displaystyle {\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 _{1}x_{i}}{\left(\beta _{2}+x_{i}\right)^{2}}}\end{pmatrix}}.}$

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

収束性

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

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

${\displaystyle {\begin{aligned}r_{1}(\beta )&=\beta +1\\r_{2}(\beta )&=\lambda \beta ^{2}+\beta -1.\end{aligned}}}$

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

ニュートン法からの導出

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

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

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

である。ここでg はS の勾配ベクトル、H はS のヘッシアンである。${\displaystyle S=\sum _{i=1}^{m}r_{i}^{2}}$であるから、勾配g は次で与えられる：

${\displaystyle {\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 をβ で微分することで計算される：

${\displaystyle 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項）を無視することでガウス・ニュートン法を得る。つまり、ヘッシアンは

${\displaystyle 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}}$

と近似される。ここで

${\displaystyle J_{ij}={\frac {\partial r_{i}}{\partial \beta _{j}}}}$

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

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

${\displaystyle {\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階微分項を無視するという近似、すなわち

${\displaystyle \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. r_i は十分小さい。少なくとも最小値付近。
2. 関数の非線形性は穏やかであり、${\displaystyle {\partial ^{2}r_{i}}/{\partial \beta _{j}\partial \beta _{k}}}$ が比較的小さくなる。

改善バージョン

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

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

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

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

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

最適な縮小因子αが 0 に近いような場合、発散を回避する別の方法はレーベンバーグ・マルカート法（信頼領域法）を使うことである^[2]。増分ベクトルが最急降下方向に向くように正規方程式は修正される。

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

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

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

関連するアルゴリズム

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

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

脚注

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