ガウス・ニュートン法
ガウス・ニュートン法(ガウス・ニュートンほう、Gauss–Newton method)は、非線形最小二乗法を解く手法の一つである。これは関数の最大・最小値を見出すニュートン法の修正とみなすことができる。ニュートン法とは違い、ガウス・ニュートン法は二乗和の最小化にしか用いることができないが、計算するのが困難な2階微分が不要という長所がある。
非線形最小二乗法は非線形回帰などで、観測データを良く表すようにモデルのパラメータを調整するために必要となる。
この手法の名称はカール・フリードリヒ・ガウスとアイザック・ニュートンにちなむ。
データフィッティングにおいて、与えられたモデル関数 y = f (x , β) がm 個のデータ点 {(xi , yi ); i = 1, ... , m } に最もよくフィットするようなn (≤ m )個のパラメータβ = (β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
の最小値を反復計算で求める。初期推測値β(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)}).
ここでは例として、ガウス・ニュートン法を使ってデータとモデルによる予測値の間の残差平方和を最小化し、データにモデルをフィットさせる。
生物実験において、酵素媒介反応における基質濃度[S] と反応率v の関係として次表のようなデータが得られたとする(右図の赤点)。
これらのデータに対し、次の形のモデル曲線(ミカエリス・メンテン式)のパラメータVmax とKM を、最小二乗の意味で最もよくフィットするように決定したい:
v=\frac{V_\mathrm{max}[\mathrm{S}]}{K_\mathrm{M}+[\mathrm{S}]}
xi とyi (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}\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}.
初期推定値としてβ1 = 0.9、β2 = 0.2から始め、ガウス・ニュートン法による5回の反復計算を行うと、最適値 \hat\beta_1=0.362, \hat\beta_2=0.556 が得られる。残差平方和は5回の反復計算で初期の1.445から0.00784まで減少する。右図はこれらの最適パラメータを用いたモデルで決まる曲線と、実験データとの比較を示す。
増分ΔがS の減少方向を向いていることは証明されている。もしこのアルゴリズムが収束すれば、その極限はS の停留点である。しかし収束については、ニュートン法では保証されている局所収束さえも保証されていない。
ガウス・ニュートン法の収束の速さは2次である。もし初期推測値が最小値から遠いか、または行列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 なら、この方法はもはや局所的にも収束しない。
後に示すように、ガウス・ニュートン法は近似関数の最適化に用いられるニュートン法から与えられる。その結果、ガウス・ニュートン法の収束の速さはほとんど2次である。
パラメータβを持つ関数S の最小化をするとき、ニュートン法による漸化式は
\boldsymbol\beta^{(s+1)} = \boldsymbol\beta^{(s)} - H^{-1} \boldsymbol{g} \,
である。ここでg はS の勾配ベクトル、H はS のヘッシアンである。 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つの条件の下であり、これらが成り立つ場合には収束が期待される:
ri は十分小さい。少なくとも最小値付近。
関数の非線形性は穏やかであり、{\partial^2 r_i}/{\partial \beta_j \partial \beta_k} が比較的小さくなる。
ガウス・ニュートン法は、初期推定値が真の解から大きく離れていたり、モデル関数の非線形性が大きい場合には安定性が悪い。また、残差平方和S は反復ごとに必ずしも減少するわけではない。そのため、実用上は安定化が必要である。
S は反復ごとに必ずしも減少するわけではないが、増分ベクトルΔはS が減少する方向を向いているから、S (βs) が停留点にない限り、任意の十分に小さなα> 0 に対して S (βs + αΔ) < S (βs) が成り立つ。したがって、発散したときに、更新方程式に縮小因子αを導入して
\boldsymbol \beta^{s+1} = \boldsymbol \beta^s+\alpha\boldsymbol{\Delta}
とすることが解決法の一つとなる。
言い換えれば、増分ベクトルΔは目的関数S の下り方向を指してはいるが長すぎるので、その道のほんの一部を行くことで、S を減少させようというアイディアである。縮小因子αの最適値は直線探索で見つけることができる。つまり、直接探索法を(通常 0 < α ≤ 1 の区間で)用いて、S を最小化する値を探すことでαの大きさは決められる。
最適な縮小因子αが 0 に近いような場合、発散を回避する別の方法はレーベンバーグ・マルカート法(信頼領域法)を使うことである。増分ベクトルが最急降下方向に向くように正規方程式は修正される。
(J^\mathrm{T} J+\lambda D)\Delta = J^\mathrm{T} \boldsymbol{r},
ここでD は正定値対角行列である。λが 0 から増大するにつれて増分ベクトルΔは長さが単調に減少し、かつ方向は最急降下方向に近づくため、λを十分大きくすれば必ずより小さいS の値を見出せることが保証されている。
いわゆるマルカートパラメータλは直線探索により最適化されるが、λが変わるたびに毎回シフトベクトルの再計算をしなければならないため非効率的である。より効率的な方法は、発散が起きた時にλをS が減少するまで増加させる。そしてその値を1回の反復から次まで維持する、しかしλを 0 に設定することができるときにもしカットオフ値に届くまで可能なら減少させる。このときS の最小値は標準のガウス・ニュートン法の最小化になる。
DFP法やBFGS法のような準ニュートン法では、ヘッシアン{\partial^2 S}/{\partial \beta_j \partial\beta_k}の推定は1階微分{\partial r_i}/{\partial\beta_j}のみを用いて数値的になされる。したがってn 回の反復計算による修正の後、この方法はパフォーマンスにおいてニュートン法を近似する。ガウス・ニュートン法やレーベンバーグ・マルカート法などは非線形最小二乗問題にのみ適用できるのに対して、準ニュートン法は一般的な実数値関数を最小化できることに注意する。
1階微分のみを使って最小化問題を解く他の方法は、最急降下法である。しかし、その方法は近似に2階微分を考慮していないので、多くの関数に対して、特にパラメータが強い相互作用を持っている場合には、計算効率が非常に悪い。