工程应用上,对于线性的方程组求解,将矩阵直接求逆的计算非常繁琐和复杂,可以使用梯度下降法来迭代求出方程组的数值解。首先任意假设一个初始值xini=[0,0]',然后根据局部线性误差的平方和最小的准则(代价函数J=min(e^2))来加入一个修正因子 。反复迭代直至相邻两次的迭代中修正因子足够小,满足误差精度时结束迭代过程。迭代的方向选取为代价函数的负梯度(逐渐减小的方向)。
在高数中,求解某一个函数的最小值,通常取函数的一阶导数为0的位置。对于线性方程组:Hx=C,H为m*n矩阵,x0为n*1矢量,C为m*1矢量。其理论上的解为:
下面探讨递推下降法求数值解,x0表示迭代过程中所得的数值解。
-------------------------------------------------------假装我是分割线-------------------------------------------------------------------
应用举例:数据拟合,具体线性方程组数值解
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%线性拟合曲线y=a+bx
%最小二乘LSE和梯度下降两种方法求解待估计得参数theta=[a,b]'
%
%最小二乘需要通过矩阵运算,梯度下降法则是数值迭代,代价函数取均方误差
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;clear ;close
y=[3 4 5 5 2 4 7 8 11 8 12 11 13 13 16 17 18 17 19 21]';%测量值
H=[ones(1,length(y));1:length(y)]'; %观测矩阵
theta_LSE = inv(H'*H)*H'*y;%参数的最小二乘解
y_LSE = H*theta_LSE;
%梯度下降法迭代,获取参数解
loop_times=1e3;%最大迭代次数
theta_gradient=zeros(loop_times,2)';
mu=0.01;%收敛步骤
theta_gradient(:,1)=[1,1]';%待估计参数初始值
for ii=1:loop_times
cost_function_derivation = H'*(H*theta_gradient(:,ii)-y)/length(y);
theta_gradient(:,ii+1)=theta_gradient(:,ii)-mu*cost_function_derivation;%沿着梯度的反方向迭代到误差曲面的最低点
if(norm(cost_function_derivation)<mu/100)%要求是代价函数最小时的参数,因此,代价函数的偏导数,在待估计参数上一拍处的取值应该接近于0,由于偏导数也是向量,取向量的模
break;
end
end
y_gradient = H*theta_gradient(:,end);%取迭代的最后一次参数的值,作为估计值
figure
plot(H(:,2),y,'b*','linewidth',2);grid on
hold on
plot(H(:,2),y_LSE,'k-','linewidth',0.5);
plot(H(:,2),y_gradient,'ro','linewidth',1);
hold off;
legend('original','LSE','gradient','location','best');
title('最小二乘解与梯度法数值解的对比');
-----------------------------------------------------------------------------------------------------------------------------------
%方程组用递推下降法的数值解法,动画展示迭代过程
% y=2*x+1
% y=3*x+5 解为:y=-7,x=-5
%以下是数值解,迭代计算,改写成方程组:Ax=b;
%其中向量x=[x,y],b=[-1,-5],A=[2 -1;3 -1];
clc;clear ;close all
x=-6:0.01:1;
y1=2*x+1;
y2=3*x+5;
b=[-1,-5]';
A=[2 -1;3 -1];
result=inv(A)*b;
x_init=[0,0]';%迭代中方程解的初始值,第一个元素是x,第二个元素是y
mu=0.08;%迭代步进,是否最终收敛与步进的大小有关
x_k(:,1)=x_init;
loop_times = 1e4;
count=0;%前后两次迭代值之差,连续多次小于某一门限的次数
figure
for ii=1:loop_times%迭代次数
err=A*x_k(:,ii)-b;
x_k(:,ii+1)=x_k(:,ii)-mu*A'*err;
plot(x,y1,'r-','linewidth',2);hold on;grid on
plot(x,y2,'b-','linewidth',2);
plot(0,0,'ko','MarkerFace','k','MarkerSize',10);
plot(result(1),result(2),'k^','MarkerFace','r','MarkerSize',10);
plot(x_k(1,1:ii),x_k(2,1:ii),'k','linewidth',1);
plot(x_k(1,end),x_k(2,end),'ks','MarkerFace','k','MarkerSize',10);
hold off;
legend('line-y1','line-y2','start-position','final-position','track','location','best');
drawnow;
recursive_err = x_k(:,ii+1)-x_k(:,ii);
if(norm(recursive_err)<mu/100)%迭代精度为迭代步长的1%
if(count>16)%连续16次迭代结果都近似相等,认为已经收敛
break;
else
count=count+1;
end
else
count=0;
end
end
clc
fprintf('梯度下降法解方程组的数值解(est_x,est_y)为:\n');
fprintf('%2.4f\n',x_k(:,end));
fprintf('方程组的真实解(x,y)为:\n');
fprintf('%2.4f\n',result);
线性方程组的数值解:
写在最后的话:
梯度下降法用处很广。在应用过程中需要注意矩阵形式的一阶导数是否正确求解;步长与迭代的速度成反比,与最终收敛精度成正比,二者需要权衡。
参考链接: