[转载]LM算法的实现

完整文章请查看这里。转载请注明出处:本文来自learnhard的博客:http://www.codelast.com/ & http://blog.csdn.net/learnhard/,并保持文章的完整性。

 

LM算法可用于解决非线性最小二乘问题。多用于曲线拟合等场合。

LM算法的实现并不难,这里不讨论使用MATLAB等工具直接得到结果的过程,使用那些工具对于算法编程能力的提高无任何益处。 

LM算法的关键是用模型函数 f 对待估参数向量p在其领域内做线性近似,忽略掉二阶以上的导数项,从而转化为线性最小二乘问题,它具有收敛速度快等优点。

LM算法需要对每一个待估参数求偏导,所以,如果你的拟合函数 f 非常复杂,或者待估参数相当地多,那么就不适合使用LM算法了,可以使用Powell算法,Powell算法不需要求导。

对于急需自己编程(VC)用LM算法解决一些问题的朋友,如果你的数学几乎都忘了,那么你还是多请教一下自己的朋友吧,要不然连函数的偏导数都不记得怎么求了,是写不出代码的。

网上有很多LM算法的示例程序,但是如果你不理解这个算法的过程,要想看懂它们,很难。而且要对自己定义的函数使用LM算法,更加应该明白该算法的原理。

有一篇很不错的文章,解释了如何实现LM算法:http://www.ics.forth.gr/~lourakis/levmar/levmar.pdf

用Google搜索“Levenberg-Marquardt”,会有很多资料可参考。有一些现成的库也可以使用,不过,到你弄明白怎么用的时候,你都能够自己写出完整的代码了。当初我对LM也是很困惑,一直没弄清它的原理,网上的示例我怎么都用不对,后来一怒之下不再看网上的sample code,重新回到理论上,后来终于弄明白了,于是自己写出了完整的LM实现代码。

需要说明的是,这是非线性无约束的问题,如果待估参数是有约束的(例如参数在某一范围内变动),要想用在LM算法中,我还不知道怎样做,但是这一个帖子或许能给你一些启示(我尚未试验):http://www.numerical-recipes.com/forum/showthread.php?threadid=179

最后,不得不说的就是,LM算法并非许多人刚接触时想像的那般难,当你了解了过程之后,你就会觉得它很有意思。希望所有在学习它的朋友们都能成功。

转载自:https://blog.csdn.net/learnhard/article/details/1733843#commentBox

  • 1
    点赞
  • 3
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
以下是一个基本的Levenberg-Marquardt (LM)算法的MATLAB实现: ```matlab function [x, cost] = lm(fun, x0, options) % fun: 目标函数 % x0: 初始参数 % options: 优化选项(可选) % 默认优化选项 default_options = struct('Display', 'iter', 'MaxIter', 100, 'TolFun', 1e-6, 'TolX', 1e-6, 'Damping', 0.01); % 合并用户选项和默认选项 if nargin < 3 options = default_options; else options = merge_options(default_options, options); end % 初始化参数 x = x0; n = length(x); cost = zeros(options.MaxIter, 1); lambda = options.Damping; J = jac(fun, x); f = feval(fun, x); cost(1) = f'*f; % 迭代优化 for k = 2:options.MaxIter % 计算增量dx A = J'*J + lambda*eye(n); b = J'*f; dx = - A\b; % 计算新的参数值和目标函数值 x_new = x + dx; f_new = feval(fun, x_new); cost(k) = f_new'*f_new; % 更新参数和Jacobian矩阵 if cost(k) < cost(k-1) lambda = lambda / 10; x = x_new; f = f_new; J = jac(fun, x); if norm(dx) < options.TolX || norm(f) < options.TolFun break; end else lambda = lambda * 10; end % 显示优化信息 if strcmp(options.Display, 'iter') fprintf('Iteration %d, Cost = %f\n', k, cost(k)); end end % 截取优化结果 cost = cost(1:k-1); end function J = jac(fun, x) % 计算Jacobian矩阵 f = feval(fun, x); n = length(x); m = length(f); J = zeros(m, n); delta = 1e-8; for i = 1:n x1 = x; x1(i) = x1(i) + delta; f1 = feval(fun, x1); J(:,i) = (f1 - f)/delta; end end function options = merge_options(defaults, options) % 合并选项 fields = fieldnames(defaults); for i = 1:length(fields) if ~isfield(options, fields{i}) options.(fields{i}) = defaults.(fields{i}); end end end ``` 其中,`fun`是目标函数,`x0`是初始参数向量,`options`是优化选项。函数返回的`x`是优化后的参数向量,`cost`是每次迭代的目标函数值。该实现使用了一个简单的Jacobian计算方法,即使用有限差分法计算Jacobian矩阵。 需要注意的是,该实现并未进行任何异常处理,例如出现数值不稳定的情况。如果需要更加稳定的LM算法实现,建议使用专业的数值优化库。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值