マルカート法 (Marquardt)法、あるいはレーベンバーグ・マルカート法（Levenberg–Marquardt, LM法) について説明します。この手法は最小二乗法の一種で、高速な収束速度と安定性を兼ね備えており、多くの場合十分な精度で最適値を求めることができるため、計算科学の分野で広く普及しています。なお、このアルゴリズムで得られる答えは、あくまで初期値近辺の局所最小値であり、大域最小値が得られる保証はないということにご注意ください。
\( n\) 個の観測点$$
(X_1, Y_1),(X_2, Y_2),\cdots,(X_n, Y_n)
$$ に対して\(m\) 個のパラメータを含む関数 (観測方程式) $$
F=F(a_1, a_2,\cdots, a_{m-1},a_m,X)
$$をフィッティングさせることを考えてみましょう。すなわち、観測値 \(Y\) をもっともよく説明できるような関数 \(F\) を探す (=パラメータ \(a\) を探す)ことがゴールです。ここで 「もっともよく説明できる」とは、観測点に関して観測値 \(Y\) と関数 \(F\) の残差の二乗を計算し、それらの和がなるべく小さくなることだとします。また各観測点の誤差を考慮する場合は、各観測点の重み (誤差のニ乗の逆数) を \(w_i\) として、残差の二乗に掛けておきます (誤差が分からない場合はとりあえず \(w_i=1\)として構いません)。すなわち、「なるべく小さくしたい量」 は以下のように定義される重み付き残差二乗の和 \(R\) です。$$
R=\sum_{i=1}^n{w_i[Y_i-F(a_1, a_2,\cdots, a_{m-1},a_m,X_i)]^2}
$$ この \(R\) を最小にするようなパラメータ \(a\) を見つけたいわけですが、そのための強力な方法の一つがマルカート法です。なお、関数 \(F\) はパラメータ \(a\) に対して微分可能であるとします1。 以下にマルカート法の手順を説明します。
STEP 1
まず適当な初期パラメータセット \(a^0 = (a^0_1,a^0_2,\cdots,a^0_{m-1},a^0_m) \)を与えてそのときの重み付き残差二乗の和 \(R\) を求めます。
$$R=\sum_{i=1}^n{w_i[Y_i-F(a^0,X)]^2}$$
STEP 2
次に\(F\)を各パラメータで偏微分し、そのときの各パラメータを入力値とするような成分を使って以下のような \(m \times m\) 対称行列 \(\alpha\) をつくります。$$
\alpha=
\begin{pmatrix}
(1+\lambda)\displaystyle\sum_{i=1}^n w_i \left(\frac{\partial F(a^0,X_i)}{\partial a_1}\right)^2 &
\displaystyle\sum_{i=1}^n w_i\frac{\partial F(a^0,X_i)}{\partial a_1}\frac{\partial F(a^0,X_i)}{\partial a_2} &
\ldots &
\displaystyle\sum_{i=1}^n w_i\frac{\partial F(a^0,X_i)}{\partial a_1}\frac{\partial F(a^0,X_i)}{\partial a_m}\\
\displaystyle\sum_{i=1}^n w_i\frac{\partial F(a^0,X_i)}{\partial a_2}\frac{\partial F(a^0,X_i)}{\partial a_1} &
(1+\lambda)\displaystyle\sum_{i=1}^n w_i \left(\frac{\partial F(a^0,X_i)}{\partial a_2}\right)^2 &
\ldots &
\displaystyle\sum_{i=1}^n w_i\frac{\partial F(a^0,X_i)}{\partial a_2}\frac{\partial F(a^0,X_i)}{\partial a_m}\\
\vdots & \vdots & \ddots & \vdots \\
\displaystyle\sum_{i=1}^n w_i\frac{\partial F(a^0,X_i)}{\partial a_m}\frac{\partial F(a^0,X_i)}{\partial a_1} &
\displaystyle\sum_{i=1}^n w_i\frac{\partial F(a^0,X_i)}{\partial a_m}\frac{\partial F(a^0,X_i)}{\partial a_2} &
\ldots &
(1+\lambda)\displaystyle\sum_{i=1}^n w_i \left( \frac{\partial F(a^0,X_i)}{\partial a_m} \right)^2
\end{pmatrix}
$$ここで \(\lambda\) は対角成分を強調させるための係数です。最初は \(\lambda=5\) くらいにしておくとよいでしょう。この行列は「誤差行列」と呼ばれています。というのも、関数 \(F\) が線形の場合、対角要素が各パラメータの分散と一致するからなのですが、詳しい説明は省略します。
STEP3
残差に偏微分をかけた\(m\) 行 1 列の行列\(\beta\)$$
\beta=
\begin{pmatrix}
\displaystyle\sum_{i=1}^n w_i{Y_i-F(a^0,X_i)}\frac{\partial F}{\partial a_1} \\
\displaystyle\sum_{i=1}^n w_i{Y_i-F(a^0,X_i)}\frac{\partial F}{\partial a_2} \\
\vdots \\
\displaystyle\sum_{i=1}^n w_i{Y_i-F(a^0,X_i)}\frac{\partial F}{\partial a_m}
\end{pmatrix}$$を計算します。
STEP4
\(\alpha\) と \(\beta\) を使って、新しいパラメータセット \((a’_1,a’_2,\cdots,a’_m)\) を以下のように求めます。$$
\begin{pmatrix} a’_1 \\ a’_2 \\ \vdots \\ a’_m \end{pmatrix}
=\begin{pmatrix} a^0_1 \\ a^0_2 \\ \vdots \\ a^0_m \end{pmatrix} + \alpha^{-1}\beta
$$
STEP 5
新しいパラメータセット \(a’\) を用いて重み付き残差二乗の和 \(R’\) を計算し、もし\(R'<R\) であれば (=残差が減少すれば) 新しいパラメータを初期値とし、\( \lambda\) を1/2 ~ 1/10 倍して、STEP 1に戻ります。一方、 もし\(R’ >= R\) であれば新しいパラメータを破棄して古いパラメータを再び初期値とし、 \( \lambda\) を2~10倍して、STEP 1に戻ります。
この作業を何回も繰り返して、 \(R\) の変化が十分に小さくなってきたところで計算を打ち切ります。
この手法の肝要な点は誤差行列の対角成分の重み付けを変化させながら次々と最適化を進めていく点です。\(R\) が順調に減っている場合では、自動的に \(\lambda\) の値は0に近づくため、対角項と非対角項はどんどん平等に扱われていきます。言い換えるとガウス・ニュートン法に近い振舞いとなります。
一方、\(R\) がなかなか減らない場合は、線形的な取り扱いができないため、ベクトル解析でいうところの \(\textrm{grad} R\) に沿ってパラメータを動かすほうが懸命です。これは \(\lambda\) の値を大きくして対角項を強調することと対応しています。言い換えると最急降下法に近い振る舞いとなります。
つまりこのマルカート法は、逐次計算に成功するとより大胆に(=線形的に)パラメータを動かし、失敗するとより慎重に勾配 (\(\textrm{grad}\)) に沿ってパラメータを動かすという手法です。失敗したらいったん引き返して慎重に道を選びなおすべし、というアルゴリズムは人生にも通じますね。
- マルカート法を適用するためには、関数 \(F\) はパラメータ \(a\) に対して非線形であってもかまいませんが、微分可能でなければいけません。 ↩︎