引入

对于系数矩阵对称正定的实线性方程组$Ax=b$,其解存在唯一，记为$x*.$由于$A$对称正定，令$|x|_A=\sqrt{x^\mathrm{T}Ax}$ ,则求$x*$等价于求如下问题

的最优解.由于$Ax_*=b$,则

记 $\varphi(x)=x^\mathrm{T}Ax-2b^\mathrm{T}x$
故求$x_*$的问题等价于求如下最优化问题

的最优解。

考虑利用迭代法求解最优化问题(4.2).首先任意给定初始向量$x_0$,沿下降方向$d_0$ ,经过步长$\alpha_0$得到一个新的向量$x_1=x_0+\alpha_0d_0$,使得

然后从$x_1$出发，再沿下降方向$d_1$前进步长$\alpha_1$,使得

此过程一直进行下去，直到求得问题(4.2)的最优解，不同的下降方向和步长确定策略可以给出来不同的迭代法.

最速下降法

原理&算法

首先考虑步长的选择策略.明显地，若下降方向$d_k$给定，则$\alpha_k$的选取应使得函数

在$x=x_k+\alpha_kd_k$处取得最小值.记$r_k=b-Ax_k$,对上式关于$\alpha$求导得

由此确定的$\alpha$即为所求的步长

其次考虑下降方向的选择策略。由微积分知识可知，若在局部要求下降最快，则下
降方向应为负梯度方向，故我们可以选择在$x_k$点的负梯度方向

作为下降方向，此时，

由$A$的正定性知，只要$rk\neq0$,有$r_k^\mathrm{T}d_k=r_k^\mathrm{T}r_k\neq0$,则$\varphi(x{k+1})<\varphi(x_k).$
此算法简单有效并且算法过程中主要是矩阵和向量的乘积及向量的内积运算，故
可以充分利用问题的稀疏性.综合上面的讨论可以得到以下算法流程：

算法 4.1 对称正定线性方程组的最速下降法

输入:矩阵$A$、向量$b$、初始向量$x_0$ 及精度$\varepsilon$

输出:向量$x$

计算残量$r_0=b-Ax_0;$

$k=0$

while $|r_k|\geqslant\varepsilon$ do

​ $k=k+1;$

​ 计算步长$\alphak-1=r{k-1}^{\mathrm{T}}r{k-1}/\left(r{k-1}^{\mathrm{T}}Ar_{k-1}\right);$

​ 计算更新后的向量$xk=x{k-1}+\alpha{k-1}r{k-1};$

​ 计算当前残量$r_k=b-Ax_k;$

end while

收敛性

引理1

定理1.18

若 $A$ 为Hermit矩阵(共轭转置等于本身),则

引理2

定理4.1

设$A$为对称正定矩阵，$p(x)$为实系数多项式，则对于任意向量$x$,有

证明 :由$A$为对称正定矩阵，则有谱分解$A=Q\Lambda Q^\mathrm{T}$,其中$Q$为正交矩阵，$\Lambda$为对角元均为正的对角矩阵.定义$A^{\frac12}=Q\Lambda ^{\frac12}Q^{\mathrm{T}}$,故对任意的$x$有：$|x|_A=\sqrt{x^\mathrm{T}Ax}=$ $\sqrt{x^T\boldsymbol{A}^{\frac12}\boldsymbol{A}^{\frac12}\boldsymbol{x}}=|\boldsymbol{A}^{\frac12}\boldsymbol{x}|_2$(因为2-范数的定义就是各个元素平方然后求和开方,等价于和自己做点乘然后开方).

我们再考虑 $p(\boldsymbol{A})$ ,设$p(\boldsymbol{A})=p_0+p_1x+\cdots+p_nx^n$ 因此 $p(\boldsymbol{A})=p_0+p_1A+\cdots+p_nA^n$ ,我们将谱分解$A=Q\Lambda Q^\mathrm{T}$带入其中,由于 $Q$ 的正定性,可以得到

所以:

最速下降法的收敛性定理:

最速下降法的收敛性定理

设系数矩阵 $A$ 的特征值为 $0<\lambda1\leq \lambda_2\cdots \lambda_n$ ,$ x*$ 为方程 $Ax=b$ 的解,则最速下降法产生的迭代序列 ${xk}{k=0}^{\infty}$ 满足:

