レーベンバーグ・マルカート法

レーベンバーグ・マルカート法（レーベンバーグ・マルカートほう、テンプレート:Lang-en-short、LM法、レーベンバーグ・マーカート法）とは、数学および計算科学における非線形最小二乗問題の解法の一つをいう。特に曲線回帰を最小二乗法により行う場合によく用いられる。LM法はガウス・ニュートン法（GN法）と最急降下法を内挿した手法といえる。GN法よりもロバストであり、初期値が解から大きく外れていた場合でも解けることが多いが、ふるまいの良い関数に対してはGN法よりも収束が遅い傾向を示す。LM法は GN法に信頼領域法を適用したものとみることもできる。

1944年、テンプレート:Ill2職員テンプレート:Ill2により初発表され、1963年にデュポン勤務の統計家、テンプレート:仮リンクにより再発見された。また、Girard, Wynne, Morrison によりそれぞれ独立に再発見されている。

LM法は、一般的な曲線回帰問題を解く必要のあるアプリケーションソフトウェアで広く用いられる。GN法を取り入れているため多くの場合で1次解法よりも収束速度が速いが、他の反復法と同様、LM法で保証されているのは局所最小値への収束のみであり、大域最小値が得られる保証はない。

問題設定

LM法の主な適用対象は、最小二乗曲線回帰問題である。 テンプレート:Mvar対の独立変数と従属変数のペア$(x_i,y_i)$が与えられたとき、テンプレート:Mvarをパラメータとする曲線テンプレート:Mathと所与のペアとの偏差の二乗和テンプレート:Mathを最小化することを考える。

$\hat{\mathit{\beta }}\in \underset{\mathit{\beta }}{\mathrm{argmin}}S\left(\mathit{\beta }\right)\equiv \underset{\mathit{\beta }}{\mathrm{argmin}}\sum _{i=1}^m{[y_i-f(x_i,\mathit{\beta })]}^2$

解法

他の数値最適化手法と同様、LM法は反復法を用いる。まず、パラメーターベクトルテンプレート:Mvarの初期推定値を与える必要がある 。極小点が1つしかない場合、事前情報に基づかない均一な初期推定値、たとえば テンプレート:Mathでも大域解に到達することができるが、複数の局所最小値が存在する場合、初期推定値が十分に大域最小点に近いときにしか大域解には収束しない。

各反復ステップにおいてパラメーターベクトルテンプレート:Mvarは新しい推定値テンプレート:Mathへ置き換えられる。テンプレート:Mvarを決めるため、テンプレート:Mathを次のように線形近似する。

$f(x_i,\mathit{\beta }+\mathit{\delta })\approx f(x_i,\mathit{\beta })+{𝑱}_i\mathit{\delta }$

ここで、

${𝑱}_i=\frac{\partial f(x_i,\mathit{\beta })}{\partial \mathit{\beta }}$

は関数テンプレート:Mvarのテンプレート:Mvarについての勾配（ここでは行ベクトルとする）である。

偏差の二乗和テンプレート:Mathは、この勾配がゼロのとき局所最小となる。上式の一次近似を用いると、テンプレート:Mathの偏差二乗和は以下のように近似される。

$S(\mathit{\beta }+\mathit{\delta })\approx \sum _{i=1}^m{[y_i-f(x_i,\mathit{\beta })-{𝑱}_i\mathit{\delta }]}^2$

ベクトル表記すると、以下のように書ける。

$\begin{array}{cc}S(\mathit{\beta }+\mathit{\delta })& \approx {\Vert 𝒚-𝒇\left(\mathit{\beta }\right)-𝑱\mathit{\delta }\Vert }^2\\ & ={[𝒚-𝒇\left(\mathit{\beta }\right)-𝑱\mathit{\delta }]}^{\mathrm{T}}[𝒚-𝒇\left(\mathit{\beta }\right)-𝑱\mathit{\delta }]\\ & ={[𝒚-𝒇\left(\mathit{\beta }\right)]}^{\mathrm{T}}[𝒚-𝒇\left(\mathit{\beta }\right)]-{[𝒚-𝒇\left(\mathit{\beta }\right)]}^{\mathrm{T}}𝑱\mathit{\delta }-{\left(𝑱\mathit{\delta }\right)}^{\mathrm{T}}[𝒚-𝒇\left(\mathit{\beta }\right)]+{\mathit{\delta }}^{\mathrm{T}}{𝑱}^{\mathrm{T}}𝑱\mathit{\delta }\\ & ={[𝒚-𝒇\left(\mathit{\beta }\right)]}^{\mathrm{T}}[𝒚-𝒇\left(\mathit{\beta }\right)]-2{[𝒚-𝒇\left(\mathit{\beta }\right)]}^{\mathrm{T}}𝑱\mathit{\delta }+{\mathit{\delta }}^{\mathrm{T}}{𝑱}^{\mathrm{T}}𝑱\mathit{\delta }\end{array}$

