非線形最小二乗法のNewton法およびGauss-Newton法について

最小二乗法という本を読んでいて、2次元での非線形最小二乗法のNewton法およびGauss-Newton法のアルゴリズムが、難しくてわかりません。
例題でのパラメータ推定をやってみたいので、ご専門の方ご教示くださいませ。

No.3 回答者： echoes_x86 回答日時： 2009/01/25 15:45

こんにちは．

経験者で恐縮ですが…
(専門家の方ほどこのようなところで「専門家」を名乗らないのでは？)
さて，2次元での非線型最小二乗法とのことですので，未知数はωとφでしょう．
これらそれぞれについて1から19の推定値が得られる，というの質問者様の誤読である可能性が高いです．
また，Taylor展開するのは「二乗誤差関数」か「誤差関数の勾配」です．
いずれにせよ誤差関数の2次微分まで考えることになります(Newton法の導出法は複数あります)．

ここでは前者の方針を取ります．
いま，誤差ベクトルを e = y - f(x)とするとき，誤差ベクトルの要素の平方和は
e^T*eと書けます(^Tは転置記号)．
後で見易くするために評価関数 g(x) = e^T*e/2を考え，評価関数をxの近傍で2次近似します．
これは一般の関数の最小値を計算するのは困難であるためです．
g(x+Δ) → g(x) + g'*Δ + Δ^T*g''*Δ ...(1)
ここで，g'とg''はそれぞれgの1次微分と2次微分です(双方共に行列です)．
また，xとΔがベクトルであることに注意してください．
次に(1)式のΔを操作して(つまりxの周りを調べて)右辺を最小化することを考えます．
幸い，これはΔの2次関数ですから，
細かいことを考えなければ(後述)微分して零の点が近くにある最小値(=極小値)と期待されます．
g' + g''*Δ = 0 ...(2)
ここで，g'について考えます．これはeの二乗ですから，
e = y - f(x)であることに注意すると，
[e^T*e/2]' = (e')^T*e = [(y-f(x))']^T*e = [-f'(x)]^T*e = -J^T*e ...(3)
ただし，Jはf'(x)，Jacobianです．
g''のために，(3)をもう一度微分します．
g'' = [(e')^T*e]' = (e'')^T*e + J^T*J = H ...(4)
この式はHessianと呼ばれます．
最終的に (3)と(4)を(2)に代入すると線型方程式
H*Δ = -J^T*e
を得ます．
これを解けば極小値の位置が分かるので x ← x + Δ と更新します．
注意しなければならないのはここでわかるのは「2次近似した関数の極小値」である点です．
近似が良くなければ x + Δ は元の関数の極小値にならないので以上の手続きを反復することになります．

ここで問題なのが，「細かいことを考えなければ」の部分です．
2変数(以上の)関数の場合，「微分して零」の点は鞍点となり極値でない場合があります．
図で書くと近似曲面が下の図の左の形をして欲しいのに，
右の形になることがあります(黒い円付近が「微分して零の点」)．
右の図の円は明らかに極小ではありませんから，
折角見つけたΔも意味をなしません．
ここで，(4)式のHが「正定値」であれば求める点は極小値となります．

Gauss-Newton法は上記の問題を解決するものです．
余程運が悪くなければJ^T*Jは正定値ですから，(e''(x))^T*e の項を無視して，
HをJ^T*Jで近似してしまえば良いのです．
近似が良くなり立つのはeかe''(x)の少なくとも一方が零に近い場合です．
要するに，
・誤差が十分に小さいとき
・e''(x)が零に近い(評価関数が線型に近い)
ときにNewton法に近い振る舞いをする(しかも扱いやすい)のがGauss-Newton法です．

この回答へのお礼

echoes_x86さま

回答および画像までいただきありがとうございました。
Newton法もGauss-Newton法も一括最小二乗法なのですね。
未知パラメータxの時間履歴が気になったので，文献を読みまして，
勉強しているところです。
ご参考にさせていただきます。ありがとうございました。

お礼日時：2009/01/26 18:22

No.2 回答者： arrysthmia 回答日時： 2009/01/23 20:17

←No.1 補足

そうなんでしょうか？ 当方、何ぶん専門家でないもので…

しかし、
与えられたデータ (t, y[t]) が 19 組。これを代入して、
方程式 y[t]　=　sin(ω[t] t + φ[t]) + e[t] が 19 連立。

それに対して、最適化するパラメータが
ω_hat[1]，ω_hat[2]，ω_hat[3], …, ω_hat[19],
φ_hat[1]，φ_hat[2]，φ_hat[3], …, φ_hat[19] の 38 個
では、最適化も何も、自由度が高すぎるように思いますが…
どうなんでしょうね。

例えば、ω[t] ≡ 0, φ[t] = arcsin(y[t]) とすれば、
誤差 e[t] ≡ 0 にすることができます。

No.1 回答者： arrysthmia 回答日時： 2009/01/23 00:15

専門家じゃないんですが…

y[t]　=　sin(ωt+φ) + e[t]　の間違いではないでしょうか。
Σ{ e[t] }^2 を最小にするような ω, φ を求めたい ということですよね？
ω, φ が t の関数では、求めようがないように思います。

極小点を求めるために、連立方程式
0　=　(∂/∂ω) Σ{ e[t] }^2　=　2 Σ{ ( y[t] - sin(ωt+φ) ) cos(ωt+φ) t },
0　=　(∂/∂φ) Σ{ e[t] }^2　=　2 Σ{ ( y[t] - sin(ωt+φ) ) cos(ωt+φ) }.
を ω, φ について解けばよいのですが、
非線型方程式なので、厳密解を得ることは期待薄です。
そこで、方程式の近似解法として、一般的な
Newton法, Gauss-Newton法 を使おうという訳です。

方程式
0　=　f(θ,φ),
0　=　g(θ,φ).　を Newton法で解くには、
まづ、何らかの近似解
0　≒　f(a,b),
0　≒　g(a,b).　を得た後、
f( ), g( ) のテーラー展開を一次で打ち切った
0　=　f(a,b) + { (∂f/∂θ)(a,b) } (θ - a) + { (∂f/∂φ)(a,b) } (φ - b),
0　=　g(a,b) + { (∂g/∂θ)(a,b) } (θ - a) + { (∂g/∂φ)(a,b) } (φ - b).
の解 θ, φ のほうが、a, b よりマシな近似だと考えます。
得られた θ, φ を次の a, b として、この操作を繰り返すと、
最初の a, b が良ければ、回数→∞ の極限で
近似解が厳密解に収束することが知られています。

今回の問題では、この漸化式は、
0　=　(θ - a) A + (φ - b) B,
0　=　(θ - a) B + (φ - b) C,
A　=　Σ{ ( 1 - y[t] sin(at+b) ) t^2 },
B　=　Σ{ ( 1 - y[t] sin(at+b) ) t },
C　=　Σ{ 1 - y[t] sin(at+b) }.　になります。
連立一次方程式を解けば、a, b から θ, φ を求めることができます。
θ, φ が適当に a, b に近くなるまで、これを繰り返せば完了です。

この回答への補足

ご解答いただきありがとうございます。

>y[t]　=　sin(ωt+φ) + e[t]　の間違いではないでしょうか。
>Σ{ e[t] }^2 を最小にするような ω, φ を求めたい
　ということですよね？
>ω, φ が t の関数では、求めようがないように思います。

上記の部分ですが，
最小二乗法による実験データ解析―プログラムSALS (UP応用数学選書 7)
というのを県立図書館で借りてきたところ
　　t＝1，2，3…19のときに
ωの推定値ω_hatとφの推定値φ_hatは
　　ω_hat[1]，ω_hat[2]，ω_hat[3]…ω_hat[19]
　　φ_hat[1]，φ_hat[2]，φ_hat[3]…φ_hat[19]
となるような解説が書かれておりました。
（あくまで、私の解釈でございます。）

また，未知の値をベクトル化して
　　θ[t]＝[ω[t] φ[t]]T
ただし， [ω[t] φ[t]]T は [ω[t] φ[t]] の転置です。
そして
　y[t]　=　g(θ[t]) + e[t]
として，g(θ[t])をテイラー展開して
　　y[t]　≒　g(θ_hat[t-1]) + (∂g/∂θ_hat[t-1])*(θ[t]-θ_hat[t-1]) + e[t]

とすればよいと解釈したのですが，
これで果して良いのか？
あっていればこれからどうやって進めていくのか？
理解できずに困っている次第です。

補足日時：2009/01/23 13:26

この回答へのお礼

arrysthmiaさま
コメントいただきありがとうございました。
参考にさせていただきますね

お礼日時：2009/01/23 14:02