非线性数据拟合（最小二乘法Gauss-Newton&Levenberg-Marquardt(LM)）matlab and C++

写在前面：其他基础的不讲，网上到处都是，贴几个链接。
1.非线性最小二乘求解方法详解
2.非线性最小二乘法原理
直接上代码：

1. 调用库函数。注意@fun处应为误差函数，不是原函数。

clc;
clear all;
%**********************   *****************************%
%          原函数：y=ax^2+exp(bx)+cx+d                 %
%**********************   *****************************%
t = [1.2, 1, 3, -8];  % 待拟合参数
% [lsqnonlin—matlab帮助文件，很详细](https://ww2.mathworks.cn/help/optim/ug/lsqnonlin.html#buul_xw-1)
options = optimset('Algorithm', 'trust-region-reflective', 'Tolfun', 1e-14, 'MaxFunEvals', 1000);
[tt, resnom, residual, exitflg, output] = lsqnonlin(@fun, t, [], [], options);

% 待拟合点
t_0 = [1, 2, 5, -10];  % 准确原参
num = 200;
rng(5); noise = 0.01*(randi(10, 1, num)-5);
x = 0.01*(1: num);
y = t_0(1)*x.^2 + exp(t_0(2)*x) + t_0(3)*x + t_0(4) + noise;

% 拟合后点 画图
xx = x;
yy = tt(1)*x.^2 + exp(tt(2)*x) + tt(3)*x + tt(4);
figure(2);
plot(xx, yy);
hold on
plot(x, y, '.');
hold on;
plot(x, y-noise, "red")

function ret = fun(t)
t_0 = [1, 2, 5, -10];  % 准确原参
num = 100;
x = 0.01*(1: num);
rng(5);
noise = 0.2*(randi(10, 1, num)-5);
y = t_0(1)*x.^2 + exp(t_0(2)*x) + t_0(3)*x + t_0(4) + noise;
% 以上为求误差函数需要
ret = t(1)*x.^2 + exp(t(2)*x) + t(3)*x + t(4) - y;  % 误差函数
end

拟合后参数 tt = [3.42619, 1.62460, 4.79166, -9.93472] VS 准确原参 t_0 = [1, 2, 5, -10];
注：红色为目标曲线；蓝色为拟合后曲线。

2. Gauss-Newton法

clc;
clear all;

ar = 1; br = 2; cr = 5; dr = -10; % real coef
ae = 1.2; be = 1; ce = 3; de = -8; % org coe and eval coef

num = 100;
iterations = 1000;
rng(5); noise = 0.2*(randi(10, 1, num)-5);

x_data = 0.01*(1: num);
y_data = ar*x_data.^2 + exp(br*x_data) + cr*x_data + dr + noise; % 待拟合数据
y = ar*x_data.^2 + exp(br*x_data) + cr*x_data + dr;

J = zeros(4,1);
cost = 0; lastcost = 0;
for iter = 1:iterations
    H = 0; b = 0;
    for i = 1:num
        xi = x_data(i);
        yi = y_data(i);
        err = yi - (ae*xi.^2 + exp(be*xi) + ce*xi + de);
        % 雅克比矩阵 求偏导
        J(1) = -xi * xi;
        J(2) = -xi * exp(be * xi);
        J(3) = -xi;
        J(4) = -1;
        H = H + J * J';
        b = b - err .* J ;
        cost = err*err;
    end
    
    % 求解线性方程 Hx=b
    dx = H\b;
    if iter > 1 &&  cost >= lastcost
        break
    end
    ae = ae + dx(1);
    be = be + dx(2);
    ce = ce + dx(3);
    de = de + dx(4);
    lastcost = cost; 
end

yy = ae*x_data.^2 + exp(be*x_data) + ce*x_data + de;
abcd = [ae, be, ce, de]; 
% 作图
figure(1);
plot(y_data, ".");
hold on;
plot(yy, "blue");
hold on;
plot(y, "red");

拟合后参数 abcd = [3.42361, 1.62509, 4.79206, -9.93481] VS 准确原参 t_0 = [1, 2, 5, -10];

3. Levenberg-Marquardt(LM)
帮助理解该算法,参考下面论文：张鸿燕，耿 征，“Levenberg-Marquardt 算法的一种新解释”。

ar = 1; br = 2; cr = 5; dr = -10; % real coef
num = 200;
x_data = 0.01*(1:num);
yr = ar*x_data.^2 + exp(br*x_data) + cr*x_data + dr;
rng(5);  noise = 0.1*(randi(10, 1, num)-5);
y_data = yr + noise;

% 初始猜测
ae = 2; be = -2; ce = -5; de = 10;
% 数据个数
Ndata = length(y_data);
% 参数维数
Nparams = 4;
% 迭代最大次数
n_iters = 300;
% 最小步长
dp_min = 1e-7;
% LM算法的阻尼系数初值
lamda=0.01;
% step1: 变量赋值
updateJ=1;
% step2: 迭代
for it=1:n_iters
    if updateJ==1
        % 根据当前估计值，计算雅克比矩阵
        J=zeros(Ndata,Nparams);
        for i=1:length(x_data)
            J(i,:)=[x_data(i)*x_data(i) x_data(i)*exp(be*x_data(i)) x_data(i) 1];
        end
        % 根据当前参数，得到函数值
        ye = ae*x_data.^2 + exp(be*x_data) + ce*x_data + de;
        % 计算误差
        dd = y_data - ye;
        % 计算（拟）海塞矩阵
        H = J'*J;
        % 若是第一次迭代，计算误差
        if it==1
            err=dot(dd,dd); % 内积 按列相乘求和
        end
    end
    % 根据阻尼系数lamda混合得到H矩阵
    % 当updataJ=0时，表示随着迭代残差e不能减小，则放大lamda再计算残差值。
    H_lm = H + (lamda*eye(Nparams,Nparams));
    % 计算步长dp，并迭代参数估计值
    dp = inv(H_lm)*(J'*dd(:));
    g = J'*dd(:);
    a_lm = ae+dp(1);
    b_lm = be+dp(2);
    c_lm = ce+dp(3);
    d_lm = de+dp(4); 
    % 计算参数迭代后对应的y和计算残差e
    ye_lm = a_lm*x_data.^2 + exp(b_lm*x_data) + c_lm*x_data + d_lm;
    dd_lm = y_data - ye_lm;
    err_lm = dot(dd_lm,dd_lm);
    % 根据新的误差err_lm，决定该次迭代是否合适，如不合适调制lamda重新迭代。
    % 当误差进一步减小时，缩小lamda,即在增量方程中相当于增大H权值，使其趋近Gauss-Newton法。
    % 当误差无法再减小时，放大lamda,即在增量方程中相当于增大lamda*I权值，使其趋近梯度下降法。
    if err_lm<err
        ae = a_lm;
        be = b_lm;
        ce = c_lm;
        de = d_lm; 
        err = err_lm;
        if  min(abs(dp)) < dp_min
            disp("小于设定步长，结束迭代！");
            break;
        end
        lamda = lamda/10;
        updateJ=1;
    else
        % 说明步距过大，需放大阻尼系数lamda，缩小步距，并重新判断误差是否能减小（e_lm<e）
        updateJ=0;  
        lamda=lamda*10;
    end
    if it == n_iters
        disp("超出迭代次数！");
    end
end

abcd = [ae be ce de];

figure(4);
plot(x_data, y_data, ".");
hold on;
plot(x_data, yr, "red");
hold on;
plot(x_data, ae*x_data.^2 + exp(be*x_data) + ce*x_data + de, "blue");

目前这个最好使。
拟合后参数 abcd = [0.95660, 1.9995, 5.13014, -10.00388] VS 准确原参 t_0 = [1, 2, 5, -10];

三种方式可以调整拟合点数，拟合点值（noise，该值很大程度决定拟合效果），迭代次数，最小步长对比效果。

c++ 代码待更新…