テンプレート:Mathをテンプレート:Mvarに関して微分した結果を0とすると、以下の式を得る。

$\left({𝑱}^{\mathrm{T}}𝑱\right)\mathit{\delta }={𝑱}^{\mathrm{T}}[𝒚-𝒇\left(\mathit{\beta }\right)]$

ここで、テンプレート:Mvarはヤコビ行列であり、そのテンプレート:Mvar行目はテンプレート:Mvar に等しい。また、テンプレート:Mathはそれぞれ、テンプレート:Mvar行目成分をテンプレート:Mathとするベクトルである。ヤコビ行列は一般的には正方行列ではなく、テンプレート:Mvarをデータ点数、テンプレート:Mvarをベクトルテンプレート:Mvarのサイズとしてテンプレート:Math長方形行列である。行列積テンプレート:Mathはテンプレート:Math正方行列となり、上式はテンプレート:Mvar連立線形方程式であるからこれを解いてテンプレート:Mvarを得ることができる。これをそのまま解くのがガウス・ニュートン法である。

LM法では、この方程式を次のように「減衰」させたものに置き換える。

$({𝑱}^{\mathrm{T}}𝑱+\lambda 𝐈)\mathit{\delta }={𝑱}^{\mathrm{T}}[𝒚-𝒇\left(\mathit{\beta }\right)](\lambda \ge 0)$

ここで、テンプレート:Mathは単位行列である。これを解いて得られるテンプレート:Mvarを用いてパラメータベクトルテンプレート:Mvarの推定値を更新する。

非負の減衰係数テンプレート:Mvarは各反復ごとに調整される。テンプレート:Mvarが急速に減少する際には小さい値が用いられ、LM法はGN法に近づく。対して、残差が十分に減少しない場合は大きい値のテンプレート:Mvarが用いられ、テンプレート:Mvarのテンプレート:Mvarについての勾配はテンプレート:Mathであることに注意すると、テンプレート:Mvarが大きいときテンプレート:Mvarは勾配の逆向きに近付きLM法は最急降下法に近づくことがわかる。計算されたテンプレート:Mvarが十分小さくなったとき、もしくは得られたパラメータ推定値テンプレート:Mathに置き換えた際の偏差二乗和の減少が十分に小さくなったときのどちらかの場合に反復は打ち切られ、解テンプレート:Mvarを得る。

減衰係数テンプレート:Mvarがテンプレート:Mathに比べて大きいときは、テンプレート:Mathの逆行列を求める必要はなく、更新ステップテンプレート:Mvarは${\lambda }^{-1}{𝑱}^{\mathrm{T}}[𝒚-𝒇\left(\mathit{\beta }\right)]$で十分に近似される。

LM法は、テンプレート:Mvarが大きい値のときテンプレート:Mathの情報がほとんど使われないという欠点を持つ。Fletcherは1971年、勾配が小さい方向への収束が遅いという問題を避けるため、勾配を曲率に応じてスケールするという考えから、単位行列テンプレート:Mathをテンプレート:Mathの対角要素で置き換え、解をスケール不変にする手法が提案された。

$[{𝑱}^{\mathrm{T}}𝑱+\lambda \mathrm{diag}\left({𝑱}^{\mathrm{T}}𝑱\right)]\mathit{\delta }={𝑱}^{\mathrm{T}}[𝒚-𝒇\left(\mathit{\beta }\right)]$

同様の減衰因子は、線形不良設定問題を解くために用いられるテンプレート:仮リンクや、リッジ回帰と呼ばれる推計法にも現れる。

減衰パラメータの選び方

最良の減衰係数テンプレート:Mvarを選ぶ方法としては、様々な議論があるが、大なり小なりヒューリスティックなものである。それらの選び方がなぜ局所最小点への収束を保証するかを示す理論的な議論はあるが、大域最小点へ収束するような選び方をすると最急降下法の望ましくない特質、とくに収束が遅いという側面が表われてしまう。

どんな選び方をしても、パラメータの大きさはもとの問題がどれほど良くスケールするかに依存する。マーカートは次のような選び方を推奨している。まず初期値テンプレート:Mathを選んで最初のステップを実行し、残差テンプレート:Mathが最初の点よりも減った場合はテンプレート:Mathなる係数を用いて次のステップはテンプレート:Mathとする。残差が増えてしまった場合は、減るようになるまで繰り返しテンプレート:Mvarを掛けテンプレート:Mathを用いて計算をする。

減衰係数テンプレート:Mathを用いた結果が二乗残差を減少させたなら、これをテンプレート:Mvarの新しい値とし（かつテンプレート:Mathを用いた結果を採用し）、プロセスを続行する。もしテンプレート:Mathを用いた残差がテンプレート:Mvarを用いた残差よりも大きくなったならば、テンプレート:Mvarの値を変えず、テンプレート:Mvarを用いた結果を採用する。

