Levenberg-Marquardt 的 MATLAB 代码

Posted 2020-08-08 菠菜僵尸

参考资料：

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))

结果也是正确的

细节和原理以后再补充