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

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

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

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

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

: \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つしかない場合、事前情報に基づかない均一な初期推定値、たとえば でも大域解に到達することができるが、複数の局所最小値が存在する場合、初期推定値が十分に大域最小点に近いときにしか大域解には収束しない。

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

: f\left (x_i, \boldsymbol\beta + \boldsymbol\delta\right ) \approx f\left (x _i, \boldsymbol\beta\right ) + \boldsymbol J_i \boldsymbol\delta

ここで、

: \boldsymbol J_i = \frac{\partial f\left (x_i, \boldsymbol\beta\right )}{\partial \boldsymbol\beta}

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

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

: 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

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

: \begin{align}
 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{align}

をに関して微分した結果を0とすると、以下の式を得る。

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

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

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

: \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)

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

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

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

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

: \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 ]

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

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

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

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

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

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

: 
\boldsymbol{v}_k + \frac{1}{2} \boldsymbol{a}_k


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

: 
\boldsymbol{J}_k \boldsymbol{a}_k = -f_{vv}(\boldsymbol x_k)


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

: 
\begin{align}
f_{vv}^i &\approx \frac{f_i(\boldsymbol{x} + h \boldsymbol{\delta}) - 2 f_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{align}


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

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

: 
\frac{2 \left\| \boldsymbol{a}_k \right\|}{\left\| \boldsymbol{v}_k \right\|} \le \alpha


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

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

## 例

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

## 脚注

## 関連文献

* 
* 
* 
* 

## 関連項目
* 信頼領域法
* ネルダー・ミード法

## 外部リンク
* 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, . Online copy
* History of the algorithm in SIAM news
* A tutorial by Ananth Ranganathan
* K. Madsen, H. B. Nielsen, O. Tingleff, *[http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf 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, .
* H. P. Gavin, ''The Levenberg-Marquardt method for nonlinear least-squares curve-fitting problems''（MATLAB内に実装されている。）