delayed gratificationテンプレート:訳語疑問点と呼ばれる減衰係数の効果的な制御方法がある。この方法では、上り坂のステップごとに係数を少しずつ増やし、下り坂のステップごとにパラメーターを大幅に減らす。この戦略は、最適化の開始時に坂を下りすぎ、後に使用できるステップが制限されて、収束が遅くなることを防ぐことを主眼においている。 ほとんどの場合、増加時には2倍、減少時には3分の1を採用すればうまくいくことが示されているが、大規模な問題の場合は、増加時は1.5倍、減少時は5分の1というより極端な値を用いる方がよいことが知られている。

測地線加速度項

レーベンバーグ・マルカート法の更新ステップテンプレート:Mathを、パラメーター空間の測地経路に沿った速度と捉えると、測地経路に沿う加速度に対応する2次の項テンプレート:Mathを次のように加える改善が考えられる。

${𝒗}_k+\frac{1}{2}{𝒂}_k$

ここで、テンプレート:Mathは次の方程式の解である。

${𝑱}_k{𝒂}_k=-f_{vv}({𝒙}_k)$

この測地線加速度項は速度テンプレート:Mathに沿う方向微分 $f_{vv}(𝒙)=\sum _{\mu \nu }{v}_{\mu }{v}_{\nu }{\partial }_{\mu }{\partial }_{\nu }f(𝒙)$のみに依存するため、完全な二次導関数行列を計算する必要はなく、計算コスト上のオーバーヘッドは比較的小さい。2次導関数はかなり複雑な式になる場合があるため、有限差分近似に置き換えると便利な場合がある。

$\begin{array}{cc}{f}_{vv}^i& \approx \frac{f_i(𝒙+h\mathit{\delta })-2f_i(𝒙)+f_i(𝒙-h\mathit{\delta })}{h^2}\\ & =\frac{2}{h}(\frac{f_i(𝒙+h\mathit{\delta })-f_i(𝒙)}{h}-{𝑱}_i\mathit{\delta })\end{array}$

ここで、テンプレート:Mathとテンプレート:Mvarは通常のアルゴリズムでもすでに計算されているため、追加で計算する必要があるのはテンプレート:Mathだけである。有限差分ステップテンプレート:Mvarの選択によってはアルゴリズムが不安定になる場合があるが、通常はおよそ0.1が妥当である。

加速度は速度と反対の方向を指す可能性があり、減衰が小さすぎて収束が失速するのを防ぐために、ステップの採用条件に加速度に関する以下のような追加の基準が追加される。

$\frac{2\Vert {𝒂}_k\Vert }{\Vert {𝒗}_k\Vert }\le \alpha$

ここでテンプレート:Mvarは通常は1未満の固定値とされる。より難しい問題ではより小さい値とする。

測地線加速度項を追加すると、収束速度が大幅に向上する可能性があり、特にアルゴリズムが目的関数ランドスケープ上の狭い峡谷を移動する場合に有用である。このような場合、2次の項の効果によりステップ幅はより小さく、精度が高くなる。

例

テンプレート:Gallery この例では、関数$y=a\cos \left(bX\right)+b\sin \left(aX\right)$ をGNU Octave上のレーベンバーグ・マルカート法実装leasqr関数をもちいてフィッティングする。グラフから、パラメーターフィッティング結果が徐々に改善し、$a=100$ 、 $b=102$ の曲線に近付いていく様子が見てとれる。初期パラメータが元のパラメータに近い場合にのみ、曲線が正確に一致する。この例は、レーベンバーグ・マルカート法が初期条件に非常に敏感であることを示す一例である。その理由の1つは、複数の最小点 (関数) が存在すること、つまりテンプレート:Mathに一致するパラメータの値はテンプレート:Mvarだけでなくテンプレート:Mathと複数あることである。

関連文献

• テンプレート:Cite journal
• テンプレート:Cite journal
• テンプレート:Cite journalテンプレート:Dead link
• テンプレート:Cite book

関連項目

• 信頼領域法
• ネルダー・ミード法

外部リンク

• Detailed description of the algorithm can be found in Numerical Recipes in C, Chapter 15.5: Nonlinear models
• C. T. Kelley, Iterative Methods for Optimization, SIAM Frontiers in Applied Mathematics, no 18, 1999, テンプレート:ISBN2. Online copy
• History of the algorithm in SIAM news
• A tutorial by Ananth Ranganathan
• K. Madsen, H. B. Nielsen, O. Tingleff, Methods for Non-Linear Least Squares Problems (nonlinear least-squares tutorial; L-M code: analytic Jacobian secant)
• T. Strutz: Data Fitting and Uncertainty (A practical introduction to weighted least squares and beyond). 2nd edition, Springer Vieweg, 2016, テンプレート:ISBN2.
• H. P. Gavin, The Levenberg-Marquardt method for nonlinear least-squares curve-fitting problems（MATLAB内に実装されている。）