首页 > 代码库 > Levenberg-Marquardt 的 MATLAB 代码
Levenberg-Marquardt 的 MATLAB 代码
参考资料:
1,《精通MATLAB最优化计算(第2版)》作者:龚纯 等 的 第9章 9.3 小节 L-M 法
2,《数值分析》 作者:Timothy Sauer 的 第4章 4.4节 非线性最小二乘的 例子
第一本书里头虽然有代码,然而有错误,修正了错误之处
% opti_LM_test1% 测试了 MATLAB最优化 书中的 L-M 的例子,结果是正确的clear all;clc;close all;syms t;f = ... [t^2+t-1; 2*t^2-3];S = transpose(f)*f;f_var = symvar(f);t_init = -5 % 自变量的初始值%%u = 2v = 1.5beta = 0.4eps = 1.0e-6x = t_init;x = transpose(x);% 删jacobian_f = jacobian(f,f_var);tol = 1;%% subs以后居然不是数值,而是符号!还要转换成double类型!!!while tol>eps fxk = double(subs(f,f_var,x)); Sxk = double(subs(S,f_var,x)); % step2: 计算 fxk Sxk delta_fxk = double(subs(jacobian_f,f_var,x)); % step3: 计算 delta_fxk delta_Sxk = transpose(delta_fxk)*fxk; % step4: 计算 delta_Sxk while 1 % step5: 计算Q,并解方程(Q+uI)delta_x = -delta_Sxk Q = transpose(delta_fxk)*delta_fxk; dx = -(Q+u*eye(size(Q)))\delta_Sxk; x1 = x + dx; fxk = double(subs(f,f_var,x1)); Sxk_new = double(subs(S,f_var,x1)); tol = norm(dx); % step6: 计算中止条件 norm(dx)<eps 是否满足,不满足转step 7 if tol<=eps break; end % step7: if Sxk_new < Sxk+beta*transpose(delta_Sxk)*dx u = u/v; break; else u = u*v; continue; end end x = x1;endt = x1minf = double(subs(S,f_var,t))
测试的结果是正确的。
参考第二本书中的例子把上述算法改成了一个多变量的程序,基本上没什么改动
% opti_LM_test2% 测试了 数值分析 Timothy Sauer 中 4.4节中的 4.19例clear all;clc;close all;x1 = -1; y1 = 0;x2 = 1; y2 = 1/2;x3 = 1; y3 = -1/2;R1 = 1; R2 = 1/2; R3 = 1/2;%syms x y;r1 = sqrt( (x-x1)^2 + (y-y1)^2 )-R1;r2 = sqrt( (x-x2)^2 + (y-y2)^2 )-R2;r3 = sqrt( (x-x3)^2 + (y-y3)^2 )-R3;r = ... [r1; r2; r3]%f = rclear r1 r2 r3 R1 R2 R3 x1 x2 x3 y1 y2 y3 x y r;%% S = transpose(f)*ff_var = symvar(f)t_init = [0 0] % 初始值,要给出u = 2v = 1.5beta = 0.4eps = 1.0e-6tol = 1%% x = t_initjacobian_f = jacobian(f,f_var)%% while tol>eps fxk = double(subs(f,f_var,x)); Sxk = double(subs(S,f_var,x)); % step2: 计算 fxk Sxk delta_fxk = double(subs(jacobian_f,f_var,x)); % step3: 计算 delta_fxk delta_Sxk = transpose(delta_fxk)*fxk; % step4: 计算 delta_Sxk while 1 % step5: 计算Q,并解方程(Q+uI)delta_x = -delta_Sxk Q = transpose(delta_fxk)*delta_fxk; dx = -(Q+u*eye(size(Q)))\delta_Sxk; x1 = x + dx‘; % 注意转置 fxk = double(subs(f,f_var,x1)); Sxk_new = double(subs(S,f_var,x1)); tol = norm(dx); % step6: 计算中止条件 norm(dx)<eps 是否满足,不满足转step 7 if tol<=eps break; end % step7: if Sxk_new < Sxk+beta*transpose(delta_Sxk)*dx u = u/v; break; else u = u*v; continue; end end x = x1;end%% format short;opti_var_value = http://www.mamicode.com/x1>
结果也是正确的
细节和原理以后再补充
Levenberg-Marquardt 的 MATLAB 代码
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。