LM算法是一种可以用于求解最小二乘问题的迭代算法.它可以看成是最速下降法和Gauss-Newton法的结合(通过调节阻尼
μ
\mu
μ切换).当当前解距离最优解较远时,算法更接近最速下降法:慢却保证下降;当当前解接近最优解,算法接近GN,快速收敛.
f
f
f为将参数向量
p
∈
R
m
p\in\Bbb R^m
p∈Rm映射为估计观测向量
x
^
=
f
(
P
)
,
x
^
∈
R
n
\hat{x}=f(P),\hat{x}\in\Bbb R^n
x^=f(P),x^∈Rn的函数.
输入: 初始的估计量
P
0
P_0
P0,观测向量
x
x
x
寻找最优参数
P
+
P^+
P+
P
+
=
a
r
g
min
p
ϵ
T
ϵ
,
ϵ
=
x
−
x
^
,
x
^
=
f
(
P
)
P^+=arg\min_{p} \epsilon^T\epsilon,\epsilon=x-\hat{x},\hat{x}=f(P)
P+=argpminϵTϵ,ϵ=x−x^,x^=f(P)在邻域泰勒展开
f
(
P
+
δ
P
)
≈
f
(
P
)
+
J
δ
P
f(P+\delta_P)\approx f(P)+J\delta_P
f(P+δP)≈f(P)+JδP其中
J
J
J是雅克比矩阵
∂
f
(
P
)
∂
P
\frac{\partial f(P)}{\partial P}
∂P∂f(P).
寻找迭代每一步的步长
δ
P
\delta_P
δP使得最小化
∣
∣
x
−
f
(
P
+
δ
P
)
∣
∣
≈
∣
∣
x
−
f
(
P
)
−
J
δ
P
∣
∣
=
∣
∣
ϵ
−
J
δ
P
∣
∣
||x-f(P+\delta_P)||\approx||x-f(P)-J\delta_P||=||\epsilon-J\delta_P||
∣∣x−f(P+δP)∣∣≈∣∣x−f(P)−JδP∣∣=∣∣ϵ−JδP∣∣,可以证明最小二乘的最优解存在于当
J
δ
P
−
ϵ
J\delta_P-\epsilon
JδP−ϵ正交于
J
J
J(展开使导数为零).
即
J
T
(
J
δ
P
−
ϵ
)
=
0
J^T(J\delta_P-\epsilon)=0
JT(JδP−ϵ)=0,得到:
J
T
J
δ
P
=
J
T
ϵ
J^TJ\delta_P=J^T\epsilon
JTJδP=JTϵ此即为GN方法的增量正规方程.
LM在此基础上引入了阻尼项
μ
\mu
μ,增量正规方程变为
(
J
T
J
+
μ
I
)
δ
P
=
J
T
ϵ
(J^TJ+\mu I)\delta_P=J^T\epsilon
(JTJ+μI)δP=JTϵ,若当前求得的
δ
P
\delta_P
δP使得误差减小,则接受该更新且减小阻尼项
μ
\mu
μ;反之,若当前增量使得函数值增大,则增大阻尼项,重新求解正规方程直到一个能使函数值减小的增量被求得.
LM算法每一步都将调整阻尼项
μ
\mu
μ以确保误差下降.当
μ
\mu
μ很大时,算法接近最速下降法,并且,步长会变小(
1
μ
\frac{1}{\mu}
μ1),另外,还加强了
J
T
J
J^TJ
JTJ的正定性;反之,接近GN.综上,LM是一种自适应的算法,当当前解远离最优解时慢速但下降,在最优解邻域快速收敛.
终止条件:
- 梯度大小,如 J T ϵ J^T\epsilon JTϵ低于阈值 ε 1 \varepsilon_1 ε1
- 步长的变化量 δ P \delta_P δP低于阈值 ε 2 \varepsilon_2 ε2
- 达到最大迭代次数 k m a x k_{max} kmax
如果观测向量 x x x的协方差矩阵 Σ x \Sigma_x Σx可以获得,问题修正为 J T Σ x − 1 J δ P = J T Σ X − 1 ϵ J^T\Sigma_x^{-1}J\delta_P=J^T\Sigma_X^{-1}\epsilon JTΣx−1JδP=JTΣX−1ϵ该式为加权正规方程.
算法
kmax=100;
eta1=eta2=1e-15;
tau=1e-3;
k:=0;v:=2;p:=p0;
A:=J.T*J;epsilon:=x-f(p);g:=J.T*epsilon;
stop:=(max(g)<eta1);
mu=tau*max(diag(A));
while ( !stop ) and ( k<kmax ):
k:=k+1;
while ( !stop ) and ( rho<0 ):
Solve (A+mu*I)*delta=g;
if(delta.norm()<=eta2*p.norm())
stop:=true;
else
pnew:=p+delta;
rho:=(epsilon.norm()-(x-f(pnew)).norm())/(delta.T*(mu*delta+g));
if rho>0
p:=pnew;
A:=J.T*J;epsilon:=x-f(p);g:=J.T*epsilon;
stop:=(max(g)<eta1);
mu:=mu*max(1/3,1-(2*rho-1)^3);v:=2;
else
mu:=mu*v;v:=2*v;
endif
endif
endwhile
endwhile