Levenberg-marquardt's MATLAB code

Resources:

1, "Proficient in MATLAB optimization calculation (2nd edition)" Shang, such as 9th Chapter 9.3 subsection L-m method

2, "Numerical analysis" Timothy Sauer of the 4th Chapter 4.4 Nonlinear Least Squares example

Although there is code in the first book, there are errors that fix the error

% opti_lm_test1% tests The example of l-m in the MATLAB optimization book, the result is the correct 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 The initial value of the argument%%u = 2v = 1.5beta = 0.4eps = 1.0e-6x = T_init;x = Transpose (x) ;% Jacobian_f = Jacobian (f,f_var); tol = 1;%% Subs It is not a value, but a symbol! Also convert to double type!!!    While tol>eps fxk = double (Subs (f,f_var,x));  SXK = Double (Subs (s,f_var,x));  % Step2: Calculates fxk sxk delta_fxk = double (Subs (jacobian_f,f_var,x));   % Step3: Calculates delta_fxk delta_sxk = Transpose (DELTA_FXK) *fxk; % STEP4: Calculate delta_sxk while 1 step5: Calculate q, and solution equation (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: Calculates the abort condition norm (DX) <eps is satisfied, does not meet the turn 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))

The result of the test is correct.

Refer to the example in the second book to change the above algorithm into a multivariable program, basically no change

% opti_lm_test2% tested the numerical analysis of 4.19 cases in section 4.4 of Timothy Sauer 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]% initial value, to give Out u = 2v = 1.5beta = 0.4eps = 1.0e-6tol = 1%% x = T_initjacobian_f = Jacobian (F,f_var) percent while tol>eps fxk = double (s    UBS (F,F_VAR,X));  SXK = Double (Subs (s,f_var,x));  % Step2: Calculates fxk sxk delta_fxk = double (Subs (jacobian_f,f_var,x));   % Step3: Calculates delta_fxk delta_sxk = Transpose (DELTA_FXK) *fxk; % STEP4: Calculate delta_sxk while 1 step5: Calculate q, and solution equation (q+ui) delta_x =-delta_sxk Q = Transpose (DELTA_FXK) *delta_        FXK;                DX =-(Q+u*eye (Size (Q))) \delta_sxk; X1 = x + dx ';       % Note transpose fxk = double (Subs (f,f_var,x1)); Sxk_new = Double (Subs (s,f_var,x1));   Tol = norm (dx);        % STEP6: Calculates the abort condition norm (DX) <eps is satisfied, does not meet the turn 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 = X1minf = Double (Subs (S,f_var,opti_var_value))

And the results are correct.

Details and principles added later

Levenberg-marquardt's MATLAB code