收敛速度

从任一向量出发,最速下降法生成的迭代向量都会收敛到问题的解,并且其收敛速度依赖于 $\frac{\lambda_n-\lambda_1}{\lambda_1+\lambda_n}$的大小,当 $\lambda_n \approx \lambda_1$ 的时候,收敛速度很快,但是当 $\lambda_n >>\lambda_1$ 的时候,收敛速度很慢.

共轭梯度法

原理

共轭梯度法是在 20 世纪 50 年代初期由 Hestenes 和 Stiefel 首先提出的，目前已成为求解大型稀疏线性方程组的一类重要算法.从最速下降法的讨论可知，每步选取最速下降方向并不是最合适的，应选择更加合理有效的下降方向，并且方向的选择不应带来太多的额外计算量.
首先考虑步长的选择策略.假设已得到了迭代向量$x_k$及下降方向$d_k$,则下一步在
直线$x_k+\alpha\boldsymbol{d}_k$上寻找极小点，与最速下降法的讨论类似，可得

其中

下面考虑如何选择迭代方向$dk.$在初始步，可用信息少，此时选择负梯度方向$d_0=r_0=b-Ax_0$作为下降方向，在后续迭代步中，从式(4.3)可以看出，在下步步长的计算中$r_k$ 是必需的；而上步方向$d{k-1}$ 是已知的，故在选取当前步的下降方向 $dk$的时候可以充分利用 $r_k$ 以及 $d{k-1}$ 的信息,这并未引入额外的计算量.

由于

​

最后一个等号用到了(4.3)

这说明: 上一步的下降方向和下一步的残差垂直.

所以我们可以在过 $xk$ 由$r_k$ 及$d{k-1}$ 张成的平面

中寻找下一步迭代点,类似前述分析, $\eta$ 和 $\xi$必须满足:

计算并且整理可以得到

根据 $\boldsymbol{d}{k-1}^{\mathrm{T}} \boldsymbol{r}{k}=0$,可以得到第二个式子为 $0$.

若$ \boldsymbol{r}{k}\neq 0$,我们观察 $\xi$ ,如果 $\xi = 0$ ,那么为了使得第二个式子 $\xi \boldsymbol{r}{k}^{\mathrm{T}} \boldsymbol{A} \boldsymbol{d}{k-1}+\eta \boldsymbol{d}{k-1}^{\mathrm{T}} \boldsymbol{A} \boldsymbol{d}{k-1} =0$ 成立,由于 $A$ 是一个正定矩阵,$ \boldsymbol{r}{k}\neq 0$,说明此时迭代没有终止,因此 $d{k-1}\neq 0$ .所以 那么$\boldsymbol{d}{k-1}^{\mathrm{T}} \boldsymbol{A} \boldsymbol{d}_{k-1}>0$ ,此时立即得到 $\eta=0$ ,只有这样才可以使得这个式子成立.这个时候观察第一个式子,很显然,第一个式子不成立.

因此,$\xi\neq 0$

上述仅仅是为了保证 $\xi$作为分母的时候不为 $0$ .

我们可以取下降方向如下:

其中,

由此得到了下面的递推公式

实际应用中，上述公式可进行进一步简化.对于每步残量 $r{k+1}$,由于在计算 $\alpha_k$时$\boldsymbol{Ad}_k$ 已计算，故$\boldsymbol{r}{k+1}$ 可写为如下形式

由关系式(4.4)有$d{k-1}^\mathrm{T}r_k=0.$进而由$d_k=r_k+\beta{k-1}d_{k-1}$有

由$r{k+1}=r_k-\alpha_kAd_k$,并结合后面定理4.3中的关系式$r{k+1}^\mathrm{T}r_k=0$,有

其中,最后一步等式用到了步长$\alpha$ 的另一种表示

算法

算法 4.2 求解对称正定线性方程组的共轭梯度法

输入: 矩阵$\boldsymbol A$、向量$\boldsymbol b$、初始向量$\boldsymbol x_0$

输出:向量$x$
计算残量及方向$r_0=b-Ax_0,p_0=r_0$ ;
$k= 1$ ;
计算步长$\alpha_0=r_0^{\mathrm{T}}r_0/(p_0^{\mathrm{T}}Ap_0)$ ;

