非線形最小二乗法とパラメータ推定 | 数値解析 I

1 残差写像と Gauss–Newton 法

定義 1.1. $m,n\in\NN$ とし、 $U\subseteq\R^n$ を開集合とする。

$r\colon U\to\R^m$ を $C^1$ 級の写像とし、全微分 $Dr(x)$ を標準基底に関する $m\times n$ 行列と同一視して $J(x):=Dr(x)$ と書く。 $\phi\colon U\to\R$ を $\phi(x):=\frac12\lVert r(x)\rVert_2^2$ で定める。 $\phi$ の局所最小点を求める問題を、残差写像 (residual map) $r$ に関する非線形最小二乗問題 (nonlinear least squares problem) という。
$h\colon U\times\R\to\R$ を、各 $t\in\R$ について $x\mapsto h(x,t)$ が $C^1$ 級である関数とする。観測時刻 $t_1,\dots,t_m\in\R$ 、観測値 $y\in\R^m$ 、正の実数 $w_1,\dots,w_m$ に対して、 $r_i(x):=\sqrt{w_i}\,\bigl(h(x,t_i)-y_i\bigr)$ （ $1\le i\le m$ ）で定まる残差写像 $r$ の非線形最小二乗問題を、モデル $h$ のデータ $(t_i,y_i)_{i=1}^m$ への重み $w$ による当てはめ (curve fitting) という。このとき $\phi(x)=\frac12\sum_{i=1}^mw_i\bigl(h(x,t_i)-y_i\bigr)^2$ である。

補題 1.2. $m$ 、 $n$ 、 $U$ 、 $r$ 、 $J$ 、 $\phi$ を定義 1.1のとおりとする。

$\phi$ は $C^1$ 級であり、任意の $x\in U$ について $\nabla\phi(x)=J(x)^{\mathsf T}r(x)$ である。
$r$ が $C^2$ 級ならば $\phi$ は $C^2$ 級であり、任意の $x\in U$ について $\nabla^2\phi(x)=J(x)^{\mathsf T}J(x)+\sum_{i=1}^mr_i(x)\nabla^2r_i(x)$ である。

証明. 演習とする（問題 8.1）。▨

定義 1.3. $m$ 、 $n$ 、 $U$ 、 $r$ 、 $J$ 、 $\phi$ を定義 1.1のとおりとし、 $D_{\mathrm{GN}}:=\{x\in U\mid J(x)\text{ の列は一次独立}\}$ と置く。 $D_{\mathrm{GN}}\ne\emptyset$ ならば $m\ge n$ である。 $x\in D_{\mathrm{GN}}$ とする。§E20.6 命題 1.2 (2)により、 $J(x)p=-r(x)$ の最小二乗解、すなわち $p\mapsto\lVert r(x)+J(x)p\rVert_2$ の最小値を与える $p\in\R^n$ はただ一つ存在し、§E20.9 命題 6.1 (3)により

s(x):=-\bigl(J(x)^{\mathsf T}J(x)\bigr)^{-1}J(x)^{\mathsf T}r(x)=-J(x)^\dagger r(x)

である。 $s(x)$ を $x$ における Gauss–Newton 修正量 (Gauss-Newton step) といい、 $G(x):=x+s(x)$ で定まる $G\colon D_{\mathrm{GN}}\to\R^n$ を Gauss–Newton 反復写像 (Gauss-Newton iteration map) という。 $x_0\in U$ とし、 $x_k\in D_{\mathrm{GN}}$ であるとき $x_{k+1}:=G(x_k)$ と置く。この反復を Gauss–Newton 法 (Gauss-Newton method) といい、すべての $k\in\N$ で $x_k\in D_{\mathrm{GN}}$ であるとき、Gauss–Newton 法の反復列 $(x_k)_{k\in\N}$ が定まるという。 $J(x)^{\mathsf T}J(x)$ は正則であるから、補題 1.2 (1)により、 $x\in D_{\mathrm{GN}}$ について $G(x)=x$ であることと $\nabla\phi(x)=0$ であることは同値である。 $m=n$ ならば $J(x)^\dagger=J(x)^{-1}$ であり、 $G$ は§E20.23 定義 1.1の $N_r$ に一致する。

注意 1.4. $x\in D_{\mathrm{GN}}$ とする。 $J(x)$ と右辺 $-r(x)$ の Householder QR 法を厳密算術で実行し、その出力を $(R,c)$ とすると、§E20.6 命題 3.6により $s(x)$ は $Rs=(c_1,\dots,c_n)$ のただ一つの解であり、§E20.6 補題 3.5により $\min_p\lVert r(x)+J(x)p\rVert_2=\lVert(c_{n+1},\dots,c_m)\rVert_2$ である。この計算は $J(x)^{\mathsf T}J(x)$ を作らない。 $J(x)^{\mathsf T}J(x)$ を係数行列とする一次方程式を解く場合、§E20.6 命題 2.3により係数行列の条件数は $\kappa_2(J(x))^2$ である。

$r$ が計算グラフの定める写像であるとき、§E20.22 命題 4.1 (3)により、 $J(x)$ の列は方向 $e_1,\dots,e_n$ の前進モードで、行は重み $e_1,\dots,e_m$ の逆モードで得られる。§E20.22 定理 3.2により、重み $r(x)$ の一回の逆モードは $J(x)^{\mathsf T}r(x)=\nabla\phi(x)$ を与える。

$r$ が定義 1.1 (2)の重み $w$ による当てはめの残差写像であるとき、 $J(x)$ の第 $i$ 行は $\sqrt{w_i}\,D_xh(x,t_i)$ である。ここで $D_xh(x,t_i)\in\R^{1\times n}$ は $x\mapsto h(x,t_i)$ の全微分である。第 $i$ 行が $D_xh(x,t_i)$ の行列を $A$ 、 $b:=\bigl(y_i-h(x,t_i)\bigr)_{i=1}^m$ 、 $C:=\operatorname{diag}(\sqrt{w_1},\dots,\sqrt{w_m})$ とすると $CA=J(x)$ 、 $Cb=-r(x)$ であり、§E20.6 命題 5.2 (2)により、 $J(x)p=-r(x)$ の最小二乗解は $Ap=b$ の重み $\operatorname{diag}(w_1,\dots,w_m)$ に関する重み付き最小二乗解である。

2 零残差での局所収束

補題 2.1. $m$ 、 $n$ 、 $U$ 、 $r$ 、 $J$ を定義 1.1のとおりとし、 $D_{\mathrm{GN}}$ 、 $G$ を定義 1.3のとおりとする。

$D_{\mathrm{GN}}$ は開集合であり、 $x\mapsto J(x)^\dagger$ は $D_{\mathrm{GN}}$ 上で連続である。 $r$ が $C^2$ 級ならば、 $x\mapsto J(x)^\dagger$ と $G$ は $D_{\mathrm{GN}}$ 上で $C^1$ 級である。
$x_*\in D_{\mathrm{GN}}$ とし、 $\sigma_*:=\sigma_n(J(x_*))$ と置く。 $x\in U$ が $\lVert J(x)-J(x_*)\rVert_2\le\sigma_*/2$ を満たすならば、 $x\in D_{\mathrm{GN}}$ であり、 $J(x)^\dagger J(x)=I_n$ 、 $\lVert J(x)^\dagger\rVert_2\le2/\sigma_*$ が成り立つ。

証明.(1)を示す。 $x\in D_{\mathrm{GN}}$ をとる。§E20.6 補題 2.1 (2)により $\sigma_n(J(x))>0$ である。§E20.6 補題 2.1 (1)により $\lVert J(y)-J(x)\rVert_2\le\lVert J(y)-J(x)\rVert_F$ であり、 $J$ は連続であるから、 $x$ の開近傍 $W\subseteq U$ が存在して、 $y\in W$ ならば $\lVert J(y)-J(x)\rVert_2<\sigma_n(J(x))$ である。§E20.6 補題 2.1 (3)により、 $y\in W$ ならば $J(y)$ の列は一次独立であり、 $W\subseteq D_{\mathrm{GN}}$ である。したがって $D_{\mathrm{GN}}$ は開集合である。 $y\in D_{\mathrm{GN}}$ について $M(y):=J(y)^{\mathsf T}J(y)$ は§E20.6 命題 1.2 (2)により正則であり、余因子行列による逆行列の表示から、 $M(y)^{-1}$ の各成分は $M(y)$ の成分の多項式を $\det M(y)\ne0$ で割ったものである。 $M(y)$ の成分は $J(y)$ の成分の多項式であるから、 $J(y)^\dagger=M(y)^{-1}J(y)^{\mathsf T}$ の各成分は、 $J(y)$ の成分の有理式で分母が $D_{\mathrm{GN}}$ 上で $0$ にならないものである。 $J$ は連続であるから $x\mapsto J(x)^\dagger$ は連続である。 $r$ が $C^2$ 級ならば $J$ は $C^1$ 級であるから、 $x\mapsto J(x)^\dagger$ は $C^1$ 級であり、 $G(x)=x-J(x)^\dagger r(x)$ も $C^1$ 級である。

