数学 および計算科学 において、レーベンバーグ・マルカート法 （レーベンバーグ・マルカートほう、英 : Levenberg–Marquardt algorithm 、LM法 、レーベンバーグ・マーカート法）とは、非線形最小二乗 問題の解法の一つをいう。特に曲線回帰 を最小二乗法 により行う場合によく用いられる。LM法はガウス・ニュートン法 （GN法）と最急降下法 を内挿した手法といえる。GN法よりもロバスト であり、初期値が解から大きく外れていた場合でも解けることが多いが、ふるまいの良い関数に対してはGN法よりも収束が遅い傾向を示す。 LM法は、GN法に信頼領域法 を適用したものとみることもできる。
1944年 、フランクフォード・アーセナル （英語版 ） ケネス・レーヴェンバーグ （英語版 ） [ 1] 1963年 にデュポン 勤務の統計家 、ドナルド・マーカート （英語版 ） [ 2] [ 3] [ 4] [ 5]
LM法は、一般的な曲線回帰問題を解く必要のあるアプリケーションソフトウェア で広く用いられる。GN法を取り入れているため多くの場合で1次解法よりも収束速度が速いが[ 6] 反復法 と同様、LM法で保証されているのは局所最小値 への収束のみであり、大域最小値 が得られる保証はない。
LM法の主な適用対象は、最小二乗曲線回帰問題である。 m 対の独立変数 と従属変数 のペア
(
x
i
,
y
i
)
{\displaystyle \left(x_{i},y_{i}\right)}
β f (x , β 偏差 の二乗和S (β
β β -->
^ ^ -->
∈ ∈ -->
argmin
β β -->
-->
S
(
β β -->
)
≡ ≡ -->
argmin
β β -->
-->
∑ ∑ -->
i
=
1
m
[
y
i
− − -->
f
(
x
i
,
β β -->
)
]
2
{\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 (xi , β δ 線形近似 する。
f
(
x
i
,
β β -->
+
δ δ -->
)
≈ ≈ -->
f
(
x
i
,
β β -->
)
+
J
i
δ δ -->
{\displaystyle f\left(x_{i},{\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)\approx f\left(x_{i},{\boldsymbol {\beta }}\right)+{\boldsymbol {J}}_{i}{\boldsymbol {\delta }}}
ここで、
J
i
=
∂ ∂ -->
f
(
x
i
,
β β -->
)
∂ ∂ -->
β β -->
{\displaystyle {\boldsymbol {J}}_{i}={\frac {\partial f\left(x_{i},{\boldsymbol {\beta }}\right)}{\partial {\boldsymbol {\beta }}}}}
は関数f のβ 勾配 （ここでは行ベクトル とする）である。
偏差の二乗和S (β f (xi , β δ
S
(
β β -->
+
δ δ -->
)
≈ ≈ -->
∑ ∑ -->
i
=
1
m
[
y
i
− − -->
f
(
x
i
,
β β -->
)
− − -->
J
i
δ δ -->
]
2
{\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}}
ベクトル 表記すると、以下のように書ける。
S
(
β β -->
+
δ δ -->
)
≈ ≈ -->
‖
y
− − -->
f
(
β β -->
)
− − -->
J
δ δ -->
‖
2
=
[
y
− − -->
f
(
β β -->
)
− − -->
J
δ δ -->
]
T
[
y
− − -->
f
(
β β -->
)
− − -->
J
δ δ -->
]
=
[
y
− − -->
f
(
β β -->
)
]
T
[
y
− − -->
f
(
β β -->
)
]
− − -->
[
y
− − -->
f
(
β β -->
)
]
T
J
δ δ -->
− − -->
(
J
δ δ -->
)
T
[
y
− − -->
f
(
β β -->
)
]
+
δ δ -->
T
J
T
J
δ δ -->
=
[
y
− − -->
f
(
β β -->
)
]
T
[
y
− − -->
f
(
β β -->
)
]
− − -->
2
[
y
− − -->
f
(
β β -->
)
]
T
J
δ δ -->
+
δ δ -->
T
J
T
J
δ δ -->
{\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 (β δ δ
(
J
T
J
)
δ δ -->
=
J
T
[
y
− − -->
f
(
β β -->
)
]
{\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 (xi ), yi 正方行列 ではなく、m をデータ点数、n をベクトルβ m × n J T J n × n n 連立線形方程式であるからこれを解いてδ
LM法では、この方程式を次のように「減衰」させたものに置き換える。
(
J
T
J
+
λ λ -->
I
)
δ δ -->
=
J
T
[
y
− − -->
f
(
β β -->
)
]
(
λ λ -->
≥ ≥ -->
0
)
{\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 であることに注意すると、λ が大きいときδ δ β δ β
減衰係数λ が‖ J T J に比べて大きいときは、J T J I δ
λ λ -->
− − -->
1
J
T
[
y
− − -->
f
(
β β -->
)
]
{\displaystyle \lambda ^{-1}{\boldsymbol {J}}^{\mathrm {T} }\left[{\boldsymbol {y}}-{\boldsymbol {f}}\left({\boldsymbol {\beta }}\right)\right]}
LM法は、λ が大きい値のときJ T J 1971年 、勾配が小さい方向への収束が遅いという問題を避けるため、勾配を曲率に応じてスケールするという考えから、単位行列I J T J 対角要素 で置き換え、解をスケール不変にする手法が提案された[ 7]
[
J
T
J
+
λ λ -->
diag
-->
(
J
T
J
)
]
δ δ -->
=
J
T
[
y
− − -->
f
(
β β -->
)
]
{\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]}
同様の減衰因子は、線形不良設定問題 を解くために用いられるティホノフ正則化 （英語版 ） リッジ回帰 と呼ばれる推計 法にも現れる。
最良の減衰係数λ を選ぶ方法としては、様々な議論があるが、大なり小なりヒューリスティック なものである。それらの選び方がなぜ局所最小点への収束を保証するかを示す理論的な議論はあるが、大域最小点へ収束するような選び方をすると最急降下法 の望ましくない特質、とくに収束が遅いという側面が表われてしまう。
どんな選び方をしても、パラメータの大きさはもとの問題がどれほど良くスケールするかに依存する。マーカートは次のような選び方を推奨している。まず初期値λ = λ 0 S (β )ν > 1λ = λ 0 / ν ν を掛けλ 0 νk
減衰係数λ /ν λ の新しい値とし（かつλ /ν λ /ν λ を用いた残差よりも大きくなったならば、λ の値を変えず、λ を用いた結果を採用する。
delayed gratification [訳語疑問点 と呼ばれる減衰係数の効果的な制御方法がある。この方法では、上り坂のステップごとに係数を少しずつ増やし、下り坂のステップごとにパラメーターを大幅に減らす。この戦略は、最適化の開始時に坂を下りすぎ、後に使用できるステップが制限されて、収束が遅くなることを防ぐことを主眼においている[ 8] [ 9]
レーベンバーグ・マルカート法の更新ステップv k 測地経路 に沿った速度と捉えると、測地経路に沿う加速度に対応する2次の項a k
v
k
+
1
2
a
k
{\displaystyle {\boldsymbol {v}}_{k}+{\frac {1}{2}}{\boldsymbol {a}}_{k}}
ここで、a k
J
k
a
k
=
− − -->
f
v
v
(
x
k
)
{\displaystyle {\boldsymbol {J}}_{k}{\boldsymbol {a}}_{k}=-f_{vv}({\boldsymbol {x}}_{k})}
この測地線加速度項は速度v 方向微分
f
v
v
(
x
)
=
∑ ∑ -->
μ μ -->
ν ν -->
v
μ μ -->
v
ν ν -->
∂ ∂ -->
μ μ -->
∂ ∂ -->
ν ν -->
f
(
x
)
{\displaystyle f_{vv}({\boldsymbol {x}})=\sum _{\mu \nu }v_{\mu }v_{\nu }\partial _{\mu }\partial _{\nu }f({\boldsymbol {x}})}
二次導関数行列 を計算する必要はなく、計算コスト上のオーバーヘッドは比較的小さい[ 10] 有限差分 近似に置き換えると便利な場合がある。
f
v
v
i
≈ ≈ -->
f
i
(
x
+
h
δ δ -->
)
− − -->
2
f
i
(
x
)
+
f
i
(
x
− − -->
h
δ δ -->
)
h
2
=
2
h
(
f
i
(
x
+
h
δ δ -->
)
− − -->
f
i
(
x
)
h
− − -->
J
i
δ δ -->
)
{\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 の選択によってはアルゴリズムが不安定になる場合があるが、通常はおよそ0.1が妥当である[ 9]
加速度は速度と反対の方向を指す可能性があり、減衰が小さすぎて収束が失速するのを防ぐために、ステップの採用条件に加速度に関する以下のような追加の基準が追加される。
2
‖
a
k
‖
‖
v
k
‖
≤ ≤ -->
α α -->
{\displaystyle {\frac {2\left\|{\boldsymbol {a}}_{k}\right\|}{\left\|{\boldsymbol {v}}_{k}\right\|}}\leq \alpha }
ここでα は通常は1未満の固定値とされる。より難しい問題ではより小さい値とする[ 9]
測地線加速度項を追加すると、収束速度が大幅に向上する可能性があり、特にアルゴリズムが目的関数ランドスケープ上の狭い峡谷を移動する場合に有用である。このような場合、2次の項の効果によりステップ幅はより小さく、精度が高くなる[ 9]
悪いフィッティング結果
ましなフィッティング結果
良いフィッティング結果
この例では、関数
y
=
a
cos
-->
(
b
X
)
+
b
sin
-->
(
a
X
)
{\displaystyle y=a\cos \left(bX\right)+b\sin \left(aX\right)}
GNU Octave 上のレーベンバーグ・マルカート法実装leasqr 関数をもちいてフィッティングする。グラフから、パラメーターフィッティング結果が徐々に改善し、
a
=
100
{\displaystyle a=100}
b
=
102
{\displaystyle b=102}
cos(βx ) に一致するパラメータの値はˆ β ˆ β nπ
^ Levenberg, Kenneth (1944). “A method for the solution of certain non-linear problems in least squares” (英語). Quarterly of Applied Mathematics 2 (2): 164–168. doi :10.1090/qam/10666 . ISSN 0033-569X . https://www.ams.org/qam/1944-02-02/S0033-569X-1944-10666-0/ . ^ Marquardt, Donald (1963). “An Algorithm for Least-Squares Estimation of Nonlinear Parameters”. SIAM Journal on Applied Mathematics 11 (2): 431–441. doi :10.1137/0111030 . hdl :10338.dmlcz/104299 ^ Girard, André (1958). “Excerpt from Revue d'optique théorique et instrumentale ”. Rev. Opt. 37 : 225–241, 397–424. ^ Wynne, C. G. (1959). “Lens Designing by Electronic Digital Computer: I”. Proc. Phys. Soc. Lond. 73 (5): 777–787. Bibcode : 1959PPS....73..777W . doi :10.1088/0370-1328/73/5/310 . ^ Morrison, David D. (1960). “Methods for nonlinear least squares problems and convergence proofs”. Proceedings of the Jet Propulsion Laboratory Seminar on Tracking Programs and Orbit Determination : 1–9. ^ Wiliamowski, Bogdan; Yu, Hao (June 2010). “Improved Computation for Levenberg–Marquardt Training” . IEEE Transactions on Neural Networks and Learning Systems 21 (6). https://www.eng.auburn.edu/~wilambm/pap/2010/Improved%20Computation%20for%20LM%20Training.pdf . ^ FLETCHER, R. (1971). “A Modified Marquardt Subroutine for Nonlinear Least Squares”. Harwell Report AERE-R (6799). NAID 10000008775 . ^ Transtrum, Mark K; Machta, Benjamin B; Sethna, James P (2011). “Geometry of nonlinear least squares with applications to sloppy models and optimization”. Physical Review E (APS) 83 (3): 036701. arXiv :1010.1449 . Bibcode : 2011PhRvE..83c6701T . doi :10.1103/PhysRevE.83.036701 . PMID 21517619 . ^ a b c d Transtrum, Mark K; Sethna, James P (2012). "Improvements to the Levenberg-Marquardt algorithm for nonlinear least-squares minimization". arXiv :1201.5885 physics.data-an ]。
^ “Nonlinear Least-Squares Fitting ”. GNU Scientific Library. 2020年4月14日時点のオリジナル よりアーカイブ。2022年9月17日 閲覧。
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, ISBN 0-89871-433-8 . Online copy
History of the algorithm in SIAM news A tutorial by Ananth Ranganathan K. Madsen, H. B. Nielsen, O. Tingleff, Methods for Non-Linear Least Squares Problems analytic Jacobian secant )
T. Strutz: Data Fitting and Uncertainty (A practical introduction to weighted least squares and beyond). 2nd edition, Springer Vieweg, 2016, ISBN 978-3-658-11455-8 .
H. P. Gavin, The Levenberg-Marquardt method for nonlinear least-squares curve-fitting problems MATLAB implementation included)
^ Kanzow, Christian; Yamashita, Nobuo; Fukushima, Masao (2004). “Levenberg–Marquardt methods with strong local convergence properties for solving nonlinear equations with convex constraints”. Journal of Computational and Applied Mathematics 172 (2): 375–397. Bibcode : 2004JCoAM.172..375K . doi :10.1016/j.cam.2004.02.013 .
![{\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}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/214db422df83e90332f678b04269388fd3162cbc)
![{\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}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/cbe6c85ed396e6f74ef38478d2b78c8d31958061)
![{\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/025fc5098589de6fcfcddd448290cd730e38cf31)
![{\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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9761bef41e277b5d395438e46cb6e57a3697b148)
![{\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)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6b85a9d81ac1dc673d1d0e82e38065cd76281e92)
![{\displaystyle \lambda ^{-1}{\boldsymbol {J}}^{\mathrm {T} }\left[{\boldsymbol {y}}-{\boldsymbol {f}}\left({\boldsymbol {\beta }}\right)\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e9b992046a9b0bc7df0c22bb794f92769dab9b50)
![{\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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b225ab30915bd00438366830a5aa9e7a3b84c6d1)
