共轭梯度法matlab代码_线性方程组(4)-变分原理与共轭梯度法

上一节讲到了静态迭代法。因为收敛太慢，更常用的还是非静态迭代法。

这里先从变分原理入手。考虑函数

若

对称，则函数的梯度

若

对称正定，则解方程

等价于求函数

的最小值。方程的残差

即为函数

的负梯度。

---

一维优化问题

假设已经得到了

步迭代的解

，在第

步中，沿

方向搜索最优解，即

代入函数表达式可得

最优的

取值应是

得到新的解后，新的残差

---

最速下降法（Steepest descent method）

搜索方向的最简单的取法就是取

，即负梯度方向，这就是优化问题中很常用的最速下降法。用MATLAB代码描述

function [x, iter] = sd_solve(A, b, x0, epsilon, iter_max)
    iter = 0;
    x = x0;
    r = b - A * x0;
    rr = dot(r, r);
    while sqrt(rr) >= epsilon && iter < iter_max
        iter = iter + 1;
        Ar = A * r;
        alpha = rr / dot(r, Ar);
        x = x + alpha * r;
        r = r - alpha * Ar;
        rr = dot(r, r);
    end
end

注意提取出了

和

两个公共表达式，避免了重复计算。

每一步迭代需要计算一次矩阵向量乘，两次向量内积，两次向量线性组合。

注意到

的计算可能有误差积累，干扰迭代终止条件的判断，使得迭代在达到要求的精度之前终止。可能的解决方法有：用

来计算更精确的残差，引入一次额外的矩阵向量乘；或用更小的

；或用其它迭代终止条件。

假设精确解为

，则有结论

---

共轭梯度法（Conjugate gradient method, CG）

假设前

步使用的搜索方向是线性无关的向量

，而且前

步得到的

满足

且希望第

步得到的结果

满足

将

和

在子空间中展开

有

对于

由前面的条件

得到

这就是共轭梯度法，即选择两两

共轭的一系列搜索方向。仅有这个约束还不足以确定

，为了得到一个确定的搜索方向，我们选择将

对

正交化

其中

是使得

最小的系数（垂线最短）。

为了下面的推导，先提出一个引理，证明留作练习（提示：数学归纳法）

由

可知

正交于

，再根据引理可以推出下面的推导中所需要的所有正交关系。

因为

所以

因为

正交于

根据归纳假设，第二项在

时取最小值，所以

又因为我们只关心

的方向，不关心其大小，可以写成

的计算式可以继续简化

因为

所以

经过上面一大段推导，共轭梯度法终于呈现在我们眼前，用MATLAB代码描述

function [x, iter] = cg_solve(A, b, x0, epsilon, iter_max)
    iter = 0;
    x = x0;
    r = b - A * x0;
    p = r;
    rr_prev = dot(r, r);            % 上一步得到的解的残差平方和
    while sqrt(rr_prev) >= epsilon && iter < iter_max
        iter = iter + 1;
        Ap = A * p;
        alpha = rr_prev / dot(p, Ap);
        x = x + alpha * p;
        r = r - alpha * Ap;
        rr_curr = dot(r, r);        % 这一步得到的解的残差平方和
        beta = rr_curr / rr_prev;
        p = r + beta * p;           % 下一步的p
        rr_prev = rr_curr;
    end
end

每一步迭代需要计算一次矩阵向量乘，两次向量内积，三次向量线性组合。

注意，因为要求

所以初值

如果能找到一个更接近精确解的非零初值

，则等价于解

中的

，最后得到原方程的解

；由于

只涉及加法，根据加法交换律和结合律，这个累加放在开头还是结尾不影响结果，于是可以取

。

按照前面的理论，共轭梯度法迭代到第

步时

即得到精确解。也就是说，共轭梯度法是一种理论上的直接解法。

但是和LU分解法一样，在实际计算中，由于浮点误差的积累，理论上的精确解没有意义。我们只需要把共轭梯度法看成迭代法即可，达到迭代终止条件之前的迭代步数小于

或者大于

都没关系。这个现象在之后还会再遇到。

共轭梯度法的收敛性

对于条件数较大的矩阵，共轭梯度法的收敛速度比最速下降法快很多。

共轭梯度法要求系数矩阵正定，因为只有正定才存在最小值；实际上对于不定矩阵，共轭梯度法也很有可能收敛到梯度为零所对应的鞍点，不过不能从理论上保证收敛性。

---

共轭梯度法的变种

将共轭梯度法推广到非对称矩阵，最简单的办法是把方程

转换为

由于

对称正定，可以使用共轭梯度法解该方程。这就是正规方程的共轭梯度法。由于

所以收敛较慢。

双共轭梯度法（Bi-conjugate gradient method）也可以应用于非对称矩阵。计算两组

和

其中

维护

的双正交性和

的双共轭性，对于

双共轭梯度法有数值不稳定的现象，在迭代一定步数后，由于浮点误差积累，双正交性和双共轭性不再保持。稳定双共轭梯度法（Bi-conjugate gradient stablilized method）通过将一维搜索扩增到高维来修正这一误差，当搜索维数高时收敛速度提升，但单步迭代计算量和消耗的存储空间也将提升。