ニュートン法の導出

この記事について

非線形最小二乗法の数値解法であるニュートン法とその改良型のLM法までを解説していきます．

ニュートン法の導出

$x\in \mathbb{R}^n$を引数とする$F(x):x\rightarrow F(x)\in \mathbb{R}^m$という関数を考えます．ちなみに、$x$と$F(x)$の元（要素）を$x_i|i= \{1,\dots,n\} $と$f_i(x)|i=\{1,\dots,m\}$とすると以下のような列ベクトルとして表せます．

\begin{eqnarray}
x &=&
\begin{bmatrix}
x_1 \\
\vdots \\
x_n
\end{bmatrix} \\
F(x) &=&
\begin{bmatrix}
f_1(x) \\
\vdots \\
f_m(x)
\end{bmatrix} \\
\end{eqnarray}

最小二乗法では以下のように定義される残差平方和$S(x)$を最小とするな$x$を求めます．

S(x)=\frac{1}{2}||F(x)||^2

$S(x)$をある$x_k$周りにテイラー展開します．すると、以下が成り立ちます．

S(x_{k+1}) = S(x_k)+ \nabla S(x_k) (x_{k+1}-x_k)+\frac{1}{2}(x_{k+1}-x_k)^T\Delta S(x_k)(x_{k+1}-x_k)+\cdots

$S(x)\rightarrow \min$となる時、$\nabla S(x_{k+1})=0$が成り立ちます．テイラー展開した式を$x_{k+1}$で偏微分することで以下の式が求まります．

\nabla ^T S(x_k)+\frac{1}{2}(\Delta S(x_k) ^T+\Delta S(x_k))(x_{k+1}-x_k)+\cdots = 0

上の式の導出にあたって$\nabla (xAx)=2(A+A^T)x$を使いました．更に、$\nabla S(x_{k+1})$は対称行列であることから以下のように式は簡略化されます．

\nabla ^T S(x_k)+\Delta S(x_k)(x_{k+1}-x_k)+\cdots = 0

さらに、3項以降を無視して

\nabla ^T S(x_k)+\Delta S(x_k)(x_{k+1}-x_k) \approx 0

これを$x_{k+1}$について解くと

x_{k+1} \approx x_k - (\Delta S(x_k))^{-1} \nabla ^T S(x_k)

ガウス・ニュートン法の導出

このままでも使える式なのですが、$S(x_k)$があるとフロベニウス和を計算しなればなりません．面倒なので$S(x)$を使わない形にできないかを考えていきます．手始めに、$\nabla S(x)$を簡略化することを考えます．

\begin{eqnarray}
\nabla ^T S(x) &=& \nabla ^T (\frac{1}{2}\sum _{i=1}^m f_i(x)^2) \\
&=& 
\begin{bmatrix}
\frac{\partial S(x)}{\partial x_1} & \cdots & \frac{\partial S(x)}{\partial x_n}
\end{bmatrix}
\end{eqnarray}

$\nabla ^T S(x)$の第$j$列成分$\frac{\partial S(x)}{\partial x_j}$は以下の式で表すことができます．

\frac{\partial S(x)}{\partial x_j} = \sum _{i=1}^n \frac{\partial f_i(x)}{\partial x_j} f_i(x)

このことより、$F(x)$の定義と$\nabla F(x)$の定義を考慮すると

\nabla ^T S(x)= \nabla F(x)F(x)=F(x)^T J_{F(x)}

となります．ここに$J_{F(x)}$は$F(x)$のヤコビアンです．

次に$\Delta S(x)$を簡略化します．$\Delta S(x)$の$j$行$l$列の成分を$\Delta S(x)_(j,k)$とすると

\begin{eqnarray}
\Delta S(x)_{(j,l)} &=& \frac{\partial}{\partial x_l}(\frac{\partial S(x)}{\partial x_j} ) \\
&=& \frac{\partial}{\partial x_l}(\sum _{i=1}^n \frac{\partial f_i(x)}{\partial x_j} f_i(x)) \\
&=& \sum _{i=1}^n (\frac{\partial f_i(x)}{\partial x_j} \frac{\partial f_i(x)}{\partial x_l} +  \frac{\partial  ^2f_i(x)}{\partial x_j \partial x_l}f_i(x)) 
\end{eqnarray}

$S(x)\rightarrow \min$の時、$F(x)\rightarrow 0$なので、$f_i\rightarrow 0$が成り立つ．したがって、

\Delta S(x)_{(j,l)} \approx \sum _{i=1}^n \frac{\partial f_i(x)}{\partial x_j} \frac{\partial f_i(x)}{\partial x_l} 

行列の形式に直すと、

\Delta S(x) \approx J_{F(x)}^TJ_{F(x)} 

さて、得られた近似式を

x_{k+1} \approx x_k - (J_{F(x)}^T J_{F(x)})^{-1} J_{F(x)}^T F(x)

さらに$J^+ {F(x)}$を$J{F(x)}$のMP逆行列（疑似逆行列）とすると

x_{k+1} \approx x_k - J_{F(x)}^+ F(x)

ちなみに、$n=m$のときに限り、$J^+ _{F(x)}=J^{-1} _{F(x)}$となります．この時、この解法はニュートン・ラフソン法と呼ばれます．ニュートン・ラフソン法の式は以下で定義されます．

x_{k+1} \approx x_k - J_{F(x)}^{-1} F(x)