写在前面:其他基础的不讲,网上到处都是,贴几个链接。
1.非线性最小二乘求解方法详解
2.非线性最小二乘法原理
直接上代码:
- 调用库函数。注意@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++ 代码待更新…