MATLAB | 求解线性方程组的最速下降法
转载请注明本文链接
及文章作者：slandarer
本期来讲求解线性方程组的最速下降法，本期主要是为了下一期内容做铺垫(下一篇大概会讲一下共轭梯度法)。
1 等价变分问题及其性质
假设 是对称正定矩阵，且 ，我们需要求解线性方程组：
此时我们会惊异的发现，线性方程组的解 与如下二次函数的极值点相同( )：
首先说明一下这个函数的性质，之后再说为啥求解线性方程组与求这个变分问题(变量分析问题)等价。
性质一
对于 ， 的偏导:
由于矩阵 的对称性( )：
因而有 的梯度:
性质二
一个关于变量线性组合的性质：
展开再约约同类项即可，在此不赘述。此处性质是为了与后面最速下降法结合(想象 为距离， 为方向)。
性质三
假设代入线性方程组的解 ，则有：
2 等价性证明
我们要证明线性方程组的解 就是使二次函数达到最小值的点，即：
对于一切 ：
而由于 为正定矩阵，则有：
当且仅当取值为 使二次函数达到最小值。由此我们已经能说明，求解线性方程组的过程可以等价为求二次函数极小值点的问题，而我们有非常多求极小值点的方法！
3 最速下降法
众所周知，函数值在梯度的反方向下降最快，因此我们可以每轮都使变量值往梯度反方向移动一定的距离，即从 出发，不断找到下降最快方向 :
我们将 叫做剩余向量，有时也用符号 来进行表示.
每次往下降最快的方向移动一段距离 :
使函数达到当前方向最小值：
那么这个 咋求呢，我们很自然的想到找偏导为0的点，由第1部分性质二得：
求偏导并令其等于0得：
可得：
即：
这样不断迭代的方法就是最速下降法。
4 收敛性
由于每轮都在找某方向上的使函数值最小的变量值，因此 为单调下降有界序列，由于最小值点 的存在性：
实时上存在如下性质：
其中 分别是对称正定矩阵最大最小特征值， ，
同时由于：
即两个相邻的搜索方向是相互正交的。
因此会产生锯齿现象，在开始几步，目标函数下降较快；但在接近极小点时，收敛速度就不理想了(明明走红色线相同距离更快，但是由于正交性要走黑色线)。
5 相关代码
这里给出一段通过最速下降法求解线性方程组的解的代码(注意要求 正定对称)：
A = [2 -2; -2 4];
b = [0; 2];
x = zeros(size(b));
epn = inf; % 初始精度
tol = 1e-6; % 要求精度
maxItr = 1e3; % 最大轮数
k = 1;
while k<=maxItr && epn>tol
p = b-A*x; % 求解梯度反方向
a = dot(p,p)./dot(A*p,p); % 求解\alpha值
x = x+a.*p; % 更新变量x值
k = k+1;
epn = norm(a.*p);
end
disp('最速下降法求解结果为:')
disp(x)
disp('线性方程组的解为:')
disp(A\b)
% 最速下降法求解结果为:
% 1.0000
% 1.0000
% 线性方程组的解为:
% 1
% 1
咱再给出一段带画图的代码：
A = [2 -1; -1 4];
b = [0; 2];
x = [-5; -5];
% 构造等价函数
vx = sym('x',size(b));
phi = vx.'*A*vx./2-b.'*vx;
phi = matlabFunction(phi);
% 坐标区域基础修饰
ax = gca; view(-71,38); grid on
ax.NextPlot = 'add';
ax.Projection = 'perspective';
ax.LineWidth = .8;
ax.XMinorTick = 'on';
ax.YMinorTick = 'on';
ax.ZMinorTick = 'on';
ax.TickLength = [.02,.01];
ax.GridLineStyle = ':';
ax.FontName = 'Cambria';
% 基本绘图
[X1,X2] = meshgrid(linspace(-5,5,50));
surf(X1,X2,phi(X1,X2),'EdgeColor','none');
contour(X1,X2,phi(X1,X2),20,'LineWidth',1)
% 创建gif图
DelayTime = .5;
F = getframe(ax);
[imind,cm] = rgb2ind(F.cdata,256);
imwrite(imind,cm,'test.gif','gif','Loopcount',inf,'DelayTime',DelayTime);
epn = inf; % 初始精度
tol = 1e-6; % 要求精度
maxItr = 1e3; % 最大轮数
k = 1;
while k<=maxItr && epn>tol
p = b-A*x; % 求解梯度反方向
a = dot(p,p)./dot(A*p,p); % 求解\alpha值
x = x+a.*p; % 更新变量x值
xx = [x-a.*p, x];
yy = [phi(xx(1,1),xx(2,1)); phi(xx(1,2),xx(2,2))];
% 绘制搜索过程
plot3(xx(1,:),xx(2,:),yy.','LineWidth',1,'LineStyle','--',...
'Color','r','Marker','o','MarkerFaceColor','r','MarkerSize',3)
plot3(xx(1,:),xx(2,:),[0,0],'LineWidth',1,'LineStyle','--',...
'Color','r','Marker','o','MarkerFaceColor','r','MarkerSize',3)
% 制作gif时至展示前几帧
if k<6
pause(.5);
F = getframe(ax);
[imind,cm] = rgb2ind(F.cdata,256);
imwrite(imind,cm,'test.gif','gif','WriteMode','append','DelayTime',DelayTime);
end
k = k+1;
epn = norm(a.*p);
end
参考
[1] 《数值分析》(第五版)李庆扬等编,清华大学出版社
以上已经是全部代码，当然不方便复 制的话我新建了一个数学建模函数gitee库，刚刚开始建立，可以去那里下载:
https://gitee.com/slandarer/math-model
