数理最適化の勉強メモ − Levenberg-Marquardt法

しばらく更新していない間にいいね/ブックマークがじわじわ増えててとても励みになります(*´ω`*)。最近は研究関連の勉強をしているので、ブログに書けないネタばかりで…。月1くらいで投稿できるようにしたいですね。

さて。このブログでは以前、数理最適化の勉強メモとして最適性条件, 最急降下法, ニュートン法についてまとめました。

今回の記事では問題設定を少し変えて、微分可能な制約なし非線形最小二乗問題の解法であるLevenberg-Marquardt法についてまとめてみます。

いつもの注意書きですが、筆者はコンピュータビジョンが専門で、数学や数値計算の専門家ではありません。間違いがある可能性も割と高いので、何か見つけたときは教えていただけると幸いです。

1. 非線形最小二乗問題とは

まずは非線形最小二乗問題の定義や性質から確認していきます。

1.1. 定式化

非線形最小二乗問題はその名の通り、複数の非線形関数 $r _ j(\mathbf{x})$ の二乗和を目的関数とする最適化問題です。

$f(\mathbf{x}) = \frac{1}{2} \sum _ {j=1} ^ m r _ j ^ 2(\mathbf{x})$

前回の記事では、一般的な微分可能な非線形関数 $f(\mathbf{x})$ を最小化することを考えましたが、今回は $f(\mathbf{x})$ がこのような特殊な形で表現できるケースだけに注目するわけです。

非線形最小二乗法は、主に関数フィッティングの目的で利用されます。たとえば、 $m$ 個のデータの組 $(\mathbf{u _ j}, \mathbf{v _ j})$ があって、これを $\mathbf{v} = g(\mathbf{x}, \mathbf{u})$ という関数で近似するとき、最も近似精度の良いパラメータ $\mathbf{x}$ を知るには、

$f(\mathbf{x}) = \frac{1}{2} \sum _ {j=1} ^ m \| g(\mathbf{x}, \mathbf{u} _ j) - \mathbf{v} _ j \| ^ 2$

を目的関数として最小二乗問題を解くのが一般的です。

このように、最小二乗法はデータを何らかのモデルに当てはめる問題に使えるため、様々な工学分野で利用されています。

1.2. 勾配ベクトルとヘッセ行列

前回の記事で紹介したように、最適化問題では勾配ベクトルとヘッセ行列が重要になるので、とりあえずこれらを計算してみます。

まず、 $r _ j(\mathbf{x})$ をベクトルの形にまとめて以下のように定義し直します。

$r(\mathbf{x}) = \begin{bmatrix}r _ 1(\mathbf{x}), \cdots, r _ m(\mathbf{x})\end{bmatrix} ^ \mathrm{T}$

すると目的関数は

$f(\mathbf{x})=\frac{1}{2} r(\mathbf{x}) ^ \mathrm{T} r(\mathbf{x})$

と表現できます。 $r(\mathbf{x})$ のヤコビアン $J(\mathbf{x}) \in \mathbb{R} ^ {m \times n}$ を計算すると以下のようになります。

$J(\mathbf{x}) = \begin{bmatrix} \nabla r _ 1(\mathbf{x}) ^ \mathrm{T} \\ \vdots \\ \nabla r _ m(\mathbf{x}) ^ \mathrm{T} \end{bmatrix} = \begin{bmatrix} \frac{\partial r _ 1}{\partial x _ 1} & \cdots & \frac{\partial r _ 1}{\partial x _ n} \\ \vdots & \ddots & \vdots \\ \frac{\partial r _ m}{\partial x _ 1} & \cdots & \frac{\partial r _ m}{\partial x _ n} \end{bmatrix}$

※ ちなみにプログラム上でどうやって微分値を計算するのか、という話題も以前の記事で解説しているのでよかったら読んでみてください→双対数を利用した自動微分のしくみ - かみのメモ。

このヤコビアンを利用すると、 $f(\mathbf{x})$ の勾配ベクトルとヘッセ行列は以下のように表されます。

$\nabla f(\mathbf{x}) = \sum ^ m _ {j=1} r _ j(\mathbf{x}) \nabla r _ j(\mathbf{x}) = J(\mathbf{x}) ^ \mathrm{T} r(\mathbf{x})$

$\begin{eqnarray}\nabla ^ 2 f(\mathbf{x}) &=& \sum ^ m _ {j=1} \nabla r _ j(\mathbf{x}) \nabla r _ j(\mathbf{x}) ^ \mathrm{T} + \sum ^ m _ {j=1} r _ j(x) \nabla ^ 2 r _ j(\mathbf{x}) \\ &=& J(\mathbf{x}) ^ \mathrm{T} J(\mathbf{x}) + \sum ^ m _ {j=1} r _ j(\mathbf{x}) \nabla ^ 2 r _ j(\mathbf{x}) \end{eqnarray}$

以降、 $\mathbf{J} _ k = J(\mathbf{x} _ k), \ \ \mathbf{r} _ k = r (\mathbf{x} _ k)$ と書くことにします。

1.3. 最小二乗問題のメリット

さて、ここでヘッセ行列の式に注目すると、以下のような考察ができます。

残差 $\mathbf{r} _ k$ のL2ノルムを小さくするような $\mathbf{x}$ を見つけることが最小二乗問題の目的なので、 $\mathbf{r} _ k$ の絶対値はそこそこ小さいはず
一般的に目的関数の凹凸はそこまで激しくないのではそこそこ小さいはず
- $f(\mathbf{x})$ の曲率変化が激しいときは大抵、関数の設定方法に難があるか、勾配法で扱うのが不適当な問題を解こうとしているとき

これらの仮定が成り立つとき、ヘッセ行列の2項目は十分小さくなるため無視でき、 $\nabla ^ 2 f(\mathbf{x} _ k) \approx \mathbf{J} _ k ^ \mathrm{T} \mathbf{J} _ k$ という近似がよく成り立つことがわかります。結果として、計算コストの高いヘッセ行列 $\nabla ^ 2 f(\mathbf{x})$ をヤコビアン $\mathbf{J} _ k$ の積で近似できることになります。

このおかげで、非線形最小二乗問題は一般的な非線形最適化問題より計算量的に解きやすいものになります。

※ ただし、変数 $\mathbf{x}$ の次元や目的関数の次元 $m$ が大きくなると（実装によるが 10⁵ ~ 10⁶ 以上）、ヤコビアンの次数も大きくなるので時間計算量・空間計算量ともに厳しくなってきます。そのため、画像データなどの次数が大きい問題の場合は、それが非線形最小二乗問題であっても準ニュートン法（ヘッセ行列を直接計算せず適当に推定する、非線形最小化問題一般に使える手法）を使うことが多いです。

2. ガウス・ニュートン法

非線形最小二乗問題を解くためのアルゴリズムの出発点がガウス・ニュートン法（Gauss-Newton method）です。

ガウス・ニュートン法は、前回の記事で紹介したニュートン法に上記の近似を適用したものです。

ニュートン法における探索方向は $-(\nabla ^ 2 f(\mathbf{x} _ k) ) ^ {-1} \nabla f(\mathbf{x} _ k)$ なので、ここに近似を適用すると $-(\mathbf{J} _ k ^ \mathrm{T} \mathbf{J} _ k) ^ {-1} \mathbf{J} _ k ^ \mathrm{T} \mathbf{r} _ k$ となります。この式を使って解を更新することで、ニュートン法の計算時間を削減できるわけです。

ただし、ガウス・ニュートン法はニュートン法と同じく、 $\mathbf{J} _ k ^ \mathrm{T} \mathbf{J} _ k$ が正定値でないときは非凸関数に近似されるため、減少方向でない更新ベクトルが選択されてしまうという問題点を抱えています。

たとえば、2変数関数 $f(x, y) = x ^ 4 -4x ^ 3 -14x ^ 2 + 12xy + 6y ^ 2$ （図中グリーン）の $(2, 4)$ （図中青点）における近似2次方程式は、鞍点 $(-0.4, 0.4)$ を停留点とする非凸関数（図中オレンジ）となります。よって、このままでは増加方向を探索することになってしまいます。

f:id:kamino-dev:20190917135557p:plain:w700 — 非凸関数に近似されてしまうケース

3. Levenberg-Marquardt法（オリジナル）

そんなガウス・ニュートン法の欠点を克服するために考案されたのがLevenberg-Marquardt（レーベンバーグ・マーカート/マルカート）法です。

LM法は文献ごとに、紹介されている内容にバラつきがあります。とりあえず、この章ではオリジナルのMarquardtの論文の記述をベースに紹介し、次の章でもう少し発展的な解釈について紹介していくことにします。

3.1. 基本的な考え方

LM法は、最急降下法やニュートン法と同じく、目的関数の勾配の情報を利用した反復法の一種です。各反復ステップでは、以下の更新式を使って解 $\mathbf{x} _ k$ を更新していきます。

$\mathbf{x} _ {k+1} = \mathbf{x} _ k -(\mathbf{J} _ k ^ \mathrm{T} \mathbf{J} _ k + \lambda \mathbf{I}) ^ {-1} \mathbf{J} _ k ^ \mathrm{T} \mathbf{r} _ k \ \ (\lambda \geq 0)$

直線探索を行っていた最急降下法やニュートン法と違い、計算されたベクトルをそのまま足し合わせている点に気をつけてください。つまりLM法では、パラメータ $\alpha _ k$ は使わず、更新方向と更新幅を同時に決定します。

それでは、この更新式の意味を考えていきましょう。

LM法のベースとなるアイデアは「凸関数に近似できたならガウス・ニュートン法を使い、そうでなければ最急降下法を使う」というものです。

上の式をよく見ると、 $\lambda = 0$ のときは更新方向が $-(\mathbf{J} _ k ^ \mathrm{T} \mathbf{J} _ k) ^ {-1} \mathbf{J} _ k ^ \mathrm{T} \mathbf{r} _ k$ でガウス・ニュートン法と同じになり、 $\lambda \gg 0$ のときは更新方向が $-\mathbf{J} _ k ^ \mathrm{T} \mathbf{r} _ k$ で最急降下法と同じになることがわかります。また、 $\lambda$ が大きくなるほど $(\mathbf{J} _ k ^ \mathrm{T} \mathbf{J} _ k + \lambda \mathbf{I}) ^ {-1}$ のトレースは小さくなるので、 $\mathbf{x} _ k$ の更新幅がどんどん小さくなることもわかります。

つまり、値が減少しない可能性はあるものの効率よく停留点に向かうガウス・ニュートン法と、必ず値が減少する最急降下法をうまく組み合わせているわけです。パラメータ $\lambda$ は2つの手法のバランスを取る役割を果たしています。

3.2. パラメータλの決め方

さて、肝心の $\lambda$ の決め方ですが、ここのアルゴリズムにはいくつかのバリエーションがあるようです。オリジナルのMarquardtの論文では、以下のような手法が提案されています。

まず、 $\Phi (\lambda) = f\left(\mathbf{x} _ k -(\mathbf{J} _ k ^ \mathrm{T} \mathbf{J} _ k + \lambda \mathbf{I}) ^ {-1} \mathbf{J} _ k ^ \mathrm{T} \mathbf{r} _ k\right)$ と定義します。つまり $\lambda$ を変えたとき $f(\mathbf{x} _ {k+1})$ の値がどう変わるのかに着目します。

そして、反復法の各ステップにおいて次の手順を実行します。

適当な初期値 $\lambda _ 0 \geq 0$ と、適当な倍率 $\nu > 1$ を決めておく
ならとする
- 値が順調に減少しているときはガウス・ニュートン法の成分を強める
かつならとする
- 現状維持
ならとなる最小のを使ってとする
- 値が減少するようになるまでガウス・ニュートン法の成分を弱める（ステップ幅が短くなり、かつ最急降下法の成分が強くなる）
$\lambda _ k$ から更新ベクトル $\mathbf{p} _ k$ を計算する

こうやって、目的関数値が減少する範囲でなるべく大きい $\lambda$ を選択するようにコントロールしていくわけです。

ここまでがオリジナルのLM法の考え方です。

4. 信頼領域法の観点からのLM法

ここからは、先ほど紹介したLM法に信頼領域法（trust region method）と呼ばれる考え方を絡めながら見ていきます。

Marquardtの論文で紹介されたLM法は「最急降下法とガウス・ニュートン法を組み合わせればいい感じになるじゃん」というある種ヒューリスティックなものでしたが（もちろん数学的に有効性が示された上でですが）、後の時代でさらに考察が進められたところ、信頼領域法という概念でも説明できることが示されました。

LM法から派生して生まれたこの信頼領域法は、最小二乗問題だけでなく一般の非線形最適化問題にも応用されています。

4.1. 信頼領域法とは

信頼領域法（Trust Region Method）は、ある意味で直線探索法と対になるような反復法の考え方です。

最急降下法やガウス・ニュートン法などの直線探索法では、「探索方向を決めてから更新幅を決める」という手順を取っていました。一方、信頼領域法では「 $f(\mathbf{x})$ をなにかの関数で近似したとき、その近似が十分に成り立っていると考えられる領域（信頼領域）を決めて、その中で近似式が最小になる点を選ぶ」という手法を取ります。

つまり、

目的関数を $\mathbf{x} _ k$ の周りで所定のモデル式に近似する
近似が十分に成り立っていると思われる領域内でモデル式を最小化する点を見つけて $\mathbf{x} _ {k+1}$ とする

という手順で解を更新するような反復解法が信頼領域法です。

4.2. 信頼領域法の基本的な考え方

信頼領域法自体は、最小二乗問題以外の一般的な非線形関数にも適用できるので、この節では $f(\mathbf{x})$ は微分可能な非線形関数全般であるものとして話を進めます。

状況設定としてはこれまでの反復法と同じで、現在の解 $\mathbf{x} _ k$ があったとき、より関数値を小さくするような更新ベクトル $\mathbf{p} _ k$ を求めることが目的です。

4.2.1. モデル関数への近似

まずは、目的関数を適当なモデル関数で近似します。多くの信頼領域法では、ニュートン法と同じく目的関数 $f(\mathbf{x})$ を二次多項式（テイラー展開の二次までの項）で近似します。

$f(\mathbf{x}) \approx f(\mathbf{x _ k}) + \nabla f(\mathbf{x} _ k)^\mathrm{T}(\mathbf{x} - \mathbf{x} _ k) + \frac{1}{2} (\mathbf{x} - \mathbf{x} _ k) ^ \mathrm{T} \nabla ^ 2 f(\mathbf{x} _ k)(\mathbf{x} - \mathbf{x} _ k)$

右辺を更新ベクトル $\mathbf{p}$ を使った式に書き直すと $\mathbf{x} = \mathbf{x} _ k + \mathbf{p}$ より、次のようになります。

$f(\mathbf{x}) \approx m _ k(\mathbf{p}) = f(\mathbf{x _ k}) + \nabla f(\mathbf{x _ k}) ^ \mathrm{T} \mathbf{p} + \frac{1}{2} \mathbf{p} ^ \mathrm{T} \nabla ^ 2 f(\mathbf{x} _ k) \mathbf{p}$

4.2.2. 信頼領域の定義

次に、 $m _ k(\mathbf{p})$ が $f(\mathbf{x})$ を十分に近似できていると思われる範囲を信頼領域として定義します。

$m _ k (\mathbf{p})$ は $\mathbf{x} _ k$ 周りのテイラー展開をベースにしているので、 $\mathbf{x} _ k$ 近傍では近似がよく成り立ち、そこから離れるほど近似精度が悪くなっていくと考えられます。そこで、 $\mathbf{x} _ k$ を中心とした球体 $\|\mathbf{p} _ k\| \leq \Delta _ k$ の領域を信頼領域とします（球体の代わりに楕円体を使うこともあります）。

半径 $\Delta _ k$ の決め方にはいくつかバリエーションがあるらしいのですが、今回参考にした書籍では以下のようなものが紹介されていました。

まず、適当な初期値 $\Delta _ 0$ と信頼領域の大きさの上限 $\hat{\Delta}$ を与えます。次に反復法の各ステップにおいて、とりあえず $\Delta = \Delta _ {k-1}$ として、後で紹介するKKT条件を満たす $\mathbf{p} ^ *$ を求めます。この $\mathbf{p} ^ *$ がちゃんと目的関数値を減少させるのかを検証するために、以下の式から $\rho$ を計算します。

$\rho = \frac{f(\mathbf{x} _ k + \mathbf{p} ^ *) - f(\mathbf{x} _ k)}{m _ k(\mathbf{p} ^ *) - m _ k(0)}$

モデル関数 $m _ k(\mathbf{p})$ が目的関数をよく近似しているなら、 $\rho$ は $1$ 前後の正の値を取るはずです。そこで以下のような操作を行います。

もし $\rho \lt 1/4$ なら信頼領域を $\Delta = \Delta _ {k-1}/4$ に縮小してもう一度 $\rho$ を計算してみる
もし $3/4 \lt \rho$ なら信頼領域を $\Delta = 2 \Delta _ {k-1}$ に拡大してもう一度 $\rho$ を計算してみる（ただし上限 $\hat{\Delta}$ に届いたら拡大をやめる）
$1/4 \leq\rho\leq 3/4$ もしくは $\Delta = \hat{\Delta}$ なら $\Delta _ k = \Delta$ 、 $\mathbf{p} _ k = \mathbf{p} ^ *$ と確定させる

このアルゴリズムで信頼領域の大きさをコントロールしていきます。

4.2.3. 部分問題の定式化

信頼領域が決まったところで、信頼領域の中で $m _ k (\mathbf{p})$ が最小になるような更新ベクトル $\mathbf{p} ^ *$ を探すことを考えます。

この問題は、以下のような制約付き非線形最適化問題として定式化できます。

$\min _ \mathbf{p} \ f(\mathbf{x _ k}) + \nabla f(\mathbf{x _ k}) ^ \mathrm{T} \mathbf{p} + \frac{1}{2} \mathbf{p} ^ \mathrm{T} \nabla ^ 2 f(\mathbf{x} _ k) \mathbf{p} \\\\ s.t. \ \ \|\mathbf{p}\| \leq \Delta _ k$

この部分問題を解けば、信頼領域の中でモデル関数が最小となる点までのベクトル、つまり最良の更新ベクトル $\mathbf{p} _ k$ が求まるわけです。

最適化問題を解くための反復ステップの中で別の最適化問題を解く、というのはなんだか変な感じがしますが、「効率的な更新ベクトルで解を更新し続ければ、より少ない反復回数で最適解に収束するだろう」という考えには納得がいくでしょう。

4.2.4. 部分問題の最適性条件

この部分問題の最適性条件（最適解であるための条件）について考えます。

今回は $\|\mathbf{p}\| \leq \Delta _ k$ という制約条件が入っているので、 $\nabla m _ k(\mathbf{p})$ がゼロベクトルになるかどうかだけでは最適解かどうかを判別できません。境界線上つまり $\|\mathbf{p}\| = \Delta _ k$ のときは、 $\nabla m _ k(\mathbf{p})$ がゼロベクトルでなくとも最適解になりうるからです。

そこで、ラグランジュの未定乗数法を使ってこの部分問題の最適性条件を探ります。この手順は、一般的な制約付き非線形最適化問題の解法と同じものです。スラック変数 $u \geq 0$ と未定乗数 $\lambda \geq 0$ を導入して、ラグランジュ関数を

$L(\mathbf{p}, \lambda) = f(\mathbf{x} _ k) + \nabla f(\mathbf{x} _ k) ^ \mathrm{T} \mathbf{p} + \frac{1}{2} \mathbf{p} ^ \mathrm{T} \nabla ^ 2 f(\mathbf{x} _ k) \mathbf{p} - \frac{\lambda}{2} (\Delta _ k ^ 2 - \mathbf{p} ^ \mathrm{T} \mathbf{p})$

と定義すると、KKT条件は

$(\nabla ^ 2f(\mathbf{x _ k}) + \lambda \mathbf{I})\mathbf{p} = - \nabla f(\mathbf{x _ k}) \\\\ \Delta _ k ^ 2 - \mathbf{p}^\mathrm{T} \mathbf{p} \geq 0 \\\\ \lambda(\Delta _ k ^ 2 - \mathbf{p}^\mathrm{T} \mathbf{p}) = 0 \\\\ \lambda \geq 0$

となります。

この条件をすべて満たすような $\lambda \geq 0$ と $\mathbf{p}$ の組み合わせが部分問題の最適解であり、それを求めることで最良の更新ベクトル $\mathbf{p} ^ *$ が得られます。 $\lambda, \mathbf{p}$ の見つけ方は、どの信頼領域法の手法を使うかによって変わってきます。

とりあえず、ここまでが信頼領域法の基本的な考え方です。

モデル関数に近似して、それがよく成り立つ信頼領域を設定した上で最良の更新ベクトルを求めようとすると、上のような部分問題が登場する、という流れが肝です。

4.3. LM法と信頼領域法の関係

さて、ここまでの説明からLM法と信頼領域法を関連付けることができます。

部分問題のKKT条件の1つ目の式は、 $(\nabla ^ 2f(\mathbf{x _ k}) + \lambda \mathbf{I})\mathbf{p} = - \nabla f(\mathbf{x _ k})$ でした。ここに最小二乗問題で使われる近似 $\nabla ^ 2 f(\mathbf{x} _ k) \approx \mathbf{J} _ k ^ \mathrm{T} \mathbf{J} _ k$ を適用すると、LM法で更新ベクトルを求めるのに使った式 $-(\mathbf{J} _ k ^ \mathrm{T} \mathbf{J} _ k + \lambda \mathbf{I}) ^ {-1} \mathbf{J} _ k ^ \mathrm{T} \mathbf{r} _ k$ と全く同じになることがわかります。

これがLM法と信頼領域法を紐付けているポイントです。つまり、LM法は最小二乗問題を二次多項式に近似して、半径 $\Delta _ k$ の球状の信頼領域を定義したときに出てくる部分問題を解いているのだと解釈できます。よって、LM法は信頼領域法に属するアルゴリズムであると見なせます。

ただ、先ほど説明したオリジナルのLM法では、直接的に信頼領域の半径 $\Delta _ k$ は登場せず、かわりに3.2. パラメータλの決め方で説明した $\lambda$ を決めるアルゴリズムが使われていました。

実は、このアルゴリズムが信頼領域の大きさを調整する作業と対応しています。 $f(\mathbf{x} _ k + \mathbf{p})$ が減少するような $\lambda _ k$ を求めて1つ目の条件式 $(\mathbf{J} _ k ^ \mathrm{T} \mathbf{J} _ k + \lambda _ k \mathbf{I})\mathbf{p} = - \mathbf{J} _ k ^ \mathrm{T} \mathbf{r} _ k$ から $\mathbf{p} _ k$ を求めたとします。すると、たしかに $f(\mathbf{x} + \mathbf{p} _ k)$ が減少しているので、後付けの解釈として信頼領域の大きさは妥当であったと言えます。
※ あくまで $m _ k(\mathbf{p})$ の値を見ながら $\lambda$ を調整しているだけなので、素早く収束するという意味で合理的な信頼領域になっている保証はありませんが…。

とまあこのように、ニュートン法と最急降下法を組み合わせただけに見えた更新ベクトルの式 $-(\mathbf{J} _ k ^ \mathrm{T} \mathbf{J} _ k + \lambda \mathbf{I}) ^ {-1} \mathbf{J} _ k ^ \mathrm{T} \mathbf{r} _ k$ は、実はある信頼領域の下での最適な更新ベクトルであったわけです。

4.4. 信頼領域法ベースのLM法

オリジナルのLM法では直接 $\lambda, \mathbf{p}$ を決定していましたが、先ほど紹介したように信頼領域 $\Delta _ k$ を決めてから部分問題を解くタイプの方法も存在します。最近は、むしろこちらの方法をLM法と呼ぶことが多いようですね。

先ほど導出した部分問題に最小二乗問題で使われる近似 $\nabla ^ 2 f(\mathbf{x} _ k) \approx \mathbf{J} _ k ^ \mathrm{T} \mathbf{J} _ k$ を当てはめると、与えられた信頼領域の半径 $\Delta _ k$ に対して

$(\mathbf{J} _ k ^ \mathrm{T} \mathbf{J} _ k + \lambda \mathbf{I})\mathbf{p} = - \mathbf{J} _ k ^ \mathrm{T} \mathbf{r} _ k \\\\ \Delta _ k ^ 2 - \|\mathbf{p} \| ^ 2 \geq 0 \\\\ \lambda(\Delta _ k ^ 2 - \|\mathbf{p} \| ^ 2) = 0 \\\\ \lambda \geq 0$

を満たす $\lambda, \mathbf{p}$ の組を求めればいいことになります。

ちなみに、楕円状の信頼領域を定義するときは、1つ目と2つ目の条件式が対角正定値行列 $\mathbf{D} _ k$ を使って

$(\mathbf{J} _ k ^ \mathrm{T} \mathbf{J} _ k + \lambda \mathbf{D} _ k ^ 2)\mathbf{p} = - \mathbf{J} _ k ^ \mathrm{T} \mathbf{r} _ k \\\\ \Delta _ k ^ 2 - \|\mathbf{D _ kp}\| ^ 2 \geq 0$

と変化します。

こうした部分問題の解き方にはいくつかのバリエーションがあるようなのですが、おおまかには以下のような考え方を利用するようです。

$\mathbf{p}(\lambda) = -(\mathbf{J} _ k ^ \mathrm{T} \mathbf{J} _ k + \lambda \mathbf{I}) ^ {-1} \mathbf{J} _ k ^ \mathrm{T} \mathbf{r} _ k$ を $\lambda$ の1変数関数に見立てる
$\Delta _ k ^ 2 = \| \mathbf{p}(\lambda) \|$ を満たす最大の $\lambda ^ *$ を探す（1次元のニュートン法）
$\lambda ^ * \geq 0$ なら $\mathbf{p}(\lambda ^ *)$ を、 $\lambda ^ * \lt 0$ なら $\mathbf{p}(0)$ を更新ベクトルとする

また、部分問題を近似的に解くバリエーションも存在します。あくまで目的は効率のよい更新ステップを見つけることなので、そこまで厳密な答えを探す必要はないだろうと割り切ってしまうわけです。

Googleが開発している最適化ライブラリCeres Solverには、こうしたいくつかのLM法アルゴリズムが実装されています（Solving Non-linear Least Squares — Ceres Solver）。どのアルゴリズムを採用するべきかはパラメータ数や $\mathbf{J} _ k ^ \mathrm{T} \mathbf{J} _ k$ の固有空間の様子に依存しますが、そこそこややこしい話になるので、実行してみてパフォーマンスのいいものを選ぶのが無難だと思います。

以上、LM法を直観的な観点と信頼領域法の観点から解説してみました。

このメモを書くに当たって以下の文献にお世話になりました。ありがとうございました。

Marquardt, Donald. An Algorithm for Least-Squares Estimation of Nonlinear Parameters. SIAM Journal on Applied Mathematics, 1963
NOCEDAL, Jorge; WRIGHT, Stephen. Numerical optimization Second Edition. Springer Science & Business Media, 2006.
Solving Non-linear Least Squares — Ceres Solver
矢部博. 工学基礎最適化とその応用. 数理工学社, 2006.

かみのメモ

コンピュータビジョン・プログラムな話題中心の勉強メモ