Levenberg-Marquardt算法

LM算法是一种可以用于求解最小二乘问题的迭代算法.它可以看成是最速下降法和Gauss-Newton法的结合(通过调节阻尼$\mu$切换).当当前解距离最优解较远时,算法更接近最速下降法:慢却保证下降;当当前解接近最优解,算法接近GN,快速收敛.
$f$为将参数向量$p\in {\mathbb{R}}^m$映射为估计观测向量$\hat{x}=f(P),\hat{x}\in {\mathbb{R}}^n$的函数.
输入: 初始的估计量$P_0$,观测向量$x$
寻找最优参数$P^{+}$
$$P^{+}=arg\underset{p}{\min }{ϵ}^Tϵ,ϵ=x-\hat{x},\hat{x}=f(P)$$在邻域泰勒展开$$f(P+{\delta }_P)\approx f(P)+J{\delta }_P$$其中$J$是雅克比矩阵$\frac{\partial f(P)}{\partial P}$.
寻找迭代每一步的步长${\delta }_P$使得最小化$||x-f(P+{\delta }_P)||\approx ||x-f(P)-J{\delta }_P||=||ϵ-J{\delta }_P||$,可以证明最小二乘的最优解存在于当$J{\delta }_P-ϵ$正交于$J$(展开使导数为零).
即$J^T(J{\delta }_P-ϵ)=0$,得到:
$$J^{T}J{\delta }_P=J^Tϵ$$此即为GN方法的增量正规方程.
LM在此基础上引入了阻尼项$\mu$,增量正规方程变为$$(J^{T}J+\mu I){\delta }_P=J^Tϵ$$,若当前求得的${\delta }_P$使得误差减小,则接受该更新且减小阻尼项$\mu$;反之,若当前增量使得函数值增大,则增大阻尼项,重新求解正规方程直到一个能使函数值减小的增量被求得.
LM算法每一步都将调整阻尼项$\mu$以确保误差下降.当$\mu$很大时,算法接近最速下降法,并且,步长会变小($\frac{1}{\mu }$),另外,还加强了$J^{T}J$的正定性;反之,接近GN.综上,LM是一种自适应的算法,当当前解远离最优解时慢速但下降,在最优解邻域快速收敛.
终止条件:

1. 梯度大小,如$J^Tϵ$低于阈值${\epsilon }_1$
2. 步长的变化量${\delta }_P$低于阈值${\epsilon }_2$
3. 达到最大迭代次数$k_{max}$

如果观测向量$x$的协方差矩阵${\Sigma }_x$可以获得,问题修正为$$J^T{\Sigma }_x^{-1}J{\delta }_P=J^T{\Sigma }_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