MATLAB | 求解线性方程组的共轭梯度法

上一期讲了求解线性方程组的最速下降法，那么自然而然的这期就是求解线性方程组的共轭梯度法，本文用到很多上一篇的知识，请自行回顾：MATLAB | 求解线性方程组的最速下降法

同时为了方便书写开(    的书写方式太占空间)本文一律使用    代替    .

共轭梯度法，顾名思义，就是在梯度下降法的定义中又引入了共轭的概念，这里的共轭指的是    -共轭。

1    -共轭

   对称正定，     为    维向量, 则可定义    -内积：

而若    则称    向量    -共轭，即    -内积正交。

若能找到一组互相    -共轭的向量作为一组基：

可将要求解的向量    用其线性表示：

则有

可将系数    求出进而求得    .

2 由线性无关向量获取    -共轭向量

其实原理还是施密特正交化，不过在操作过程中不使用内积而是使用    -内积，假设    是一组线性无关的向量，则令：

经此方法生成的    -共轭基向量依旧线性无关：

3 共轭梯度法(CG)

给定初值为0向量    ，计算残差:

取     计算:

(此处的    的含义在最速下降法中讲过，就还是找到    方向使函数最小的距离)

由    得到新的残量：    ，对    做    -投影得到    ：

重复进行下面操作：

这里为啥    只减去1项而不是减去k+1项就能    -共轭呢？？？

4 正交性证明

在这略微证明一下为啥这样得到的向量有如下性质：

实际上前面的递推公式可以改写为：

同时可以推出一个后面会用到的玩意：

由数学归纳法：

现假设    相互正交，    相互    -共轭,则有：

当    时，由    的表达式，    显然成立。当    时，由于归纳假设    

由：

得到：

因而有：

因此    与    均正交。再看    ，先看其与    的共轭性：

再看    时：

证毕。

5 收敛速度

由于残量    相互正交，当    时由于向量维数，在    中必有0向量，即残量为0，即必定能在    时求出精确解。

6 代码实现

共轭梯度法代码

% 求解线性方程组的共轭梯度法
% @author : slandarer
A = [ 3 -1  0  1 -1;
     -1  6  0  0  2;
      0  0  2  1  0;
      1  0  1  4  1;
     -1  2  0  1  5];
b = [9; 18; 1; 2; 6];
x = zeros(size(b));

epn    = inf;       % 初始精度
tol    = 1e-6;      % 要求精度
maxItr = size(A,1); % 最大轮数

k=1;
tic;
r = b-A*x;  % 计算r_0
P = r;      % 计算P_0 
while k<=maxItr && epn>tol
    alp = dot(r,P)/dot(A*P,P);    % 求解\alpha值
    x   = x+alp.*P;               % 更新变量x值
    r   = b-A*x;                  % 计算残量
    bet = dot(A*P,r)./dot(A*P,P); % 计算\beta值
    P   = r-bet.*P;               % 计算基向量
    
    k = k+1;
    epn = norm(r);
end
toc;

disp(['共进行[',num2str(k-1),']次迭代'])
disp('共轭梯度法求解结果为:')
disp(x)
disp('线性方程组的解为:')
disp(A\b)
% 历时 0.000023 秒。
% 共进行[5]次迭代
% 共轭梯度法求解结果为:
%     4.9367
%     3.4709
%     1.1418
%    -1.2835
%     1.0557
% 
% 线性方程组的解为:
%     4.9367
%     3.4709
%     1.1418
%    -1.2835
%     1.0557

与最速下降法对比

% 求解线性方程组的最速下降法
% @author : slandarer
A = [ 3 -1  0  1 -1;
     -1  6  0  0  2;
      0  0  2  1  0;
      1  0  1  4  1;
     -1  2  0  1  5];
b = [9; 18; 1; 2; 6];
x = zeros(size(b));

epn    = inf;     % 初始精度
tol    = 1e-6;    % 要求精度
maxItr = 1e3;     % 最大轮数

k=1;
tic;
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
toc;

disp(['共进行[',num2str(k-1),']次迭代'])
disp('最速下降法求解结果为:')
disp(x)
disp('线性方程组的解为:')
disp(A\b)
% 历时 0.000099 秒。
% 共进行[46]次迭代
% 最速下降法求解结果为:
%     4.9367
%     3.4709
%     1.1418
%    -1.2835
%     1.0557
% 线性方程组的解为:
%     4.9367
%     3.4709
%     1.1418
%    -1.2835
%     1.0557

可以看到同样的线性方程组，最速下降法需要46轮，而梯度下降法只用了5轮，轮次远远小于最速下降法，同时也小于矩阵阶数(虽然n轮肯定是精确解，但是不需要n轮就能得到误差范围内解)。

参考

[1] 《数值分析》(第五版)李庆扬等编,清华大学出版社
[2] 数值分析6(3共轭梯度法)苏州大学(https://www.bilibili.com/video/BV16a  4 y1t76z/)

以上已经是全部代码，当然不方便复 制的话我新建了一个数学建模函数gitee库，刚刚开始建立，可以去那里下载:

https://gitee.com/slandarer/math-model