标签：
参考资料：
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 = 2
v = 1.5
beta = 0.4
eps = 1.0e-6
x = 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;
end
t = x1
minf = 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 = r
clear r1 r2 r3 R1 R2 R3 x1 x2 x3 y1 y2 y3 x y r;
%%
S = transpose(f)*f
f_var = symvar(f)
t_init = [0 0] % 初始值，要给出
u = 2
v = 1.5
beta = 0.4
eps = 1.0e-6
tol = 1
%%
x = t_init
jacobian_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 = x1
minf = double(subs(S,f_var,opti_var_value))
结果也是正确的
细节和原理以后再补充
Levenberg-Marquardt 的 MATLAB 代码
标签：
原文地址：http://www.cnblogs.com/shepherd2015/p/5854836.html