计算更新后向量$x_1=x_0+\alpha_0p_0;$

计算当前残量$r1=b-Ax_1$ ;
$\mathbf{while} $ $r_k\neq 0$ $\mathbf{do}$
$k= k+ 1$ ;
计算参数$\beta{k-2}=\dfrac{r{k-1}^\mathrm{T}r{k-1}}{(r{k-2}^\mathrm{T}r{k-2})};$
计算方向$pk-1=r{k-1}+\beta{k-2}p{k-2};$
计算步长$\alpha{k-1}=\dfrac{r{k-1}^\mathrm{T}r{k-1}}{(p{k-1}^\mathrm{T}Ap_{k-1})};$

​ 计算更新后的向量$xk=x{k-1}+\alpha{k-1}p{k-1};$

​ 计算当前残量$rk=r{k-1}-\alpha{k-1}Ap{k-1};$
$\mathbf{endwhile}$

共轭梯度法具有良好的性质，可以归纳为如下定理：

性质

定理4.3 共轭梯度法的性质

对于任意给定的初始向量$x0$,由共轭梯度法迭代$m$步产生的残量序列${\boldsymbol{r}_k}{k=0}^m$及方向序列${\boldsymbol{d}_k}^m_k=0$满足：

  • 当前步的残量与前面任一步的下降方向均正交，即

​

  • 任意两步的残量均正交，即

  • 任意两步的下降方向均关于矩阵$A$正交，即

  • 对于给定的非奇异矩阵$A\in\mathbb{R}^n\times n$及向量$r_0\in\mathbb{R}^n$,定义

    为$\mathbb{R}^n$关于$A$及$r_0$的 Krylov 子空间，其维数为

    其中$\mathrm{grad}(\boldsymbol{A},\boldsymbol{r}0)$为使得$p(\boldsymbol{A})\boldsymbol{r}_0=0$ 成立的所有首项系数为$1$W的多项式的最低次数.由残量的正交性可知，共轭梯度法最多经过$n$ 步 一 定 会 终 止 , 故 共 轭 梯 度 法 可 以 看作是直接法.但是由于实际使用中舍入误差的存在使得正交性很快会丢失，其有限步终止的性质也会不成立，故实际中一般作为迭代法使用.另外，在算法4.2中，部分运算是重复性的，例如 $r{k-1}^\mathrm{T}r{k-1},Ap{k-1}$ 等，这些重复计算可以通过存储已计算结果避免. 综合上述讨论，可得如下实用的算法：

    实用算法

    算法 4.3 求解对称正定线性方程组的共轭梯度法
    输入：矩阵$\boldsymbol A$、向量$\boldsymbol b$、初始向量$x$、精度$\varepsilon$ 及最大迭代步数$k_\mathrm{max}$
    输出：向量$x$
    $r=b-Ax,\gamma=r^{\mathrm{T}}r,p=r;$

    $r=r-\alpha\omega,\widetilde{\gamma}=\gamma,\gamma=r^{\mathrm{T}}r;$
    $\mathbf{while}\sqrt{\gamma}\geqslant\varepsilon\left|b\right|\&k<k_{\max}\mathbf{do}$

    $$k=k+1,\beta=\gamma/\widetilde{\gamma},p=r+\beta p;$$
$$\omega=Ap,\alpha=\gamma/(p^\mathrm{T}\omega),x=x+\alpha p;$$
$r=r-\alpha\omega,\widetilde{\gamma}=\gamma,\gamma=r^{T}r;$

    $\textbf{endwhile}$

定理4.4

定理4.4

令$x_k$为利用共轭梯度法得到的第$k$个迭代向量，则有

即:

在空间逐步增大的过程中. $\varphi(x_{k})$始终最优

定理4.5

关于共轭梯度法的收敛性估计，有以下定理成立

定理4.5

由共轭梯度法产生的迭代序列${xk}{k=0}^\infty$满足：

由定理4.5可知，共轭梯度法的收敛性依赖于系数矩阵的条件数.当系数矩阵的条件数接近于1或者大部分的特征值集中在一点附近时，共轭梯度法具有较快的收敛速度.