(2)を示す。§E20.6 補題 2.1 (2)により $\sigma_*>0$ であり、 $\lVert J(x)-J(x_*)\rVert_2\le\sigma_*/2<\sigma_*$ であるから、§E20.6 補題 2.1 (3)により $J(x)$ の列は一次独立であって $\sigma_n(J(x))\ge\sigma_*/2$ である。 $J(x)^\dagger J(x)=(J(x)^{\mathsf T}J(x))^{-1}J(x)^{\mathsf T}J(x)=I_n$ であり、§E20.9 命題 6.1 (3)により $\lVert J(x)^\dagger\rVert_2=1/\sigma_n(J(x))\le2/\sigma_*$ である。▨

定理 2.2. $m$ 、 $n$ 、 $U$ 、 $r$ 、 $J$ を定義 1.1のとおりとし、 $D_{\mathrm{GN}}$ 、 $G$ を定義 1.3のとおりとする。 $\R^n$ と $\R^m$ には Euclid ノルムを、行列にはそれに関する作用素ノルムを用いる。 $x_*\in U$ は $r(x_*)=0$ を満たし、 $J(x_*)$ の列は一次独立であるとし、 $\sigma_*:=\sigma_n(J(x_*))$ と置く。実数 $R>0$ 、 $\gamma\ge0$ について、閉球 $\overline B(x_*,R):=\{x\in\R^n\mid\lVert x-x_*\rVert_2\le R\}$ が $U$ に含まれ、任意の $x,y\in\overline B(x_*,R)$ に対して $\lVert J(x)-J(y)\rVert_2\le\gamma\lVert x-y\rVert_2$ が成り立つとする。 $\gamma>0$ ならば $\rho:=\min\{R,\sigma_*/(2\gamma)\}$ 、 $\gamma=0$ ならば $\rho:=R$ と置く。

任意の $x\in\overline B(x_*,\rho)$ について $x\in D_{\mathrm{GN}}$ であり、 $\lVert G(x)-x_*\rVert_2\le\frac{\gamma}{\sigma_*}\lVert x-x_*\rVert_2^2\le\frac12\lVert x-x_*\rVert_2$ が成り立つ。
任意の $x_0\in\overline B(x_*,\rho)$ について Gauss–Newton 法の反復列 $(x_k)$ が定まり、すべての $k\in\N$ で $x_k\in\overline B(x_*,\rho)$ であって、 $e_k:=x_k-x_*$ は $\lVert e_{k+1}\rVert_2\le\frac{\gamma}{\sigma_*}\lVert e_k\rVert_2^2,\qquad\lVert e_{k+1}\rVert_2\le\frac12\lVert e_k\rVert_2$ を満たす。 $(x_k)$ は $x_*$ へ少なくとも二次で収束する。 $\gamma=0$ ならば $x_1=x_*$ であり、ある $k$ で $x_k=x_*$ ならば、 $j\ge k$ を満たすすべての $j$ で $x_j=x_*$ である。

証明.(1)を示す。 $x\in\overline B(x_*,\rho)$ をとり、 $e:=x-x_*$ と置く。Lipschitz 条件と $\gamma\rho\le\sigma_*/2$ により $\lVert J(x)-J(x_*)\rVert_2\le\gamma\lVert e\rVert_2\le\sigma_*/2$ であるから、補題 2.1 (2)により $x\in D_{\mathrm{GN}}$ 、 $J(x)^\dagger J(x)=I_n$ 、 $\lVert J(x)^\dagger\rVert_2\le2/\sigma_*$ である。 $r(x_*)=0$ であるから

G(x)-x_*=e-J(x)^\dagger r(x)=J(x)^\dagger\bigl(J(x)e-r(x)\bigr)=J(x)^\dagger\bigl(r(x_*)-r(x)-J(x)(x_*-x)\bigr)

である。 $\overline B(x_*,\rho)$ は凸であるから、 $t\in[0,1]$ について $x+t(x_*-x)\in\overline B(x_*,\rho)\subseteq U$ であり、Lipschitz 条件から $\lVert J(x+t(x_*-x))-J(x)\rVert_2\le\gamma t\lVert e\rVert_2$ である。§E20.23 補題 2.2を $F=r$ 、 $a=x$ 、 $h=x_*-x$ に適用すると、右辺の括弧内のベクトルのノルムは $\frac\gamma2\lVert e\rVert_2^2$ 以下であり、

\lVert G(x)-x_*\rVert_2\le\frac2{\sigma_*}\cdot\frac\gamma2\lVert e\rVert_2^2=\frac\gamma{\sigma_*}\lVert e\rVert_2^2

を得る。 $\gamma\lVert e\rVert_2\le\gamma\rho\le\sigma_*/2$ であるから、右辺は $\frac12\lVert e\rVert_2$ 以下である。

(2)を示す。 $x_k\in\overline B(x_*,\rho)$ ならば、(1)により $x_k\in D_{\mathrm{GN}}$ であって $x_{k+1}=G(x_k)$ が定まり、二つの不等式が成り立ち、 $\lVert e_{k+1}\rVert_2\le\rho/2$ である。 $x_0\in\overline B(x_*,\rho)$ から、 $k$ に関する帰納法により、反復列が定まり、すべての $k$ で $x_k\in\overline B(x_*,\rho)$ と二つの不等式が成り立つ。 $\lVert e_k\rVert_2\le2^{-k}\lVert e_0\rVert_2$ であるから $x_k\to x_*$ であり、第一の不等式は§E20.4 定義 1.5の不等式を $C=\gamma/\sigma_*$ 、 $k_0=0$ として与える。 $\gamma=0$ ならば第一の不等式から $\lVert e_1\rVert_2\le0$ である。 $e_k=0$ ならば第一の不等式から $e_{k+1}=0$ であり、 $j$ に関する帰納法により $j\ge k$ について $e_j=0$ である。▨

注意 2.3. $m=n$ のとき、§E20.9 命題 6.1 (3)により $\lVert J(x_*)^{-1}\rVert_2=1/\sigma_*$ であるから、定理 2.2の $\rho$ と定数 $\gamma/\sigma_*$ は、 $F=r$ についての§E20.23 定義 2.1の $\rho$ と§E20.23 定理 2.5の定数 $\beta\gamma$ に一致する。

3 非零残差と反復写像の微分

命題 3.1. $m$ 、 $n$ 、 $U$ 、 $r$ 、 $J$ 、 $\phi$ を定義 1.1のとおりとし、 $D_{\mathrm{GN}}$ 、 $G$ を定義 1.3のとおりとする。 $r$ は $C^2$ 級であるとし、 $x_*\in D_{\mathrm{GN}}$ は $\nabla\phi(x_*)=0$ を満たすとする。 $J_*:=J(x_*)$ 、 $S_*:=\sum_{i=1}^mr_i(x_*)\nabla^2r_i(x_*)$ と置く。このとき $G$ は $x_*$ を含む開集合 $D_{\mathrm{GN}}$ 上で $C^1$ 級であり、 $G(x_*)=x_*$ であって

DG(x_*)=-(J_*^{\mathsf T}J_*)^{-1}S_*

が成り立つ。特に $r(x_*)=0$ ならば $DG(x_*)=0$ である。

証明.補題 2.1 (1)により $D_{\mathrm{GN}}$ は開集合であり、 $G$ はその上で $C^1$ 級である。定義 1.3により $G(x_*)=x_*$ である。 $x\in D_{\mathrm{GN}}$ について $M(x):=J(x)^{\mathsf T}J(x)$ と置くと、 $G$ の定義と補題 1.2 (1)により

M(x)\bigl(G(x)-x\bigr)=-J(x)^{\mathsf T}r(x)=-\nabla\phi(x)

である。 $r$ は $C^2$ 級であるから $M$ と $\nabla\phi$ は $C^1$ 級であり、両辺を $x_*$ で方向 $h\in\R^n$ に微分すると、 $G(x_*)-x_*=0$ と補題 1.2 (2)により

M(x_*)\bigl(DG(x_*)h-h\bigr)=-\nabla^2\phi(x_*)h=-M(x_*)h-S_*h

である。 $M(x_*)$ は正則であるから $DG(x_*)h=-M(x_*)^{-1}S_*h$ を得る。 $r(x_*)=0$ ならば $S_*=0$ である。▨

定理 3.2. $m$ 、 $n$ 、 $U$ 、 $r$ 、 $\phi$ 、 $D_{\mathrm{GN}}$ 、 $G$ 、 $x_*$ を命題 3.1のとおりとする。 $\R^n$ にノルム $\lVert\cdot\rVert$ を一つ固定し、 $n$ 次正方行列にはそれに関する作用素ノルムを用いる。 $\lVert DG(x_*)\rVert<1$ とし、実数 $q$ が $\lVert DG(x_*)\rVert<q<1$ を満たすとする。

実数 $\delta>0$ が存在して、 $\lVert x_0-x_*\rVert\le\delta$ を満たす任意の $x_0\in\R^n$ について、Gauss–Newton 法の反復列 $(x_k)$ が定まり、すべての $k\in\N$ で $\lVert x_k-x_*\rVert\le\delta$ かつ $\lVert x_{k+1}-x_*\rVert\le q\lVert x_k-x_*\rVert$ が成り立ち、 $x_k\to x_*$ である。
(1)の反復列がすべての $k$ で $x_k\ne x_*$ を満たすならば、 $\limsup_{k\to\infty}\lVert x_{k+1}-x_*\rVert/\lVert x_k-x_*\rVert\le\lVert DG(x_*)\rVert$ である。

証明.

