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

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

1944年、フランクフォード・アーセナル（英語版）職員ケネス・レーヴェンバーグ（英語版）により初発表され^[1]、1963年にデュポン勤務の統計家、ドナルド・マーカート（英語版）により再発見された^[2]。また、Girard^[3], Wynne^[4], Morrison^[5] によりそれぞれ独立に再発見されている。

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

問題設定

LM法の主な適用対象は、最小二乗曲線回帰問題である。 m対の独立変数と従属変数のペア${\displaystyle \left(x_{i},y_{i}\right)}$が与えられたとき、βをパラメータとする曲線f(x, β)と所与のペアとの偏差の二乗和S(β)を最小化することを考える。

${\displaystyle {\hat {\boldsymbol {\beta }}}\in \operatorname {argmin} \limits _{\boldsymbol {\beta }}S\left({\boldsymbol {\beta }}\right)\equiv \operatorname {argmin} \limits _{\boldsymbol {\beta }}\sum _{i=1}^{m}\left[y_{i}-f\left(x_{i},{\boldsymbol {\beta }}\right)\right]^{2}}$

解法

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

各反復ステップにおいてパラメーターベクトルβは新しい推定値β + δへ置き換えられる。δを決めるため、f(x_i, β + δ)を次のように線形近似する。

${\displaystyle f\left(x_{i},{\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)\approx f\left(x_{i},{\boldsymbol {\beta }}\right)+{\boldsymbol {J}}_{i}{\boldsymbol {\delta }}}$

ここで、

${\displaystyle {\boldsymbol {J}}_{i}={\frac {\partial f\left(x_{i},{\boldsymbol {\beta }}\right)}{\partial {\boldsymbol {\beta }}}}}$

は関数fのβについての勾配（ここでは行ベクトルとする）である。

偏差の二乗和S(β)は、この勾配がゼロのとき局所最小となる。上式の一次近似を用いると、f(x_i, β + δ)の偏差二乗和は以下のように近似される。

${\displaystyle S\left({\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)\approx \sum _{i=1}^{m}\left[y_{i}-f\left(x_{i},{\boldsymbol {\beta }}\right)-{\boldsymbol {J}}_{i}{\boldsymbol {\delta }}\right]^{2}}$

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

${\displaystyle {\begin{aligned}S\left({\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)&\approx \left\|{\boldsymbol {y}}-{\boldsymbol {f}}\left({\boldsymbol {\beta }}\right)-{\boldsymbol {J}}{\boldsymbol {\delta }}\right\|^{2}\\&=\left[{\boldsymbol {y}}-{\boldsymbol {f}}\left({\boldsymbol {\beta }}\right)-{\boldsymbol {J}}{\boldsymbol {\delta }}\right]^{\mathrm {T} }\left[{\boldsymbol {y}}-{\boldsymbol {f}}\left({\boldsymbol {\beta }}\right)-{\boldsymbol {J}}{\boldsymbol {\delta }}\right]\\&=\left[{\boldsymbol {y}}-{\boldsymbol {f}}\left({\boldsymbol {\beta }}\right)\right]^{\mathrm {T} }\left[{\boldsymbol {y}}-{\boldsymbol {f}}\left({\boldsymbol {\beta }}\right)\right]-\left[{\boldsymbol {y}}-{\boldsymbol {f}}\left({\boldsymbol {\beta }}\right)\right]^{\mathrm {T} }{\boldsymbol {J}}{\boldsymbol {\delta }}-\left({\boldsymbol {J}}{\boldsymbol {\delta }}\right)^{\mathrm {T} }\left[{\boldsymbol {y}}-{\boldsymbol {f}}\left({\boldsymbol {\beta }}\right)\right]+{\boldsymbol {\delta }}^{\mathrm {T} }{\boldsymbol {J}}^{\mathrm {T} }{\boldsymbol {J}}{\boldsymbol {\delta }}\\&=\left[{\boldsymbol {y}}-{\boldsymbol {f}}\left({\boldsymbol {\beta }}\right)\right]^{\mathrm {T} }\left[{\boldsymbol {y}}-{\boldsymbol {f}}\left({\boldsymbol {\beta }}\right)\right]-2\left[{\boldsymbol {y}}-{\boldsymbol {f}}\left({\boldsymbol {\beta }}\right)\right]^{\mathrm {T} }{\boldsymbol {J}}{\boldsymbol {\delta }}+{\boldsymbol {\delta }}^{\mathrm {T} }{\boldsymbol {J}}^{\mathrm {T} }{\boldsymbol {J}}{\boldsymbol {\delta }}\end{aligned}}}$

S(β + δ)をδに関して微分した結果を0とすると、以下の式を得る。

${\displaystyle \left({\boldsymbol {J}}^{\mathrm {T} }{\boldsymbol {J}}\right){\boldsymbol {\delta }}={\boldsymbol {J}}^{\mathrm {T} }\left[{\boldsymbol {y}}-{\boldsymbol {f}}\left({\boldsymbol {\beta }}\right)\right]}$

ここで、Jはヤコビ行列であり、そのi行目はJ_i に等しい。また、f(β), yはそれぞれ、i行目成分をf(x_i), y_iとするベクトルである。ヤコビ行列は一般的には正方行列ではなく、mをデータ点数、nをベクトルβのサイズとしてm × n長方形行列である。行列積J^TJはn × n正方行列となり、上式はn連立線形方程式であるからこれを解いてδを得ることができる。これをそのまま解くのがガウス・ニュートン法である。

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

${\displaystyle \left({\boldsymbol {J}}^{\mathrm {T} }{\boldsymbol {J}}+\lambda \mathbf {I} \right){\boldsymbol {\delta }}={\boldsymbol {J}}^{\mathrm {T} }\left[{\boldsymbol {y}}-{\boldsymbol {f}}\left({\boldsymbol {\beta }}\right)\right]\quad (\lambda \geq 0)}$

ここで、Iは単位行列である。これを解いて得られるδを用いてパラメータベクトルβの推定値を更新する。

非負の減衰係数λは各反復ごとに調整される。Sが急速に減少する際には小さい値が用いられ、LM法はGN法に近づく。対して、残差が十分に減少しない場合は大きい値のλが用いられ、Sのβについての勾配は−2 (J^T[y-f(β)])^Tであることに注意すると、λが大きいときδは勾配の逆向きに近付きLM法は最急降下法に近づくことがわかる。計算されたδが十分小さくなったとき、もしくは得られたパラメータ推定値β + δに置き換えた際の偏差二乗和の減少が十分に小さくなったときのどちらかの場合に反復は打ち切られ、解βを得る。

減衰係数λが‖ J^TJ ‖に比べて大きいときは、J^TJ+λIの逆行列を求める必要はなく、更新ステップδは${\displaystyle \lambda ^{-1}{\boldsymbol {J}}^{\mathrm {T} }\left[{\boldsymbol {y}}-{\boldsymbol {f}}\left({\boldsymbol {\beta }}\right)\right]}$で十分に近似される。

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

${\displaystyle \left[{\boldsymbol {J}}^{\mathrm {T} }{\boldsymbol {J}}+\lambda \operatorname {diag} \left({\boldsymbol {J}}^{\mathrm {T} }{\boldsymbol {J}}\right)\right]{\boldsymbol {\delta }}={\boldsymbol {J}}^{\mathrm {T} }\left[{\boldsymbol {y}}-{\boldsymbol {f}}\left({\boldsymbol {\beta }}\right)\right]}$

同様の減衰因子は、線形不良設定問題を解くために用いられるティホノフ正則化（英語版）や、リッジ回帰と呼ばれる推計法にも現れる。

減衰パラメータの選び方

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

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

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

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

測地線加速度項

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

${\displaystyle {\boldsymbol {v}}_{k}+{\frac {1}{2}}{\boldsymbol {a}}_{k}}$

ここで、a_kは次の方程式の解である。

${\displaystyle {\boldsymbol {J}}_{k}{\boldsymbol {a}}_{k}=-f_{vv}({\boldsymbol {x}}_{k})}$

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

${\displaystyle {\begin{aligned}f_{vv}^{i}&\approx {\frac {f_{i}({\boldsymbol {x}}+h{\boldsymbol {\delta }})-2f_{i}({\boldsymbol {x}})+f_{i}({\boldsymbol {x}}-h{\boldsymbol {\delta }})}{h^{2}}}\\&={\frac {2}{h}}\left({\frac {f_{i}({\boldsymbol {x}}+h{\boldsymbol {\delta }})-f_{i}({\boldsymbol {x}})}{h}}-{\boldsymbol {J}}_{i}{\boldsymbol {\delta }}\right)\end{aligned}}}$

ここで、f(x)とJは通常のアルゴリズムでもすでに計算されているため、追加で計算する必要があるのはf(x+hδ)だけである。有限差分ステップhの選択によってはアルゴリズムが不安定になる場合があるが、通常はおよそ0.1が妥当である^[9]。

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

${\displaystyle {\frac {2\left\|{\boldsymbol {a}}_{k}\right\|}{\left\|{\boldsymbol {v}}_{k}\right\|}}\leq \alpha }$

ここでαは通常は1未満の固定値とされる。より難しい問題ではより小さい値とする^[9]。

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

例

• 
  悪いフィッティング結果
• 
  ましなフィッティング結果
• 
  良いフィッティング結果

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