マルカール(Marquardt)法

マルカール(Marquardt)法について説明します。この手法は高速な収束速度と安定性を兼ね備えており、十分な精度で最適値を求めることができます。
マルカール法の要点を以下に述べます。 $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)$
をフィッティングさせることを考えます。すなわち最も観測値に合うパラメータ $a$ の組を探すことが目的です。関数 $F$ はパラメータ $a$ に対して非線形であってもかまいません。まず適当な初期パラメータ

$(a^0)=(a^0_1,a^0_2,\cdots,a^0_{m-1},a^0_m)$
を与えてそのときの二乗残差 $R$ を求めます。

$R=\sum_{i=1}^n{(Y_i-F(a^0,X))^2}$
次に $F$ を各パラメータで偏微分し、そのときの各パラメータを入力値とするような成分を使って以下のような $m$ 行 $m$ 列対称行列をつくります。

$\alpha= \left( \begin{array}{cccc} (1+\lambda)\displaystyle\sum_{i=1}^n w_i(\cfrac{\partial F(a^0,X_i)}{\partial a_1})^2 & \displaystyle\sum_{i=1}^n w_i\cfrac{\partial F(a^0,X_i)}{\partial a_1}\cfrac{\partial F(a^0,X_i)}{\partial a_2} & \ldots & \displaystyle\sum_{i=1}^n w_i\cfrac{\partial F(a^0,X_i)}{\partial a_1}\cfrac{\partial F(a^0,X_i)}{\partial a_m}\\ \displaystyle\sum_{i=1}^n w_i\cfrac{\partial F(a^0,X_i)}{\partial a_2}\cfrac{\partial F(a^0,X_i)}{\partial a_1} & (1+\lambda)\displaystyle\sum_{i=1}^n w_i(\cfrac{\partial F(a^0,X_i)}{\partial a_2})^2 & \ldots & \displaystyle\sum_{i=1}^n w_i\cfrac{\partial F(a^0,X_i)}{\partial a_2}\cfrac{\partial F(a^0,X_i)}{\partial a_m}\\ \displaystyle \vdots & \vdots & \ddots & \vdots \\ \displaystyle\sum_{i=1}^n w_i\cfrac{\partial F(a^0,X_i)}{\partial a_m}\cfrac{\partial F(a^0,X_i)}{\partial a_1} & \displaystyle\sum_{i=1}^n w_i\cfrac{\partial F(a^0,X_i)}{\partial a_m}\cfrac{\partial F(a^0,X_i)}{\partial a_2} & \ldots & (1+\lambda)\displaystyle\sum_{i=1}^n w_i(\cfrac{\partial F(a^0,X_i)}{\partial a_m})^2 \end{array} \right)$

ここで $\lambda$ は対角成分を強調させるための係数です。また $w_i$ は観測点ごとの重み(誤差の二乗の逆数)です。重みが分からない場合はすべて1でもかまいません。この行列は「誤差行列」と呼ばれています。その理由は線形の場合、対角要素がパラメータの分散になるからなのですが、詳しい説明は省略します。さらに残差に偏微分をかけたような $m$ 行 1 列の行列

$\beta= \left( \begin{array}{c} \displaystyle\sum_{i=1}^n w_i\{Y_i-F(a^0,X_i)\}\cfrac{\partial F}{\partial a_1} \\ \displaystyle\sum_{i=1}^n w_i\{Y_i-F(a^0,X_i)\}\cfrac{\partial F}{\partial a_2} \\ \vdots \\ \displaystyle\sum_{i=1}^n w_i\{Y_i-F(a^0,X_i)\}\cfrac{\partial F}{\partial a_m} \end{array} \right)$

を計算します。この二つの行列を使って新しいパラメータ

$a^0=(a'_1,a'_2,\cdots,a'_m)=(a')$

を以下のように求めます。

$\left(\begin{array}{c} a'_1 \\ a'_2 \\ \vdots \\ a'_m \end{array} \right) =\left(\begin{array}{c} a^0_1 \\ a^0_2 \\ \vdots \\ a^0_m \end{array} \right) + \alpha^{-1}\beta$

次にこの新しいパラメータを用いて二乗残差 $R'$ 　を計算します。
もし、 $R'<R$ であれば $a'$ を新しい初期値 $a_0$ とし、また $R'$ を新しい $R$ として採用し、さらに $\lambda$ を0.1~0.5倍します。もし、 $R' > R$ であれば新しいパラメータを破棄し、 $\lambda$ を2~10倍します。
そして再び、 $\alpha$ を計算するところに戻ります。この作業を何回も繰り返して、 $R$ の変化が適当に小さくなってきたところで計算を打ち切れば終了です。

この手法の肝要な点は誤差行列の対角成分の重み付けを変化させながら次々と最適化を進めていく点です。
$F$ が非線形な関数でも、初期値が真の値に近いときには $F$ のテイラー展開が一次で十分なため、線形的な取り扱いができます。このときは自動的にの値が0に近くなるように設定され、対角項と非対角項を平等に扱います。
一方初期値が真の値から遠いときは、線形的な取り扱いができないため、ベクトル解析でいうところの $grad R$ に沿ってパラメータを動かすほうが懸命です。これはの $\lambda$ は値を大きくして対角項を強調することを意味しています。
つまり、上で述べたアルゴリズムは逐次計算に成功するとより大胆に(=線形的に)パラメータを動かし、失敗するとより小さく勾配( $grad$ )に沿ってパラメータを動かす手法であるといえます。適切な初期値を選べばかなり高速に収束します。