主張 3.2.1. $\lVert DG(x_*)\rVert<q'<1$ を満たす任意の実数 $q'$ について、実数 $\delta'>0$ が存在して、 $B':=\{x\in\R^n\mid\lVert x-x_*\rVert\le\delta'\}$ は $D_{\mathrm{GN}}$ に含まれ、任意の $x\in B'$ について $\lVert G(x)-x_*\rVert\le q'\lVert x-x_*\rVert$ が成り立つ。

証明.命題 3.1により $D_{\mathrm{GN}}$ は $x_*$ を含む開集合であり、 $DG$ は $D_{\mathrm{GN}}$ 上で連続である。作用素ノルムは $n$ 次正方行列の空間のノルムであり、有限次元のノルムはすべて同値であるから、 $x\mapsto\lVert DG(x)-DG(x_*)\rVert$ は $x_*$ で連続である。この連続性と§E20.2 補題 1.2 (2)により、 $\delta'>0$ が存在して、 $B'\subseteq D_{\mathrm{GN}}$ であり、任意の $x\in B'$ について $\lVert DG(x)\rVert\le\lVert DG(x_*)\rVert+\lVert DG(x)-DG(x_*)\rVert\le q'$ である。 $x\in B'$ をとる。 $B'$ は凸であるから、 $t\in[0,1]$ について $x_*+t(x-x_*)\in B'$ であり、§E20.2 命題 4.1を開集合 $D_{\mathrm{GN}}$ 上の $G$ 、点 $x_*$ 、 $\Delta x=x-x_*$ 、 $L=q'$ に適用すると $\lVert G(x)-G(x_*)\rVert\le q'\lVert x-x_*\rVert$ である。 $G(x_*)=x_*$ である。▨

(1)を示す。主張 3.2.1を $q'=q$ に適用して $\delta:=\delta'$ と置く。 $\lVert x_k-x_*\rVert\le\delta$ ならば、 $x_k\in D_{\mathrm{GN}}$ であって $x_{k+1}=G(x_k)$ が定まり、 $\lVert x_{k+1}-x_*\rVert\le q\lVert x_k-x_*\rVert\le\delta$ である。 $k$ に関する帰納法により、反復列が定まり、すべての $k$ で二つの不等式が成り立つ。 $\lVert x_k-x_*\rVert\le q^k\delta$ であるから $x_k\to x_*$ である。

(2)を示す。 $\lVert DG(x_*)\rVert<q'<1$ を満たす $q'$ をとり、主張 3.2.1の $\delta'$ をとる。 $x_k\to x_*$ であるから、 $K\in\N$ が存在して、 $k\ge K$ ならば $\lVert x_k-x_*\rVert\le\delta'$ であり、 $\lVert x_{k+1}-x_*\rVert/\lVert x_k-x_*\rVert\le q'$ である。したがって上極限は $q'$ 以下であり、 $q'$ は $\lVert DG(x_*)\rVert$ より大きい任意の実数であるから、上極限は $\lVert DG(x_*)\rVert$ 以下である。▨

例 3.3. $\lambda\in\R$ とし、 $r\colon\R\to\R^2$ を $r(x):=(x+1,\ \lambda x^2+x-1)$ で定める。 $J(x)=(1,\ 2\lambda x+1)^{\mathsf T}$ は $0$ でないから $D_{\mathrm{GN}}=\R$ であり、

G(x)=x-\frac{(x+1)+(\lambda x^2+x-1)(2\lambda x+1)}{1+(2\lambda x+1)^2}

である。 $r(0)=(1,-1)$ 、 $J(0)=(1,1)^{\mathsf T}$ から $\nabla\phi(0)=0$ 、 $J(0)^{\mathsf T}J(0)=2$ であり、残差のノルム $\lVert r(0)\rVert_2=\sqrt2$ は $\lambda$ によらない。 $\nabla^2r_1=0$ 、 $\nabla^2r_2=2\lambda$ から $S_*=-2\lambda$ であり、命題 3.1により $G'(0)=\lambda$ 、補題 1.2 (2)により $\phi''(0)=2-2\lambda$ である。

$\lvert\lambda\rvert<1$ ならば、定理 3.2 (1)により、実数 $\delta>0$ が存在して、 $\lvert x_0\rvert\le\delta$ を満たす任意の初期値 $x_0$ からの反復列は $0$ に収束する。そのような反復列 $(x_k)$ がすべての $k\in\N$ で $x_k\ne0$ を満たすならば、定理 3.2 (2)により $\lvert x_{k+1}\rvert/\lvert x_k\rvert$ の上極限は $\lvert\lambda\rvert$ 以下である。 $\lambda=\frac12$ 、 $x_0=\frac1{10}$ では $x_1=\frac{211}{4420}\approx4.774\times10^{-2}$ 、 $x_2\approx2.333\times10^{-2}$ 、 $x_3\approx1.153\times10^{-2}$ であり、 $x_{k+1}/x_k$ （ $k=0,1,2$ ）は $0.477$ 、 $0.489$ 、 $0.494$ である。
$\lambda=0$ ならば $G(x)=x-\frac{(x+1)+(x-1)}2=0$ であるから、任意の初期値について $x_1=0$ である。
$\lvert\lambda\rvert>1$ とする。 $D_{\mathrm{GN}}=\R$ であるから、任意の初期値 $x_0\in\R$ について Gauss–Newton 法の反復列は $x_0$ から始まる $G$ の不動点反復である。 $G$ は $0$ で微分可能であり、 $G(0)=0$ 、 $\lvert G'(0)\rvert=\lvert\lambda\rvert>1$ であるから、§E20.4 命題 3.4を $D=\R$ 、 $g=G$ 、 $x^*=0$ に適用すると、 $0$ に収束する反復列については、 $K\in\N$ が存在して、 $k\ge K$ を満たすすべての $k$ で $x_k=0$ である。 $\lambda<-1$ ならば $\phi''(0)=2-2\lambda>0$ であり、 $0$ は $\phi$ の狭義の局所最小点である。 $\lambda=-2$ 、 $x_0=\frac1{10}$ では $x_1=-\frac{103}{340}\approx-0.3029$ である。

4 Levenberg–Marquardt 法

定義 4.1. $m$ 、 $n$ 、 $U$ 、 $r$ 、 $J$ 、 $\phi$ を定義 1.1のとおりとし、 $x\in U$ 、実数 $\mu>0$ をとる。§E20.25 命題 2.2 (1)を $A=J(x)$ 、 $L=I_n$ 、 $c=-r(x)$ 、 $\lambda=\sqrt\mu$ に適用すると、 $\lVert r(x)+J(x)p\rVert_2^2+\mu\lVert p\rVert_2^2$ の最小値を与える $p\in\R^n$ はただ一つ存在し、それは

\bigl(J(x)^{\mathsf T}J(x)+\mu I_n\bigr)p=-J(x)^{\mathsf T}r(x)

のただ一つの解であって、 $\begin{pmatrix}J(x)\\\sqrt\mu\,I_n\end{pmatrix}p=\begin{pmatrix}-r(x)\\0\end{pmatrix}$ のただ一つの最小二乗解である。この $p$ を $p_\mu(x)$ と書き、 $x$ における減衰係数 $\mu$ の Levenberg–Marquardt 修正量 (Levenberg-Marquardt step) といい、 $\mu$ を減衰係数 (damping parameter) という。§E20.25 定義 2.3の行列 $R_\lambda$ を $A=J(x)$ について作ると $p_\mu(x)=R_{\sqrt\mu}\bigl(-r(x)\bigr)$ であり、 $p_\mu(x)$ はデータ $-r(x)$ の、正則化係数 $\sqrt\mu$ の Tikhonov 正則化解である。拡大行列の列は下の $n$ 行 $\sqrt\mu\,I_n$ により一次独立であるから、§E20.6 命題 3.6により、 $p_\mu(x)$ はこの拡大行列と右辺 $(-r(x),0)$ の Householder QR 法で計算することができ、 $J(x)^{\mathsf T}J(x)+\mu I_n$ を作る必要はない。

命題 4.2. $m$ 、 $n$ 、 $U$ 、 $r$ 、 $J$ 、 $\phi$ を定義 1.1のとおりとし、 $D_{\mathrm{GN}}$ 、 $s$ を定義 1.3のとおりとする。 $x\in U$ 、実数 $\mu>0$ をとり、 $g:=\nabla\phi(x)$ 、 $p:=p_\mu(x)$ と置く。

$g\ne0$ ならば $p\ne0$ であり、 $g^{\mathsf T}p=-\lVert J(x)p\rVert_2^2-\mu\lVert p\rVert_2^2<0,\qquad\phi(x)-\frac12\lVert r(x)+J(x)p\rVert_2^2=\frac12\lVert J(x)p\rVert_2^2+\mu\lVert p\rVert_2^2>0$ が成り立つ。
$\lVert p\rVert_2\le\lVert g\rVert_2/\mu$ であり、 $-g^{\mathsf T}p\ge\frac{\mu}{\lVert J(x)\rVert_2^2+\mu}\lVert g\rVert_2\lVert p\rVert_2$ が成り立つ。
$\mu\to\infty$ のとき $\mu p_\mu(x)\to-g$ である。 $x\in D_{\mathrm{GN}}$ ならば、 $\mu\to+0$ のとき $p_\mu(x)\to s(x)$ である。

証明. $J:=J(x)$ 、 $M:=J^{\mathsf T}J+\mu I_n$ と置く。補題 1.2 (1)により $g=J^{\mathsf T}r(x)$ であり、 $Mp=-g$ 、 $p^{\mathsf T}Mp=\lVert Jp\rVert_2^2+\mu\lVert p\rVert_2^2$ である。

(1)を示す。 $g\ne0$ ならば $Mp=-g$ から $p\ne0$ であり、 $g^{\mathsf T}p=-p^{\mathsf T}Mp$ から第一の等式を得る。 $\frac12\lVert r(x)+Jp\rVert_2^2=\phi(x)+g^{\mathsf T}p+\frac12\lVert Jp\rVert_2^2$ に第一の等式を代入して第二の等式を得る。 $p\ne0$ と $\mu>0$ から二つの量の符号を得る。

(2)を示す。 Cauchy–Schwarz の不等式により $\mu\lVert p\rVert_2^2\le p^{\mathsf T}Mp=-g^{\mathsf T}p\le\lVert g\rVert_2\lVert p\rVert_2$ であり、第一の不等式を得る。§E20.6 補題 2.1 (1)により $\lVert J^{\mathsf T}\rVert_2=\lVert J\rVert_2$ であるから、 $\lVert g\rVert_2=\lVert Mp\rVert_2\le(\lVert J\rVert_2^2+\mu)\lVert p\rVert_2$ であり、

-g^{\mathsf T}p\ge\mu\lVert p\rVert_2^2\ge\frac{\mu}{\lVert J\rVert_2^2+\mu}\lVert g\rVert_2\lVert p\rVert_2

である。

(3)を示す。 $\mu p+g=-J^{\mathsf T}Jp$ であるから、(2)により $\lVert\mu p+g\rVert_2\le\lVert J\rVert_2^2\lVert g\rVert_2/\mu$ であり、右辺は $\mu\to\infty$ のとき $0$ に収束する。 $x\in D_{\mathrm{GN}}$ ならば $J\ne0$ であり、§E20.25 命題 2.4 (3)を $A=J$ 、 $c=-r(x)$ に適用すると、 $\mu\to+0$ のとき $p_\mu(x)=R_{\sqrt\mu}(-r(x))\to J^\dagger(-r(x))=s(x)$ である。▨

命題 4.3. $m$ 、 $n$ 、 $U$ 、 $r$ 、 $J$ 、 $\phi$ を定義 1.1のとおりとし、 $T\in\R^{n\times n}$ を正則行列とする。 $\tilde U:=\{z\in\R^n\mid Tz\in U\}$ と置き、 $\tilde r\colon\tilde U\to\R^m$ を $\tilde r(z):=r(Tz)$ で定める。 $\tilde U$ は開集合であり、 $\tilde r$ は $C^1$ 級である。残差写像 $\tilde r$ について定義 1.3と定義 4.1で定まる $D_{\mathrm{GN}}$ 、 $G$ 、 $p_\mu$ を $\tilde D_{\mathrm{GN}}$ 、 $\tilde G$ 、 $\tilde p_\mu$ と書く。 $z\in\tilde U$ とし、 $x:=Tz$ と置く。

$z\in\tilde D_{\mathrm{GN}}$ であることと $x\in D_{\mathrm{GN}}$ であることは同値であり、そのとき $T\tilde G(z)=G(x)$ である。
実数 $\mu>0$ について、 $T\tilde p_\mu(z)$ は $q\mapsto\lVert r(x)+J(x)q\rVert_2^2+\mu\lVert T^{-1}q\rVert_2^2$ の最小値を与えるただ一つの $q\in\R^n$ であり、 $\bigl(J(x)^{\mathsf T}J(x)+\mu\,T^{-\mathsf T}T^{-1}\bigr)T\tilde p_\mu(z)=-J(x)^{\mathsf T}r(x)$ を満たす。
$T=dI_n$ （ $d\in\R$ 、 $d\ne0$ ）ならば、任意の実数 $\mu>0$ について $T\tilde p_\mu(z)=p_{\mu/d^2}(x)$ である。さらに $\nabla\phi(x)\ne0$ かつ $d^2\ne1$ ならば $T\tilde p_\mu(z)\ne p_\mu(x)$ である。

証明. 演習とする（問題 8.2）。▨

定義 4.4. $m$ 、 $n$ 、 $U$ 、 $r$ 、 $J$ 、 $\phi$ を定義 1.1のとおりとする。初期値 $x_0\in U$ 、初期減衰係数 $\mu_0>0$ 、受理の定数 $\eta\in(0,1)$ 、更新の倍率 $\nu>1$ 、許容量 $\varepsilon_g\ge0$ 、 $\varepsilon_x\ge0$ 、試行の上限 $K\in\NN$ を固定する。 $k=0,1,\dots,K-1$ の順に、 $x_k\in U$ と $\mu_k>0$ から次の試行 $k$ を行う手続きを Levenberg–Marquardt 法 (Levenberg-Marquardt method) という。

$\lVert\nabla\phi(x_k)\rVert_2\le\varepsilon_g$ ならば、 $x_k$ を出力して停止する。この停止理由を「勾配」という。
$p_k:=p_{\mu_k}(x_k)$ 、 $\operatorname{pred}_k:=\phi(x_k)-\frac12\lVert r(x_k)+J(x_k)p_k\rVert_2^2$ と置く。 $x_k+p_k\in U$ であり $\phi(x_k)-\phi(x_k+p_k)\ge\eta\operatorname{pred}_k$ であるならば、試行 $k$ を受理し、 $x_{k+1}:=x_k+p_k$ 、 $\mu_{k+1}:=\mu_k/\nu$ と置く。このとき $\lVert p_k\rVert_2\le\varepsilon_x(1+\lVert x_{k+1}\rVert_2)$ ならば $x_{k+1}$ を出力して停止し、この停止理由を「修正量」という。
試行 $k$ を受理しないならば、試行 $k$ を棄却し、 $x_{k+1}:=x_k$ 、 $\mu_{k+1}:=\nu\mu_k$ と置く。

試行 $K-1$ の後に停止していなければ、 $x_K$ を出力して停止し、この停止理由を「反復上限」という。

命題 4.5. 記号を定義 4.4のとおりとする。

試行 $k$ で定義 4.4 (2)に達したならば $\operatorname{pred}_k>0$ であり、試行 $k$ を受理したならば $\phi(x_{k+1})<\phi(x_k)$ である。したがって、定まった $x_0,x_1,\dots$ について $\phi(x_0)\ge\phi(x_1)\ge\cdots$ である。
$x\in U$ が $\nabla\phi(x)\ne0$ を満たすとする。実数 $\bar\mu>0$ が存在して、任意の実数 $\mu\ge\bar\mu$ について $x+p_\mu(x)\in U$ であり、 $\phi(x)-\phi\bigl(x+p_\mu(x)\bigr)\ge\eta\Bigl(\phi(x)-\frac12\bigl\lVert r(x)+J(x)p_\mu(x)\bigr\rVert_2^2\Bigr)$ が成り立つ。特に、 $j\in\NN$ とし、 $x_k=x$ であって試行 $k,\dots,k+j-1$ がすべて棄却されたならば、 $\nu^{j-1}\mu_k<\bar\mu$ である。

証明.(1)を示す。試行 $k$ で定義 4.4 (2)に達したならば、定義 4.4 (1)で停止していないので $\lVert\nabla\phi(x_k)\rVert_2>\varepsilon_g\ge0$ であり、命題 4.2 (1)により $\operatorname{pred}_k>0$ である。試行 $k$ を受理したならば $\phi(x_k)-\phi(x_{k+1})\ge\eta\operatorname{pred}_k>0$ であり、棄却したならば $x_{k+1}=x_k$ である。

(2)を示す。 $g:=\nabla\phi(x)$ 、 $\varepsilon:=(1-\eta)\lVert g\rVert_2/2$ と置くと $\varepsilon>0$ である。補題 1.2 (1)により $\phi$ は $x$ で全微分可能であり、§E20.2 補題 1.2 (3)により、実数 $\rho>0$ が存在して、 $\lVert h\rVert_2<\rho$ を満たす任意の $h\in\R^n$ について $x+h\in U$ かつ $\lvert\phi(x+h)-\phi(x)-g^{\mathsf T}h\rvert\le\varepsilon\lVert h\rVert_2$ である。 $\bar\mu:=\max\{\lVert J(x)\rVert_2^2,\ 2\lVert g\rVert_2/\rho\}$ と置くと $\bar\mu>0$ である。 $\mu\ge\bar\mu$ とし、 $p:=p_\mu(x)$ 、 $\pi:=\phi(x)-\frac12\lVert r(x)+J(x)p\rVert_2^2$ と置く。命題 4.2 (2)により $\lVert p\rVert_2\le\lVert g\rVert_2/\mu\le\rho/2<\rho$ であるから $x+p\in U$ であり、 $\mu\ge\lVert J(x)\rVert_2^2$ から $\mu/(\lVert J(x)\rVert_2^2+\mu)\ge1/2$ であるので $-g^{\mathsf T}p\ge\frac12\lVert g\rVert_2\lVert p\rVert_2$ である。命題 4.2 (1)の二つの等式から $\pi=-g^{\mathsf T}p-\frac12\lVert J(x)p\rVert_2^2\le-g^{\mathsf T}p$ である。したがって

\phi(x)-\phi(x+p)\ge-g^{\mathsf T}p-\varepsilon\lVert p\rVert_2=(1-\eta)(-g^{\mathsf T}p)+\eta(-g^{\mathsf T}p)-\varepsilon\lVert p\rVert_2\ge\Bigl(\frac{1-\eta}2\lVert g\rVert_2-\varepsilon\Bigr)\lVert p\rVert_2+\eta\pi=\eta\pi

である。

$j\in\NN$ とし、 $x_k=x$ であって試行 $k,\dots,k+j-1$ がすべて棄却されたとする。定義 4.4 (3)により $x_{k+j-1}=x$ 、 $\mu_{k+j-1}=\nu^{j-1}\mu_k$ である。 $\mu_{k+j-1}\ge\bar\mu$ ならば、上で示したことにより試行 $k+j-1$ は受理され、棄却されたことと両立しない。したがって $\nu^{j-1}\mu_k<\bar\mu$ である。▨

5 停止と局所感度

系 5.1. $m$ 、 $n$ 、 $U$ 、 $r$ 、 $J$ 、 $\phi$ を定義 1.1のとおりとし、 $r$ は $C^2$ 級であるとする。 $x_*\in U$ は $\nabla\phi(x_*)=0$ を満たし、 $H_*:=\nabla^2\phi(x_*)$ は正則であるとする。実数 $R>0$ 、 $\gamma\ge0$ について $\overline B(x_*,R)\subseteq U$ であり、任意の $x,y\in\overline B(x_*,R)$ に対して $\lVert\nabla^2\phi(x)-\nabla^2\phi(y)\rVert_2\le\gamma\lVert x-y\rVert_2$ が成り立つとする。 $\beta:=\lVert H_*^{-1}\rVert_2$ と置き、 $\gamma>0$ ならば $\rho:=\min\{R,1/(2\beta\gamma)\}$ 、 $\gamma=0$ ならば $\rho:=R$ と置く。

任意の $x\in\overline B(x_*,\rho)$ について $\lVert x-x_*\rVert_2\le\frac43\beta\lVert\nabla\phi(x)\rVert_2$ であり、 $\overline B(x_*,\rho)$ に属する $\phi$ の停留点は $x_*$ だけである。特に、定義 4.4の手続きが停止理由「勾配」で出力した点 $x$ が $\overline B(x_*,\rho)$ に属するならば、 $\lVert x-x_*\rVert_2\le\frac43\beta\varepsilon_g$ である。
$r(x_*)=0$ ならば、 $H_*=J(x_*)^{\mathsf T}J(x_*)$ であり、 $H_*$ が正則であることは $J(x_*)$ の列が一次独立であることと同値であって、 $\beta=\sigma_n(J(x_*))^{-2}$ である。

証明.(1)を示す。 $F:=\nabla\phi\colon U\to\R^n$ は補題 1.2 (2)により $C^1$ 級であり、 $DF(x)=\nabla^2\phi(x)$ である。仮定により $x_*$ は $F$ の正則零点で、半径 $R$ 、定数 $\gamma$ の Lipschitz 条件を満たし、§E20.23 定義 2.1の $\beta$ 、 $\rho$ はここでの $\beta$ 、 $\rho$ に一致する。§E20.23 命題 3.1 (1)と§E20.23 命題 3.1 (3)により前半の二つの主張を得る。停止理由「勾配」で出力された点 $x$ は $\lVert\nabla\phi(x)\rVert_2\le\varepsilon_g$ を満たす。

(2)を示す。 $r(x_*)=0$ ならば補題 1.2 (2)の和は $0$ であり、 $H_*=J(x_*)^{\mathsf T}J(x_*)$ である。§E20.6 命題 1.2 (2)により、 $J(x_*)$ の列が一次独立ならば $H_*$ は正則である。一次従属ならば、 $J(x_*)h=0$ を満たす $h\ne0$ について $H_*h=0$ であり、 $H_*$ は正則でない。列が一次独立なとき、§E20.6 補題 2.1 (2)により $\beta=\sigma_n(J(x_*))^{-2}$ である。▨

命題 5.2. $m,n\in\NN$ 、 $U\subseteq\R^n$ を開集合、 $f\colon U\to\R^m$ を $C^2$ 級の写像とし、 $J(x):=Df(x)$ と置く。 $y\in\R^m$ について、残差写像 $x\mapsto f(x)-y$ の目的関数を $\phi_y(x):=\frac12\lVert f(x)-y\rVert_2^2$ と書く。 $y_0\in\R^m$ と $x_*\in U$ が $J(x_*)^{\mathsf T}\bigl(f(x_*)-y_0\bigr)=0$ を満たし、

H_*:=J(x_*)^{\mathsf T}J(x_*)+\sum_{i=1}^m\bigl(f_i(x_*)-y_{0,i}\bigr)\nabla^2f_i(x_*)

が正則であるとする。

$y_0$ の開近傍 $V\subseteq\R^m$ 、 $x_*$ の開近傍 $B\subseteq U$ 、 $C^1$ 級写像 $\xi\colon V\to B$ が存在して、 $\xi(y_0)=x_*$ であり、任意の $y\in V$ について $\nabla\phi_y(\xi(y))=0$ が成り立ち、 $x\in B$ 、 $y\in V$ 、 $\nabla\phi_y(x)=0$ ならば $x=\xi(y)$ である。
$D\xi(y_0)=H_*^{-1}J(x_*)^{\mathsf T}$ である。
$f(x_*)=y_0$ ならば、 $J(x_*)$ の列は一次独立であり、 $D\xi(y_0)=J(x_*)^\dagger$ 、 $\lVert D\xi(y_0)\rVert_2=1/\sigma_n(J(x_*))$ である。さらに、正規直交基底 $u_1,\dots,u_m$ 、 $v_1,\dots,v_n$ と実数 $\sigma_1\ge\dots\ge\sigma_n\ge0$ が $J(x_*)=\sum_{i=1}^n\sigma_iu_iv_i^{\mathsf T}$ を満たすならば、 $1\le i\le n$ について $D\xi(y_0)u_i=\sigma_i^{-1}v_i$ である。

証明. $Q:=U\times\R^m$ は開集合である。 $\Gamma\colon Q\to\R^n$ を $\Gamma(x,y):=J(x)^{\mathsf T}\bigl(f(x)-y\bigr)$ で定めると、補題 1.2 (1)を残差写像 $f-y$ に適用して $\Gamma(x,y)=\nabla\phi_y(x)$ である。 $f$ は $C^2$ 級であるから $\Gamma$ は $C^1$ 級であり、補題 1.2 (2)により $x$ に関する微分は $\nabla^2\phi_y(x)=J(x)^{\mathsf T}J(x)+\sum_i(f_i(x)-y_i)\nabla^2f_i(x)$ 、 $\Gamma$ は $y$ についてアフィンであって $y$ に関する微分は $-J(x)^{\mathsf T}$ である。 $\Gamma(x_*,y_0)=0$ であり、 $x$ に関する微分の $(x_*,y_0)$ での値 $H_*$ は正則である。

(1)を示す。§E20.22 命題 5.1を $G=\Gamma$ 、 $(u_0,p_0)=(x_*,y_0)$ に適用し、§E20.22 命題 5.1 (1)の $V$ 、 $B$ 、 $u$ をとる。 $p\in V$ について $(u(p),p)\in Q$ であるから $u(p)\in U$ であり、 $B$ を $B\cap U$ に、 $u$ を $\xi$ に置き換えて主張を得る。

(2)を示す。§E20.22 命題 5.1 (2)により、任意の $\dot y\in\R^m$ について $H_*D\xi(y_0)\dot y=J(x_*)^{\mathsf T}\dot y$ であり、 $H_*$ は正則である。

(3)を示す。 $f(x_*)=y_0$ ならば $H_*=J(x_*)^{\mathsf T}J(x_*)$ であり、 $J(x_*)h=0$ ならば $H_*h=0$ であるから、 $H_*$ の正則性により $h=0$ である。したがって $J(x_*)$ の列は一次独立であり、§E20.9 命題 6.1 (3)により $D\xi(y_0)=H_*^{-1}J(x_*)^{\mathsf T}=J(x_*)^\dagger$ 、 $\lVert D\xi(y_0)\rVert_2=1/\sigma_n(J(x_*))$ である。 $J(x_*)\ne0$ であるから $\sigma_1>0$ であり、§E20.25 補題 1.1 (1)により $\sigma_i>0$ を満たす $i$ の個数は $\operatorname{rank}J(x_*)=n$ であって、 $J(x_*)^\dagger=\sum_{i=1}^n\sigma_i^{-1}v_iu_i^{\mathsf T}$ である。 $u_1,\dots,u_m$ の正規直交性から最後の等式を得る。▨

例 5.3. $h(x,t):=a_1e^{-b_1t}+a_2e^{-b_2t}$ （ $x=(a_1,b_1,a_2,b_2)\in\R^4$ 、 $t\in\R$ ）とし、 $P(a_1,b_1,a_2,b_2):=(a_2,b_2,a_1,b_1)$ と置く。任意の $x$ と $t$ について $h(Px,t)=h(x,t)$ であるから、任意のデータと重みについて、当てはめの目的関数は $\phi(Px)=\phi(x)$ を満たす。 $x_*:=(1,\log2,1,\log3)$ とし、 $t_i:=i-1$ 、 $y_i:=h(x_*,t_i)=2^{-t_i}+3^{-t_i}$ 、 $w_i:=1$ （ $1\le i\le4$ ）とする。 $\phi(x_*)=\phi(Px_*)=0$ であり、 $x_*\ne Px_*$ である。 $J(x_*)$ の第 $i$ 行は $(2^{-t_i},\ -t_i2^{-t_i},\ 3^{-t_i},\ -t_i3^{-t_i})$ であり、有理数の計算により $\det J(x_*)=1/7776$ である。 $J(Px_*)$ は $J(x_*)$ の列を並べ替えた行列であるから正則である。したがって命題 5.2は $(x_*,y)$ と $(Px_*,y)$ のいずれにも適用することができ、データから $x_*$ の近くの停留点への写像と $Px_*$ の近くの停留点への写像はいずれも $C^1$ 級であるが、 $\phi$ の最小値 $0$ を与える点は $x_*$ と $Px_*$ の二つを含む。

6 減衰曲線の当てはめ

例 6.1. $h(x,t):=ae^{-bt}+c$ （ $x=(a,b,c)\in U:=\R^3$ ）とし、重み $w_i:=1$ の当てはめを考える。 $J(x)$ の第 $i$ 行は $(e^{-bt_i},\ -at_ie^{-bt_i},\ 1)$ である。観測は次の二つとする。

長い区間: $t_i:=(i-1)/2$ （ $1\le i\le9$ ）、 $y=(2.510,1.831,1.409,1.092,0.914,0.761,0.691,0.612,0.592)$ 。
短い区間: $t_i:=(i-1)/20$ （ $1\le i\le9$ ）、 $y=(2.510,2.412,2.356,2.264,2.214,2.127,2.083,2.002,1.962)$ 。

いずれも $y_i$ は $2e^{-0.8t_i}+0.5+0.01\cdot(-1)^{i-1}$ を小数第 $3$ 位に丸めた値であり、 $x^\circ:=(2,0.8,0.5)$ について $\lVert y-(h(x^\circ,t_i))_i\rVert_2$ は長い区間で $2.999\times10^{-2}$ 、短い区間で $2.947\times10^{-2}$ である。

定義 4.4の手続きを $\mu_0=10^{-3}$ 、 $\nu=10$ 、 $\eta=1/4$ 、 $\varepsilon_g=\varepsilon_x=10^{-10}$ 、 $K=100$ で実行した。演算は binary64 で行い、 $e^{-bt}$ は NumPy 2.0.2 の numpy.exp で計算し、 $p_k$ は拡大行列の QR 分解を numpy.linalg.qr（LAPACK の Householder QR 法）で求めて上三角系を後退代入で解いた。以下の値はこの実行の観察であり、丸めて示す。試行数は停止までに計算した $p_k$ の個数である。

長い区間で $x_0=(1.5,1,0.3)$ とした実行 A の各試行は次のとおりである。

$k$	$x_k$	$\phi(x_k)$	$\lVert\nabla\phi(x_k)\rVert_2$	$\mu_k$	判定
0	$(1.5,\ 1,\ 0.3)$	$9.669\times10^{-1}$	$4.396$	$10^{-3}$	受理
1	$(1.98034,\ 0.73755,\ 0.52468)$	$1.671\times10^{-2}$	$6.462\times10^{-1}$	$10^{-4}$	受理
2	$(1.99549,\ 0.80971,\ 0.51126)$	$4.577\times10^{-4}$	$3.124\times10^{-2}$	$10^{-5}$	受理
3	$(2.00051,\ 0.80972,\ 0.50670)$	$4.085\times10^{-4}$	$8.836\times10^{-8}$	$10^{-6}$	受理
4	$(2.00051,\ 0.80972,\ 0.50670)$	$4.085\times10^{-4}$	$6.994\times10^{-11}$		停止（勾配）

四つの実行の結果は次のとおりである。 $\sigma_i$ は出力点での $J$ の特異値である。

実行	区間	$x_0$	停止理由	試行数	棄却数	出力	$\lVert r\rVert_2$	$(\sigma_1,\sigma_2,\sigma_3)$
A	長い	$(1.5,1,0.3)$	勾配	4	0	$(2.00051,\ 0.80972,\ 0.50670)$	$2.858\times10^{-2}$	$(3.615,\ 1.011,\ 0.5961)$
B	長い	$(1,5,0)$	勾配	10	3	$(2.00051,\ 0.80972,\ 0.50670)$	$2.858\times10^{-2}$	$(3.615,\ 1.011,\ 0.5961)$
C	長い	$(1,-1,0)$	反復上限	100	50	$(-9.2047,\ -0.042584,\ 11.195)$	$0.7764$	$(75.44,\ 2.380,\ 2.813\times10^{-3})$
D	短い	$(1.5,1,0.3)$	勾配	5	0	$(1.58549,\ 1.06121,\ 0.92008)$	$2.823\times10^{-2}$	$(3.939,\ 0.5051,\ 1.510\times10^{-2})$

実行 B では、試行 $0$ の $p_0$ のノルムは $13.58$ であり、 $\phi(x_0+p_0)$ は $\phi(x_0)=5.557$ を上回って棄却された。試行 $1$ を $\mu_1=10^{-2}$ で受理して $x_2=(1.6142,\ 2.0923,\ 0.9447)$ 、 $\phi(x_2)=0.2613$ となり、試行 $2$ 、 $3$ を棄却し、試行 $4$ を $\mu_4=10^{-1}$ で受理した。同じ $x_0$ で $J(x_0)$ の特異値は $(3.024,\ 0.9279,\ 0.03821)$ であり、Gauss–Newton 法の一段を同じ演算条件で計算すると $G(x_0)=(1.6764,\ -17.800,\ 0.83057)$ 、 $\phi(G(x_0))\approx9.8\times10^{61}$ である。
実行 C では、命題 4.5 (1)により $\phi(x_k)$ は $\phi(x_0)=2269$ から増えずに出力で $0.3014$ となり、この値は A の出力の $4.085\times10^{-4}$ より大きい。出力で $\lVert\nabla\phi\rVert_2=7.780\times10^{-2}>\varepsilon_g$ である。
実行 A と D の出力点で、残差のノルムは $2.858\times10^{-2}$ と $2.823\times10^{-2}$ であり、 $\lVert x-x^\circ\rVert_2$ は $1.182\times10^{-2}$ と $0.6454$ である。出力点 $x$ で命題 5.2 (2)の $H_*^{-1}J^{\mathsf T}$ を $y_0=y$ 、 $x_*=x$ として計算すると、その 2 ノルムは A で $1.680$ 、D で $66.28$ であり、 $1/\sigma_3$ は A で $1.678$ 、D で $66.24$ である。二つの値の差は $H_*$ の残差の項による。
観測時刻は、長い区間で $h=1/2$ 、短い区間で $h=1/20$ として $t_i=(i-1)h$ である。 $q:=e^{-bh}$ と置くと、 $J(x)$ の最初の $3$ 行からなる行列は $\begin{pmatrix}1&0&1\\q&-ahq&1\\q^2&-2ahq^2&1\end{pmatrix}$ であり、その行列式は $-ahq(q-1)^2$ である。実行 A と D の出力はいずれも $a\ne0$ 、 $b\ne0$ を満たし、 $q>0$ 、 $q\ne1$ であるから、この行列式は $0$ でなく、出力点で $J$ の列は一次独立であって $\operatorname{rank}J=3$ である。許容量 $0.02$ に関する§E20.9 定義 5.1の数値的階数は、表の $(\sigma_1,\sigma_2,\sigma_3)$ のうち $0.02$ より大きいものの個数であり、A で $3$ 、D で $2$ である。この判定は、観測時刻、パラメータ $(a,b,c)$ の尺度、許容量 $0.02$ を固定したものである。A と D はいずれも停止理由「勾配」で出力したが、許容量 $0.02$ に関する数値的階数は異なり、出力点での $H_*^{-1}J^{\mathsf T}$ の 2 ノルムは A で $1.680$ 、D で $66.28$ である。
出力点 $x$ で命題 3.1の $-(J^{\mathsf T}J)^{-1}S_*$ を計算すると、その 2 ノルムは A で $1.787\times10^{-3}$ 、D で $1.184\times10^{-3}$ である。

7 統計的分散との比較

命題 7.1. $m\ge n$ とし、 $X\in\R^{m\times n}$ の列は一次独立であるとする。正規直交基底 $u_1,\dots,u_m$ 、 $v_1,\dots,v_n$ が $X=\sum_{i=1}^n\sigma_i(X)u_iv_i^{\mathsf T}$ を満たすとする。

任意の $e\in\R^m$ と $1\le i\le n$ について $v_i^{\mathsf T}X^\dagger e=\sigma_i(X)^{-1}u_i^{\mathsf T}e$ である。実数 $\delta\ge0$ について、 $\lVert e\rVert_2\le\delta$ を満たす $e\in\R^m$ の全体での $\lVert X^\dagger e\rVert_2$ の最大値は $\delta/\sigma_n(X)$ であり、 $e=\delta u_n$ で達する。
$y=X\beta+\varepsilon$ を§E14.19 定義 1.1の線形回帰モデルとし、Gauss–Markov 仮定 $E[\varepsilon]=0$ 、 $\operatorname{Cov}(\varepsilon)=\sigma^2I_m$ を置く。 $\hat\beta:=X^\dagger y$ は $E[\hat\beta]=\beta$ 、 $\operatorname{Cov}(\hat\beta)=\sigma^2\sum_{i=1}^n\sigma_i(X)^{-2}v_iv_i^{\mathsf T}$ を満たし、 $1\le i\le n$ について $\operatorname{Var}(v_i^{\mathsf T}\hat\beta)=\sigma^2/\sigma_i(X)^2$ 、また $E\lVert\hat\beta-\beta\rVert_2^2=\sigma^2\sum_{i=1}^n\sigma_i(X)^{-2}$ である。

証明.(1)を示す。 $\sigma_i(X)>0$ （ $1\le i\le n$ ）であるから、§E20.25 補題 1.1 (1)により $X^\dagger=\sum_{i=1}^n\sigma_i(X)^{-1}v_iu_i^{\mathsf T}$ であり、 $v_1,\dots,v_n$ の正規直交性から第一の等式を得る。§E20.25 補題 1.1 (2)を $g_i=\sigma_i(X)^{-1}$ に適用すると、 $\lvert g_i\rvert$ の最大値は $i=n$ でとられるから、第二の主張を得る。

(2)を示す。§E20.9 命題 6.1 (3)により $\hat\beta=(X^{\mathsf T}X)^{-1}X^{\mathsf T}y$ は§E14.19 系 2.2の最小二乗推定量であり、§E14.19 系 4.1により $E[\hat\beta]=\beta$ 、 $\operatorname{Cov}(\hat\beta)=\sigma^2(X^{\mathsf T}X)^{-1}$ である。 $X^{\mathsf T}X=\sum_{i=1}^n\sigma_i(X)^2v_iv_i^{\mathsf T}$ であり、 $v_1,\dots,v_n$ は $\R^n$ の正規直交基底であるから、 $(X^{\mathsf T}X)^{-1}=\sum_{i=1}^n\sigma_i(X)^{-2}v_iv_i^{\mathsf T}$ である。 $\operatorname{Var}(v_i^{\mathsf T}\hat\beta)=v_i^{\mathsf T}\operatorname{Cov}(\hat\beta)v_i=\sigma^2\sigma_i(X)^{-2}$ であり、 $E[\hat\beta]=\beta$ から $E\lVert\hat\beta-\beta\rVert_2^2=\operatorname{tr}\operatorname{Cov}(\hat\beta)=\sigma^2\sum_{i=1}^n\sigma_i(X)^{-2}$ である。▨

例 7.2. 直線 $\beta_1+\beta_2t$ の当てはめで、 $X$ の第 $i$ 行を $(1,t_i)$ （ $1\le i\le9$ ）とし、例 6.1と同じ長い区間 $t_i=(i-1)/2$ と短い区間 $t_i=(i-1)/20$ を比べる。

長い区間では $X^{\mathsf T}X=\begin{pmatrix}9&18\\18&51\end{pmatrix}$ 、 $(X^{\mathsf T}X)^{-1}=\frac1{135}\begin{pmatrix}51&-18\\-18&9\end{pmatrix}$ であり、 $X^{\mathsf T}X$ の固有値 $30\pm3\sqrt{85}$ から $\sigma_1(X)\approx7.593$ 、 $\sigma_2(X)\approx1.530$ である。
短い区間では $X^{\mathsf T}X=\begin{pmatrix}9&9/5\\9/5&51/100\end{pmatrix}$ 、 $(X^{\mathsf T}X)^{-1}=\begin{pmatrix}17/45&-4/3\\-4/3&20/3\end{pmatrix}$ であり、 $\sigma_1(X)\approx3.060$ 、 $\sigma_2(X)\approx0.3797$ である。

命題 7.1 (2)により、傾きの分散 $\operatorname{Var}(\hat\beta_2)$ は長い区間で $\sigma^2/15$ 、短い区間で $20\sigma^2/3$ であり、 $E\lVert\hat\beta-\beta\rVert_2^2$ は $4\sigma^2/9$ と $317\sigma^2/45\approx7.044\sigma^2$ である。命題 7.1 (1)により、 $\lVert e\rVert_2\le\delta$ の範囲での $\lVert X^\dagger e\rVert_2$ の最大値は $0.6535\,\delta$ と $2.634\,\delta$ である。 $E\lVert\varepsilon\rVert_2^2=9\sigma^2$ であり、 $\delta^2=9\sigma^2$ とすると最大値の二乗は $3.844\sigma^2$ と $62.44\sigma^2$ であって、 $E\lVert\hat\beta-\beta\rVert_2^2$ のそれぞれ約 $8.65$ 倍と $8.86$ 倍である。

注意 7.3. $f$ 、 $J$ 、 $\phi_y$ を命題 5.2のとおりとし、 $x^\circ\in U$ で $J(x^\circ)$ の列が一次独立であるとして、 $y_0:=f(x^\circ)$ 、 $x_*:=x^\circ$ とする。 $H_*=J(x^\circ)^{\mathsf T}J(x^\circ)$ は正則であり、同命題の $V$ 、 $\xi$ が定まる。観測 $y=y_0+\varepsilon$ の $\varepsilon$ に Gauss–Markov 仮定を置く。命題 5.2 (3)の一次の項 $D\xi(y_0)\varepsilon=J(x^\circ)^\dagger\varepsilon$ は、線形回帰モデル $y'=J(x^\circ)x^\circ+\varepsilon$ の $\hat\beta=J(x^\circ)^\dagger y'$ について $\hat\beta-x^\circ$ に等しいので、命題 7.1 (2)を $X=J(x^\circ)$ に適用して、平均 $0$ 、共分散 $\sigma^2(J(x^\circ)^{\mathsf T}J(x^\circ))^{-1}$ をもつ。 $y_0+\varepsilon$ は $V$ に属するとは限らず、 $\xi(y_0+\varepsilon)$ が定まるとは限らない。

$\varepsilon$ の成分が独立に正規分布 $N(0,\sigma^2)$ に従うならば、 $x$ における尤度は $(2\pi\sigma^2)^{-m/2}\exp(-\phi_y(x)/\sigma^2)$ であり、尤度の最大値を与える $x$ と $\phi_y$ の最小値を与える $x$ は一致する。重み $1$ の当てはめで $y_i$ の分布は $N(h(x,t_i),\sigma^2)$ であり、 $h(x,t_1),\dots,h(x,t_m)$ が一定でなければ $y_1,\dots,y_m$ は同分布でないので、§E14.13 定理 3.1の独立同分布の仮定を満たさない。例 5.3の $P$ について、尤度は $x$ と $Px$ で等しい。

8 演習

問題 8.1.補題 1.2の証明を完成させよ。

解答.

$\phi=\frac12\sum_{i=1}^mr_i^2$ であり、各 $r_i$ は $C^1$ 級であるから、 $1\le j\le n$ について $\partial_j\phi=\sum_{i=1}^mr_i\,\partial_jr_i$ は $U$ 上で連続である。したがって $\phi$ は $C^1$ 級であり、 $\nabla\phi(x)$ の第 $j$ 成分 $\sum_i\partial_jr_i(x)\,r_i(x)$ は $J(x)^{\mathsf T}r(x)$ の第 $j$ 成分である。これで補題 1.2 (1)は示された。

$r$ が $C^2$ 級ならば、 $1\le j,k\le n$ について

\partial_k\partial_j\phi=\sum_{i=1}^m\bigl(\partial_kr_i\,\partial_jr_i+r_i\,\partial_k\partial_jr_i\bigr)

は $U$ 上で連続であり、 $\phi$ は $C^2$ 級である。右辺の第一項の和は $J(x)^{\mathsf T}J(x)$ の $(k,j)$ 成分であり、第二項の和は $\sum_ir_i(x)\nabla^2r_i(x)$ の $(k,j)$ 成分である。これで補題 1.2 (2)は示された。▨

問題 8.2.命題 4.3の証明を完成させよ。

解答.

$z\mapsto Tz$ は連続であるから $\tilde U$ は開集合であり、全微分の連鎖律により $\tilde r$ は $C^1$ 級であって $\tilde J(z):=D\tilde r(z)=J(x)T$ である。

命題 4.3 (1)を示す。 $T$ は正則であるから、 $h\in\R^n$ について $J(x)Th=0$ であることと $Th\in\ker J(x)$ であることは同値であり、 $\ker\tilde J(z)=T^{-1}\ker J(x)$ である。したがって $\tilde J(z)$ の列が一次独立であることと $J(x)$ の列が一次独立であることは同値である。そのとき

\tilde s(z)=-(T^{\mathsf T}J(x)^{\mathsf T}J(x)T)^{-1}T^{\mathsf T}J(x)^{\mathsf T}r(x)=-T^{-1}(J(x)^{\mathsf T}J(x))^{-1}J(x)^{\mathsf T}r(x)=T^{-1}s(x)

であり、 $T\tilde G(z)=Tz+s(x)=G(x)$ である。

命題 4.3 (2)を示す。 $p\in\R^n$ について $\lVert\tilde r(z)+\tilde J(z)p\rVert_2^2+\mu\lVert p\rVert_2^2=\lVert r(x)+J(x)Tp\rVert_2^2+\mu\lVert p\rVert_2^2$ である。定義 4.1により $\tilde p_\mu(z)$ は左辺の最小値を与えるただ一つの $p$ であり、 $p\mapsto Tp$ は $\R^n$ の全単射であるから、 $T\tilde p_\mu(z)$ は $\lVert r(x)+J(x)q\rVert_2^2+\mu\lVert T^{-1}q\rVert_2^2$ の最小値を与えるただ一つの $q$ である。定義 4.1の一次方程式 $(T^{\mathsf T}J(x)^{\mathsf T}J(x)T+\mu I_n)\tilde p_\mu(z)=-T^{\mathsf T}J(x)^{\mathsf T}r(x)$ の両辺に左から $T^{-\mathsf T}$ を掛けると、 $T^{-\mathsf T}\tilde p_\mu(z)=T^{-\mathsf T}T^{-1}T\tilde p_\mu(z)$ から主張の等式を得る。

命題 4.3 (3)を示す。 $T=dI_n$ ならば $T^{-\mathsf T}T^{-1}=d^{-2}I_n$ であり、命題 4.3 (2)の等式は $(J(x)^{\mathsf T}J(x)+(\mu/d^2)I_n)T\tilde p_\mu(z)=-J(x)^{\mathsf T}r(x)$ となる。定義 4.1によりこの一次方程式のただ一つの解は $p_{\mu/d^2}(x)$ である。 $\nabla\phi(x)\ne0$ 、 $d^2\ne1$ とし、 $p:=p_{\mu/d^2}(x)=p_\mu(x)$ と仮定する。二つの一次方程式の差をとると $(\mu/d^2-\mu)p=0$ であり、 $\mu/d^2\ne\mu$ から $p=0$ である。一方命題 4.2 (1)により $p_\mu(x)\ne0$ であり、 $p=0$ と両立しない。したがって $T\tilde p_\mu(z)=p_{\mu/d^2}(x)\ne p_\mu(x)$ である。▨

§E20.26非線形最小二乗法とパラメータ推定