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

結果也是正確的

細節和原理以后再補充