首页 > 代码库 > 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 代码