Chapter3加强版：解线性方程组的迭代法

概述#

这一章的核心是：

对大规模、稀疏线性方程组 $Ax=b$ ，直接法往往代价太高，迭代法通过构造一个序列 $x^{(0)},x^{(1)},\dots$ ，让它逐步逼近真解 $x^*$ 。

迭代法要回答两个问题：

收不收敛：当 $k\to\infty$ 时，是否有 $x^{(k)}\to x^*$ ？
收敛多快：为了达到给定精度，要迭代多少步？

本章主要方法：

Jacobi 迭代法：每一步完全使用上一轮旧值
Gauss-Seidel 迭代法：一旦算出新分量，立刻使用新值
SOR 松弛法：在 Gauss-Seidel 的基础上引入松弛因子 $\omega$
最速下降法：把对称正定线性方程组转成二次函数极小值问题
共轭梯度法：对对称正定稀疏系统特别重要的迭代法

这一章的主线可以概括成：

1
Ax=b
2
  ↓ 改写成不动点问题
3
x=Mx+g
4
  ↓ 迭代
5
x^{(k+1)}=Mx^{(k)}+g
6
  ↓ 分析误差
7
x^{(k)}-x^* = M^k(x^{(0)}-x^*)
8
  ↓ 判断收敛
9
ρ(M)<1

课前过渡：条件数与迭代法动机#

条件数的意义#

上节课最后得到过扰动估计：若

A(x+\delta x)=b+\delta b,

则在合适范数下有

\frac{\|\delta x\|}{\|x\|} \leq \|A\|\,\|A^{-1}\|\frac{\|\delta b\|}{\|b\|}.

定义

\mathrm{cond}(A)=\|A\|\,\|A^{-1}\|.

它表示：右端项 $b$ 的相对扰动可能被放大多少倍后传递到解 $x$ 上。

直观理解：

$b$ 往往来自测量，例如外力、载荷、边界条件
$A$ 是模型离散后得到的系数矩阵
$x$ 是真正想求的未知量
如果 $\mathrm{cond}(A)$ 很大，测量误差会被严重放大

所以条件数大的问题，本身就很难算准。

[Slides 图占位：插入【lesson4】上课手写 PPT 演示第 2-3 页，内容为 $A(x+\delta x)=b+\delta b$ 与 $\mathrm{cond}(A)=\|A\|\|A^{-1}\|$ 的手写推导。]

WARNING
条件数大时，不要指望单纯换一个求解程序就能彻底解决问题。
可能原因有两类：

模型 / 离散 / 变量选择导致矩阵写得不好

原问题本身就是病态问题，任何合理建模都会导致接近奇异的矩阵

第一类可以尝试换模型、换变量、重新尺度化；第二类需要正则化、约束、先验信息或更换问题表述。

Hilbert 矩阵：病态矩阵的经典例子#

Hilbert 矩阵定义为

H_{ij}=\frac{1}{i+j-1},\qquad 1\le i,j\le n.

三阶 Hilbert 矩阵为

H_3=\begin{bmatrix} 1 & \frac12 & \frac13\\ \frac12 & \frac13 & \frac14\\ \frac13 & \frac14 & \frac15 \end{bmatrix}.

四阶 Hilbert 矩阵为

H_4=\begin{bmatrix} 1 & \frac12 & \frac13 & \frac14\\ \frac12 & \frac13 & \frac14 & \frac15\\ \frac13 & \frac14 & \frac15 & \frac16\\ \frac14 & \frac15 & \frac16 & \frac17 \end{bmatrix}.

MATLAB 中可以直接用：

1
hilb(3)
2
hilb(4)
3
cond(hilb(4))
4
cond(hilb(10))

课堂演示中的数量级：

$\mathrm{cond}(H_3)\approx 5.24\times 10^2$
$\mathrm{cond}(H_4)\approx 1.55\times 10^4$
$\mathrm{cond}(H_{10})\approx 1.60\times 10^{13}$

这说明条件数增长极快。双精度浮点数大约只有 $15\sim16$ 位有效数字，若误差被放大 $10^{13}$ 倍，最后可能只剩下很少几位可信数字。

[Slides 图占位：插入【lesson4】上课手写 PPT 演示第 7、11、12 页，内容为 Hilbert 矩阵定义、 $H_3$ 书写、MATLAB cond(hilb(4)) 与 cond(hilb(10)) 输出。]

非方阵方程组与线性最小二乘#

前面主要讨论 $A$ 是方阵时的 $Ax=b$ 。实际中也常见非方阵：

方程数多于未知数：超定方程组
方程数少于未知数：欠定方程组

课堂举了线性回归的例子。

假设有 12 个月销量数据，希望拟合直线

y_i=ax_i+b, \qquad i=1,2,\dots,12.

未知量只有两个：

\begin{bmatrix}a\\ b\end{bmatrix},

但方程有 12 个：

\begin{cases} y_1=a\cdot 1+b,\\ y_2=a\cdot 2+b,\\ \cdots\\ y_{12}=a\cdot 12+b. \end{cases}

写成矩阵形式：

A=\begin{bmatrix} 1&1\\ 2&1\\ \vdots&\vdots\\ 12&1 \end{bmatrix},\qquad x=\begin{bmatrix}a\\ b\end{bmatrix},\qquad b=\begin{bmatrix}y_1\\ y_2\\ \vdots\\ y_{12}\end{bmatrix}.

通常没有精确解，于是求最小二乘解：

\min_x\|Ax-b\|_2^2.

正规方程为

A^TAx=A^Tb.

若 $A^TA$ 可逆，则

x=(A^TA)^{-1}A^Tb.

这也可以看成广义逆的一种形式。

WARNING
正规方程是最小二乘的直接做法，思路简单，但稳定性不够好。
原因是 $A^TA$ 会放大条件数问题。实际计算中更常用 QR 分解或 SVD，这部分后续章节会继续讨论。

为什么要学迭代法#

第二章的直接法，例如 Gaussian 消去、LU 分解，理论上可以给出精确解。问题在于：

计算量大：一般规模下复杂度约为 $O(n^3)$
会破坏稀疏性：原来大量为零的元素，在消去过程中可能变成非零元素
大规模稀疏系统无法承受填充：例如流体计算、结构力学、偏微分方程离散后，矩阵规模可能极大

稀疏矩阵的核心价值在于“只存非零元”。如果直接法造成大量 fill-in，就会让存储量和计算量暴涨。

迭代法的思想是：

不对矩阵做大规模分解，尽量保留 $A$ 的稀疏结构，只通过矩阵-向量运算逐步逼近解。

课前补充：Gram-Schmidt 正交化#

Lesson5 开头先回顾了线性代数中的 Gram-Schmidt 正交化。它不是本章迭代法本身的算法，但它是理解共轭梯度法的预备知识。

核心问题是：已知一组线性无关向量

v_1,v_2,\dots,v_n,

希望把它们改造成一组两两正交的向量

u_1,u_2,\dots,u_n,

并且保持它们张成的空间不变。

课堂直观解释是：

新向量中“不听话”的部分，就是它在已有方向上的投影。把投影减掉，剩下的部分就和已有方向正交。

[Slides 图占位：插入【lesson5】手写 PPT Slide 10-28，展示二维、三维情况下“减去投影”的几何图。]

二维情形：减去一个投影#

先令

u_1=v_1.

要让第二个方向与 $u_1$ 正交，就从 $v_2$ 中减去它在 $u_1$ 方向上的投影。

$v_2$ 在 $u_1$ 上的投影为

\operatorname{proj}_{u_1}v_2 = \frac{(v_2,u_1)}{(u_1,u_1)}u_1.

于是

u_2 = v_2- \frac{(v_2,u_1)}{(u_1,u_1)}u_1.

验证正交性：

(u_2,u_1) = (v_2,u_1)- \frac{(v_2,u_1)}{(u_1,u_1)}(u_1,u_1) =0.

这就是课堂上反复强调的“把投影部分剪掉”。

三维与一般情形#

三维时，先得到 $u_1,u_2$ ，再处理 $v_3$ 。此时 $v_3$ 需要同时去掉在 $u_1,u_2$ 上的投影：

u_3 = v_3- \frac{(v_3,u_1)}{(u_1,u_1)}u_1 - \frac{(v_3,u_2)}{(u_2,u_2)}u_2.

一般地，若已经得到两两正交的

u_1,u_2,\dots,u_k,

则

u_{k+1} = v_{k+1} - \sum_{i=1}^{k} \frac{(v_{k+1},u_i)}{(u_i,u_i)}u_i.

如果还要得到标准正交基，只需单位化：

q_i=\frac{u_i}{\|u_i\|}.

于是

(q_i,q_j)=0\quad(i\ne j), \qquad \|q_i\|=1.

TIP
这个公式不建议死背。忘记时按几何图像重新推：

保留当前向量 $v_{k+1}$

找出它在每个旧方向 $u_i$ 上的投影

把这些投影全部减掉

为什么它和共轭梯度法有关#

Gram-Schmidt 做的是普通内积下的正交化：

(u_i,u_j)=u_i^Tu_j=0.

共轭梯度法要做的事情很像，只是把普通正交换成关于 $A$ 的正交：

(d_i,Ad_j)=d_i^TAd_j=0.

如果 $A$ 对称正定，那么

(x,y)_A=x^TAy

本身就是一种内积。于是 $A$ -共轭可以理解为：

在被 $A$ 拉伸 / 压扁后的几何空间里正交。

所以 lesson5 先回顾 Gram-Schmidt，是为了后面自然引出：

1
普通正交化：减掉普通投影
2
A-共轭化：减掉 A-内积意义下的投影

这正是共轭梯度法搜索方向构造的核心思想。

迭代法的一般形式#

从方程组到不动点迭代#

设线性方程组

Ax=b.

若能将它等价改写为

x=Mx+g,

则可以从任意初值 $x^{(0)}$ 出发，构造迭代序列：

x^{(k+1)}=Mx^{(k)}+g,\qquad k=0,1,2,\dots

这里：

$M$ 称为 迭代矩阵
$g$ 是常向量
$x^{(0)}$ 是初始猜测

若 $x^{(k)}\to x^*$ ，并且 $x^*=Mx^*+g$ ，那么 $x^*$ 就是原方程组的解。

TIP
初值 $x^{(0)}$ 理论上可以任取。
实际计算中：

若有物理经验，尽量取接近真实解的初值

全零初值很常见，但有时过于特殊

随机初值或全一初值也可以作为尝试

误差递推#

设真解为 $x^*$ ，满足

x^*=Mx^*+g.

迭代式减去不动点方程：

x^{(k+1)}-x^*=M(x^{(k)}-x^*).

递推得到

x^{(k)}-x^*=M^k(x^{(0)}-x^*).

所以迭代能否收敛，关键看

M^k\to 0.

如果存在某个矩阵范数使

\|M\|<1,

则

\|x^{(k)}-x^*\|\le \|M\|^k\|x^{(0)}-x^*\|\to 0.

这给出一个很直观的充分条件。

Jacobi 迭代法#

分量形式#

设

A=(a_{ij})_{n\times n},\qquad b=(b_1,b_2,\dots,b_n)^T.

第 $i$ 个方程为

a_{i1}x_1+\cdots+a_{ii}x_i+\cdots+a_{in}x_n=b_i.

若 $a_{ii}\ne 0$ ，把 $x_i$ 解出来：

x_i=\frac{1}{a_{ii}}\left(b_i-\sum_{j\ne i}a_{ij}x_j\right).

Jacobi 迭代的做法是：右端全部使用上一轮旧值。

x_i^{(k+1)}=\frac{1}{a_{ii}}\left(b_i-\sum_{j\ne i}a_{ij}x_j^{(k)}\right), \qquad i=1,2,\dots,n.

也可以写成

x_i^{(k+1)}=-\frac{1}{a_{ii}}\sum_{j\ne i}a_{ij}x_j^{(k)}+\frac{b_i}{a_{ii}}.

矩阵形式#

记

D=\mathrm{diag}(a_{11},a_{22},\dots,a_{nn}).

则 Jacobi 迭代为

x^{(k+1)}=(I-D^{-1}A)x^{(k)}+D^{-1}b.

也就是

M_J=I-D^{-1}A, \qquad g_J=D^{-1}b.

若把矩阵分裂写成

A=D-L-U,

其中 $L$ 是严格下三角部分取负， $U$ 是严格上三角部分取负，则

x^{(k+1)}=D^{-1}(L+U)x^{(k)}+D^{-1}b.

NOTE
Jacobi 的特点：

每个分量 $x_i^{(k+1)}$ 都只依赖上一轮 $x^{(k)}$

各分量之间可以并行计算

需要同时保存旧向量和新向量

收敛速度通常较慢

算法步骤#

输入： $A,b,x^{(0)}$ ，最大迭代次数 $N$ ，容忍精度 tol。

对 $k=0,1,2,\dots,N-1$
对 $i=1,2,\dots,n$ ，计算

x_i^{(k+1)}=\frac{1}{a_{ii}}\left(b_i-\sum_{j\ne i}a_{ij}x_j^{(k)}\right).

计算残差

r^{(k+1)}=b-Ax^{(k+1)}.

若 $\|r^{(k+1)}\|<\text{tol}$ ，停止。

WARNING
使用 Jacobi 前至少要保证 $a_{ii}\ne 0$ 。若对角元为零，可以尝试交换方程顺序。方程顺序改变后，迭代矩阵也会改变，收敛性可能随之改变。

例 1：Jacobi 迭代求三元方程组#

求解

\begin{cases} 10x_1-x_2-2x_3=72,\\ -x_1+10x_2-2x_3=83,\\ -x_1-x_2+5x_3=42. \end{cases}

解出各分量：

\begin{cases} x_1^{(k+1)}=0.1x_2^{(k)}+0.2x_3^{(k)}+7.2,\\ x_2^{(k+1)}=0.1x_1^{(k)}+0.2x_3^{(k)}+8.3,\\ x_3^{(k+1)}=0.2x_1^{(k)}+0.2x_2^{(k)}+8.4. \end{cases}

取

x^{(0)}=(0,0,0)^T.

第一步：

x^{(1)}=(7.2,8.3,8.4)^T.

第二步：

\begin{aligned} x_1^{(2)}&=0.1\times 8.3+0.2\times 8.4+7.2=9.71,\\ x_2^{(2)}&=0.1\times 7.2+0.2\times 8.4+8.3=10.70,\\ x_3^{(2)}&=0.2\times 7.2+0.2\times 8.3+8.4=11.50. \end{aligned}

继续迭代可得：

$k$	$x_1^{(k)}$	$x_2^{(k)}$	$x_3^{(k)}$
0	0.0000	0.0000	0.0000
1	7.2000	8.3000	8.4000
2	9.7100	10.7000	11.5000
3	10.5700	11.5700	12.4820
4	10.8535	11.8534	12.8282
5	10.9510	11.9510	12.9414
6	10.9834	11.9834	12.9804
7	10.9944	11.9981	12.9934
8	10.9981	11.9941	12.9978
9	10.9994	11.9994	12.9992

精确解为

x^*=(11,12,13)^T.

可以看到，Jacobi 迭代逐渐靠近精确解。

Gauss-Seidel 迭代法#

Jacobi 的问题在于：计算 $x_2^{(k+1)}$ 时，明明 $x_1^{(k+1)}$ 已经算出来了，却仍然使用旧值 $x_1^{(k)}$ 。

Gauss-Seidel 的思想是：

只要新值已经算出，后面的分量立刻使用新值。

分量形式#

Gauss-Seidel 迭代为

x_i^{(k+1)}=\frac{1}{a_{ii}}\left( b_i- \sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)}- \sum_{j=i+1}^{n}a_{ij}x_j^{(k)} \right).

这里：

$j<i$ 的分量已经更新，用 $x_j^{(k+1)}$
$j>i$ 的分量尚未更新，用 $x_j^{(k)}$

矩阵形式#

仍记

A=D-L-U.

则

(D-L)x^{(k+1)}=Ux^{(k)}+b.

因此

x^{(k+1)}=(D-L)^{-1}Ux^{(k)}+(D-L)^{-1}b.

即

M_{GS}=(D-L)^{-1}U, \qquad g_{GS}=(D-L)^{-1}b.

TIP
Gauss-Seidel 的两个优势：

通常比 Jacobi 快

可以原地更新，只需要保存一个向量 $x$

但它不保证总是比 Jacobi 收敛，也存在 Jacobi 收敛而 Gauss-Seidel 发散的例子。

例 2：Gauss-Seidel 迭代继续求例 1#

仍解

\begin{cases} 10x_1-x_2-2x_3=72,\\ -x_1+10x_2-2x_3=83,\\ -x_1-x_2+5x_3=42. \end{cases}

取 $x^{(0)}=(0,0,0)^T$ 。

第一步：

x_1^{(1)}=\frac{1}{10}(72)=7.2000.

计算 $x_2^{(1)}$ 时已经使用新值 $x_1^{(1)}$ ：

x_2^{(1)}=\frac{1}{10}(x_1^{(1)}+2x_3^{(0)}+83) =\frac{1}{10}(7.2000+83)=9.0200.

计算 $x_3^{(1)}$ 时继续使用新值 $x_1^{(1)},x_2^{(1)}$ ：

x_3^{(1)}=\frac{1}{5}(x_1^{(1)}+x_2^{(1)}+42) =\frac{1}{5}(7.2000+9.0200+42)=11.6440.

迭代结果：

$k$	$x_1^{(k)}$	$x_2^{(k)}$	$x_3^{(k)}$
0	0.0000	0.0000	0.0000
1	7.2000	9.0200	11.6440
2	10.4308	11.6719	12.8050
3	10.9313	11.9572	12.9778
4	10.9913	11.9947	12.9972
5	10.9989	11.9993	12.9996
6	10.9999	11.9999	13.0000

和 Jacobi 相比，Gauss-Seidel 在这个例子中明显更快。

松弛法与 SOR 迭代#

从 Gauss-Seidel 到松弛修正#

Gauss-Seidel 每一步会给出一个修正量

\Delta x=x_{GS}^{(k+1)}-x^{(k)}.

松弛法把修正量乘上一个参数 $\omega$ ：

x^{(k+1)}=x^{(k)}+\omega\Delta x.

也就是：

$\omega=1$ ：普通 Gauss-Seidel
$0<\omega<1$ ：低松弛，步子变小
$1<\omega<2$ ：超松弛，即 SOR

SOR 全称为 Successive Over-Relaxation。

SOR 的分量形式#

x_i^{(k+1)}=(1-\omega)x_i^{(k)}+ \frac{\omega}{a_{ii}} \left( b_i- \sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)}- \sum_{j=i+1}^{n}a_{ij}x_j^{(k)} \right).

这可以理解成：

1
新值 = 旧值 + ω × Gauss-Seidel 修正量

SOR 的矩阵形式#

仍设

A=D-L-U.

SOR 迭代可以写成

(D-\omega L)x^{(k+1)}=[(1-\omega)D+\omega U]x^{(k)}+\omega b.

所以迭代矩阵为

M_\omega=(D-\omega L)^{-1}[(1-\omega)D+\omega U].

右端项为

g_\omega=\omega(D-\omega L)^{-1}b.

例 3：SOR 迭代#

取

\omega=1.4, \qquad x^{(0)}=(1,1,1)^T,

求解

\begin{cases} 2x_1-x_2=1,\\ -x_1+2x_2-x_3=0,\\ -x_2+2x_3=1.8. \end{cases}

由 SOR 公式得

\begin{cases} x_1^{(k+1)}=-0.4x_1^{(k)}+0.7(1+x_2^{(k)}),\\ x_2^{(k+1)}=-0.4x_2^{(k)}+0.7(x_1^{(k+1)}+x_3^{(k)}),\\ x_3^{(k+1)}=-0.4x_3^{(k)}+0.7(1.8+x_2^{(k+1)}). \end{cases}

第一步：

\begin{aligned} x_1^{(1)}&=-0.4+0.7(1+1)=1,\\ x_2^{(1)}&=-0.4+0.7(1+1)=1,\\ x_3^{(1)}&=-0.4+0.7(1.8+1)=1.56. \end{aligned}

继续迭代：

$k$	$x_1^{(k)}$	$x_2^{(k)}$	$x_3^{(k)}$
0	1.0000	1.0000	1.0000
1	1.0000	1.0000	1.5600
2	1.0000	1.3920	1.6184
3	1.2744	1.4682	1.6404
4	1.2180	1.4136	1.5934
5	1.2023	1.3916	1.6068
6	1.1932	1.4034	1.6007
7	1.2051	1.4027	1.6016
8	1.1999	1.4000	1.5994
9	1.2000	1.3996	1.6001

精确解为

x^*=(1.2,1.4,1.6)^T.

第 9 步已经非常接近真解。

松弛因子如何选#

SOR 的效果高度依赖 $\omega$ 。

一般事实：

收敛必要条件： $0<\omega<2$
当 $A$ 为对称正定矩阵时，SOR 在 $0<\omega<2$ 下收敛
$\omega=1$ 退化为 Gauss-Seidel
$\omega$ 太接近 $2$ 可能导致振荡或发散

当 $A$ 是某些标准三对角矩阵时，可以估计最佳松弛因子：

\omega_{\mathrm{opt}} = \frac{2}{1+\sqrt{1-\rho(B)^2}},

其中 $B$ 是 Jacobi 迭代矩阵。

NOTE
实际计算中， $\omega$ 往往靠问题结构和试算选择。
老师课堂给出的经验是：大规模问题中 $\omega$ 常常会取得比较接近 $2$ ，例如 $1.9$ 以上，但必须小于 $2$ ，并且要结合具体矩阵结构试算。

迭代法的收敛条件#

谱半径#

设 $A$ 的特征值为

\lambda_1,\lambda_2,\dots,\lambda_n.

定义矩阵 $A$ 的 谱半径 为

\rho(A)=\max_{1\le i\le n}|\lambda_i|.

谱半径和矩阵范数有重要关系：

\rho(A)\le \|A\|.

并且对任意 $\varepsilon>0$ ，存在某种矩阵范数，使

\|A\|\le \rho(A)+\varepsilon.

因此，谱半径是判断迭代矩阵长期作用效果的核心量。

一般迭代法的充要条件#

对迭代

x^{(k+1)}=Mx^{(k)}+g,

若对任意初值 $x^{(0)}$ 都收敛到唯一不动点，则充要条件为

\rho(M)<1.

理由很直接：

x^{(k)}-x^*=M^k(x^{(0)}-x^*).

所以必须有

M^k\to 0.

而矩阵幂 $M^k\to 0$ 的充要条件就是

\rho(M)<1.

TIP
实用判断层次：

若能算出 $\rho(M)$ ，看 $\rho(M)<1$

若不方便算谱半径，可以找一个范数验证 $\|M\|<1$

$\|M\|<1$ 是充分条件， $\rho(M)<1$ 是充要条件

常用判别条件#

设 $Ax=b$ 。

严格对角占优#

若对每一行都有

|a_{ii}|>\sum_{j\ne i}|a_{ij}|,

则称 $A$ 严格对角占优。

若 $A$ 严格对角占优，则 Jacobi 与 Gauss-Seidel 迭代均收敛。

不可约弱对角占优#

若

|a_{ii}|\ge \sum_{j\ne i}|a_{ij}|

对所有行成立，并且至少有一行严格成立，同时矩阵不可约，则称为不可约弱对角占优。

在这种情况下，Jacobi 与 Gauss-Seidel 也有收敛性保证。

对称正定#

若 $A$ 是对称正定矩阵，则：

Gauss-Seidel 迭代收敛
SOR 迭代在 $0<\omega<2$ 时收敛

WARNING
对称正定不能直接保证 Jacobi 一定收敛。
判断 Jacobi 时仍要看对应迭代矩阵 $M_J$ 的谱半径，或者使用对角占优等充分条件。

方程顺序会影响收敛#

同一个方程组，如果交换方程顺序，解不变，但迭代矩阵会改变。

所以：

直接法通常不太关心方程顺序对收敛的影响
迭代法的收敛性可能因方程顺序改变而改变

实际写程序前，常常会通过重排、预处理、尺度化来改善迭代性质。

停机准则：不要只看前后两步差#

教材给出了误差估计：若 $\|M\|<1$ ，则

\|x^{(k)}-x^*\|\le \frac{\|M\|}{1-\|M\|}\|x^{(k)}-x^{(k-1)}\|.

因此在 $\|M\|$ 不太接近 1 时，可以用

\|x^{(k)}-x^{(k-1)}\|<\varepsilon

作为停止依据。

但老师课堂特别强调：实际编程时，不应只用前后两步差作为停机准则。

更可靠的做法是看残差：

r^{(k)}=b-Ax^{(k)}.

常用停机准则：

\|r^{(k)}\|<\text{tol},

或相对残差：

\frac{\|b-Ax^{(k)}\|}{\|b\|}<\text{tol}.

WARNING
前后两步差很小，只说明迭代序列局部变化很慢。
如果收敛非常慢， $x^{(k)}$ 和 $x^{(k-1)}$ 可能已经很接近，但它们离真解还很远。
实际程序建议同时设置：

残差阈值 tol

最大迭代次数 max_iter

必要时记录残差曲线

最速下降法#

对称正定方程组与二次函数#

设 $A$ 对称正定。

考虑二次函数

\varphi(x)=\frac12x^TAx-b^Tx.

它的梯度为

\nabla\varphi(x)=Ax-b.

令梯度为零：

\nabla\varphi(x)=0 \quad\Longleftrightarrow\quad Ax=b.

由于 $A$ 对称正定， $\varphi(x)$ 是严格凸二次函数，存在唯一全局极小值点。

所以求解 $Ax=b$ 等价于求

\min_x \varphi(x).

负梯度方向#

在点 $x^{(k)}$ ，下降最快方向是负梯度方向：

-\nabla\varphi(x^{(k)})=b-Ax^{(k)}.

记残差

r^{(k)}=b-Ax^{(k)}.

则最速下降法沿 $r^{(k)}$ 前进：

x^{(k+1)}=x^{(k)}+\alpha_k r^{(k)}.

这里 $\alpha_k$ 是步长。

[Slides 图占位：插入课堂手写示意图：对称正定二次函数像一个碗，从 $x^{(0)}$ 沿负梯度方向下降。]

最优步长#

最速下降法每一步选择当前方向上的最佳步长：

\alpha_k=\arg\min_{\alpha}\varphi(x^{(k)}+\alpha r^{(k)}).

代入并对 $\alpha$ 求导，可得

\alpha_k= \frac{(r^{(k)},r^{(k)})}{(Ar^{(k)},r^{(k)})}.

因此迭代公式为

x^{(k+1)}= x^{(k)}+ \frac{(r^{(k)},r^{(k)})}{(Ar^{(k)},r^{(k)})}r^{(k)}.

NOTE
这里的内积为
$(x,y)=x^Ty=\sum_{i=1}^n x_iy_i.$

从一维线搜索重新推导步长#

课堂中对最优步长又完整推了一遍。关键是把多元问题临时变成一元问题。

已知 $x^{(k)}$ ，记

r^{(k)}=b-Ax^{(k)}=-\nabla\varphi(x^{(k)}).

沿负梯度方向走：

x^{(k+1)}=x^{(k)}+\alpha r^{(k)}.

此时 $x^{(k)}$ 和 $r^{(k)}$ 都已知，唯一未知量是 $\alpha$ 。令

g(\alpha)=\varphi(x^{(k)}+\alpha r^{(k)}).

把

\varphi(x)=\frac12x^TAx-b^Tx

代入，得到

g(\alpha) = \frac12(x^{(k)}+\alpha r^{(k)})^TA(x^{(k)}+\alpha r^{(k)}) -b^T(x^{(k)}+\alpha r^{(k)}).

对 $\alpha$ 求导。利用 $A=A^T$ ，交叉项可以合并：

g'(\alpha) = (r^{(k)})^T(Ax^{(k)}-b) + \alpha (r^{(k)})^TA r^{(k)}.

又因为

Ax^{(k)}-b=-r^{(k)},

所以

g'(\alpha) = -(r^{(k)})^Tr^{(k)} + \alpha (r^{(k)})^TA r^{(k)}.

令 $g'(\alpha)=0$ ，得

\alpha_k = \frac{(r^{(k)})^Tr^{(k)}}{(r^{(k)})^TA r^{(k)}}.

这就是前面最速下降法中的步长公式。

[Slides 图占位：插入【lesson5】手写 PPT Slide 30-72，展示二次型、线搜索和 $\alpha_k$ 推导。]

NOTE
这里 $A$ 的对称性很关键。若 $A$ 不对称，二次型梯度会出现 $\frac12(A+A^T)x-b$ ，不能直接得到 $Ax=b$ 。

最速下降中相邻残差正交#

由更新公式

x^{(k+1)}=x^{(k)}+\alpha_k r^{(k)},

可得残差更新：

r^{(k+1)} =b-Ax^{(k+1)} =r^{(k)}-\alpha_kA r^{(k)}.

于是

(r^{(k+1)},r^{(k)}) =(r^{(k)},r^{(k)})- \alpha_k(Ar^{(k)},r^{(k)}).

代入

\alpha_k= \frac{(r^{(k)},r^{(k)})}{(Ar^{(k)},r^{(k)})},

得到

(r^{(k+1)},r^{(k)})=0.

因此，在精确线搜索下，最速下降法相邻两步的残差方向正交。

WARNING
这里说的是相邻残差正交。不要把它误读成所有历史残差都自动两两正交。共轭梯度法通过构造 $A$ -共轭方向，能获得更强的正交结构。

为什么最速下降可能很慢#

若 $A$ 的最大特征值和最小特征值分别为

\lambda_{\max},\lambda_{\min},

则最速下降的收敛速度受

\frac{\lambda_{\max}-\lambda_{\min}}{\lambda_{\max}+\lambda_{\min}}

控制。

若 $\lambda_{\max}\gg\lambda_{\min}$ ，这个比例接近 1，收敛会很慢。

几何图像：

如果等高线接近圆形，负梯度方向几乎直指极小值
如果等高线是很扁的椭圆，最速下降会在谷底两侧来回“之”字形摆动

[Slides 图占位：插入课堂手写示意图：椭圆等高线、负梯度方向与最速下降的锯齿形路径。]

这解释了为什么最速下降法虽然简单，但实际大规模计算中常常不够快。

补充写成条件数形式：

\kappa(A)=\frac{\lambda_{\max}}{\lambda_{\min}}.

则

\frac{\lambda_{\max}-\lambda_{\min}} {\lambda_{\max}+\lambda_{\min}} = \frac{\kappa(A)-1}{\kappa(A)+1}.

所以：

$\kappa(A)\approx 1$ 时，等高线接近圆，最速下降非常快
$\kappa(A)\gg 1$ 时，等高线很扁，下降路径容易来回折返

如果 $A$ 的所有特征值几乎相同，那么 $A$ 接近某个常数倍单位阵，二次函数的等高线接近圆形；若特征值差异很大，等高线被拉成长椭圆，最速下降会出现明显的“锯齿形”。

[Slides 图占位：插入【lesson5】手写 PPT Slide 88-98，展示圆形 / 扁椭圆等高线与最速下降锯齿形路径。]

共轭梯度法#

$A$ -共轭方向#

最速下降法每一步都沿当前负梯度方向走，容易在狭长椭圆等高线中来回摆动。

共轭梯度法的改进思想是：选择一组关于 $A$ 共轭的搜索方向。

若

d_i^TAd_j=0, \qquad i\ne j,

则称 $d_i,d_j$ 关于 $A$ 共轭，也叫 $A$ -正交。

当 $A$ 对称正定时，

(x,y)_A=x^TAy

可以看成一种新的内积。

所以 $A$ -共轭本质上是在被 $A$ 改变尺度后的空间里正交。

TIP
普通正交：
$d_i^Td_j=0.$
$A$ -共轭：
$d_i^TAd_j=0.$
可以把 $A$ 理解成把空间“压扁 / 拉伸”的尺度矩阵。

算法公式#

共轭梯度法适用于 $A$ 对称正定。

取初值 $x^{(0)}$ ，令

r^{(0)}=b-Ax^{(0)}, \qquad d^{(0)}=r^{(0)}.

第 $k$ 步：

\alpha_k=\frac{(r^{(k)},r^{(k)})}{(d^{(k)},Ad^{(k)})},

x^{(k+1)}=x^{(k)}+\alpha_k d^{(k)},

r^{(k+1)}=r^{(k)}-\alpha_kAd^{(k)},

\beta_k=\frac{(r^{(k+1)},r^{(k+1)})}{(r^{(k)},r^{(k)})},

d^{(k+1)}=r^{(k+1)}+\beta_kd^{(k)}.

这是常用实现形式。

课本也写成等价形式：

\beta_{k-1}=-\frac{(r^{(k)},Ad^{(k-1)})}{(d^{(k-1)},Ad^{(k-1)})},

d^{(k)}=r^{(k)}+\beta_{k-1}d^{(k-1)},

\lambda_k=\frac{(r^{(k)},d^{(k)})}{(d^{(k)},Ad^{(k)})},

x^{(k+1)}=x^{(k)}+\lambda_kd^{(k)}.

NOTE
共轭梯度法的重要性质：

精确算术下， $n$ 维对称正定系统最多 $n$ 步得到精确解

实际浮点计算中受舍入误差影响，通常作为迭代法使用

每步主要需要矩阵-向量乘 $Ad^{(k)}$ ，适合大规模稀疏矩阵

从 Gram-Schmidt 到 $A$ -共轭化#

最速下降法每次取

d^{(k)}=r^{(k)}.

共轭梯度法第一步仍然这样做：

d^{(0)}=r^{(0)}=b-Ax^{(0)}.

第一步就是沿负梯度方向做一次最优线搜索：

x^{(1)}=x^{(0)}+ \alpha_0d^{(0)}.

但从第二步开始，不再直接取 $d^{(1)}=r^{(1)}$ ，而是把 $r^{(1)}$ 改造成与旧方向 $d^{(0)}$ 关于 $A$ 共轭的方向。

设

d^{(1)}=r^{(1)}+\beta_0d^{(0)}.

要求

(d^{(1)},Ad^{(0)})=0.

代入：

(r^{(1)}+\beta_0d^{(0)},Ad^{(0)})=0.

于是

\beta_0 =-\frac{(r^{(1)},Ad^{(0)})}{(d^{(0)},Ad^{(0)})}.

这和 Gram-Schmidt 的形式完全对应：

1
新方向 = 当前残差 - 当前残差在旧方向上的“坏分量”

只是这里的“坏分量”按 $A$ -内积来判断。

[Slides 图占位：插入【lesson5】手写 PPT Slide 99-145，展示从 $d^{(0)}=r^{(0)}$ 到 $d^{(k)}=r^{(k)}+\beta_{k-1}d^{(k-1)}$ 的推导。]

一般地，若已有

d^{(0)},d^{(1)},\dots,d^{(k-1)}

两两 $A$ -共轭，则取

d^{(k)}=r^{(k)}+\beta_{k-1}d^{(k-1)},

其中

\beta_{k-1} =-\frac{(r^{(k)},Ad^{(k-1)})}{(d^{(k-1)},Ad^{(k-1)})}.

然后沿 $d^{(k)}$ 做线搜索：

x^{(k+1)}=x^{(k)}+ \alpha_kd^{(k)}.

对

\varphi(x^{(k)}+\alpha d^{(k)})

关于 $\alpha$ 求导，得

\alpha_k = \frac{(r^{(k)},d^{(k)})}{(d^{(k)},Ad^{(k)})}.

在常用实现中，因为 $r^{(k)}$ 与旧方向正交，且 $d^{(k)}=r^{(k)}+$ 旧方向组合，所以

(r^{(k)},d^{(k)})=(r^{(k)},r^{(k)}).

于是得到前面算法公式里的

\alpha_k = \frac{(r^{(k)},r^{(k)})}{(d^{(k)},Ad^{(k)})}.

为什么只需要减去前一个方向#

如果照搬 Gram-Schmidt，似乎应该写成

d^{(k)}=r^{(k)}- \sum_{i=0}^{k-1} \frac{(r^{(k)},Ad^{(i)})}{(d^{(i)},Ad^{(i)})}d^{(i)}.

但共轭梯度法实际只需要处理 $d^{(k-1)}$ 。原因在于：

前面每一步都做了精确线搜索
残差 $r^{(k)}$ 已经和旧的搜索子空间满足正交关系
在对称正定二次问题中，新残差的结构会自动消去更早方向上的分量

所以只剩下最近的方向 $d^{(k-1)}$ 需要修正。

这就是共轭梯度法比直接做完整 Gram-Schmidt 更省的地方：

1
完整 Gram-Schmidt：每次减掉所有旧方向的投影
2
共轭梯度法：每次只用前一个方向递推

NOTE
这条结论依赖 $A$ 对称正定和精确算术。实际浮点计算会破坏完全正交 / 完全共轭，所以工程计算中常用残差阈值判断停止，必要时会重启或使用预条件技术。

有限步收敛与条件数估计#

共轭梯度法有两个层次的收敛结论。

第一，精确算术下，若 $A$ 有 $r$ 个互异特征值，则共轭梯度法最多 $r$ 步得到精确解。因此也有粗略结论：

r\le n,

所以 $n$ 维问题最多 $n$ 步精确收敛。

这个结论解释了课堂中的直观例子：

若所有特征值都相同，等高线是圆 / 球，一步到达极小值
若只有两个互异特征值，例如某些旋转椭球结构，两步即可精确收敛

第二，实际计算更关心误差下降估计。令

\kappa(A)=\frac{\lambda_{\max}(A)}{\lambda_{\min}(A)}.

则常见估计为

\frac{\|x^{(k)}-x^*\|_A}{\|x^{(0)}-x^*\|_A} \le 2\left( \frac{\sqrt{\kappa(A)}-1}{\sqrt{\kappa(A)}+1} \right)^k.

这里

\|e\|_A=\sqrt{e^TAe}.

和最速下降法相比：

\text{最速下降因子：} \frac{\kappa(A)-1}{\kappa(A)+1},

\text{共轭梯度因子：} \frac{\sqrt{\kappa(A)}-1}{\sqrt{\kappa(A)}+1}.

共轭梯度法相当于把条件数开了根号，通常明显更快。但当 $\kappa(A)$ 极大时，这个因子仍然接近 1，收敛仍可能偏慢。

[Slides 图占位：插入【lesson5】手写 PPT Slide 150-156，展示 CG 收敛因子与条件数 $\kappa(A)$ 的关系。]

TIP
课堂中的判断逻辑是：

小规模矩阵：直接用 CG 往往已经够用

大规模稀疏矩阵：只靠原始 CG 可能还不够，需要改善条件数

改善条件数的核心工具就是预条件

课堂补充：预条件共轭梯度的基本想法#

lesson5 后半段在共轭梯度法后继续引出了 预条件（preconditioning） 思想。它仍然属于求解线性方程组迭代法的内容，与 Chapter4 特征值计算相接。

CG 的收敛速度受 $\kappa(A)$ 控制。理想情况下，希望把原方程

Ax=b

变成一个条件数更小的等价问题。

如果能找到一个容易计算的矩阵 $\widehat A^{-1}$ ，使

\widehat A^{-1}\approx A^{-1},

并且

\kappa(\widehat A^{-1}A)\ll \kappa(A),

那么就可以考虑求

\widehat A^{-1}Ax= \widehat A^{-1}b.

如果 $\widehat A^{-1}=A^{-1}$ ，那么左端矩阵就是单位阵，一步得到解。但这没有实际意义，因为求出 $A^{-1}$ 本身已经比解方程更贵。

所以预条件矩阵要同时满足两点：

足够接近 $A^{-1}$ 或能明显降低条件数
计算代价足够低，不能比原问题还难

WARNING
直接左乘 $\widehat A^{-1}$ 可能破坏对称正定结构。CG 需要对称正定矩阵，所以实际常用对称预条件。

一种常见想法是找 $C\approx A$ ，并令 $C$ 对称正定，且容易求解

Cz=r.

若 $C$ 有 Cholesky 型分解

C=LL^T,

则可把原方程等价变换为

L^{-1}AL^{-T}y=L^{-1}b, \qquad x=L^{-T}y.

如果

\kappa(L^{-1}AL^{-T})\ll \kappa(A),

CG 的收敛就会明显改善。

[Slides 图占位：插入【lesson5】手写 PPT Slide 156-259，展示 $\widehat A^{-1}A$ 降低条件数、Cholesky 型预处理与 MATLAB 演示。]

预条件共轭梯度法常用实现写成：

计算残差

r^{(k)}=b-Ax^{(k)}.

解预条件方程

Cz^{(k)}=r^{(k)}.

用 $z^{(k)}$ 代替普通 CG 中的残差方向来构造搜索方向。

常见递推为：

\alpha_k= \frac{(r^{(k)},z^{(k)})}{(d^{(k)},Ad^{(k)})},

x^{(k+1)}=x^{(k)}+ \alpha_kd^{(k)},

r^{(k+1)}=r^{(k)}- \alpha_kAd^{(k)},

解

Cz^{(k+1)}=r^{(k+1)},

再取

\beta_k= \frac{(r^{(k+1)},z^{(k+1)})}{(r^{(k)},z^{(k)})},

d^{(k+1)}=z^{(k+1)}+\beta_kd^{(k)}.

常见预条件选择包括：

Jacobi / 对角预条件：取 $C=D$ ，便宜但改善有限
Cholesky / 不完全 Cholesky 预条件：更接近 $A$ ，但需要控制计算和存储代价

NOTE
预条件的本质不是追求精确逆矩阵，而是在“便宜”和“有效降低条件数”之间折中。

例 4：共轭梯度法两步得到精确解#

求解

\begin{cases} 3x_1-x_2=2,\\ -x_1+x_2=0. \end{cases}

即

A=\begin{bmatrix}3&-1\\-1&1\end{bmatrix}, \qquad b=\begin{bmatrix}2\\0\end{bmatrix}.

取

x^{(0)}=(0,0)^T.

于是

r^{(0)}=b-Ax^{(0)}=(2,0)^T, \qquad d^{(0)}=(2,0)^T.

第一步：

\alpha_0=\frac{(r^{(0)},d^{(0)})}{(d^{(0)},Ad^{(0)})}=\frac{4}{12}=\frac13.

x^{(1)}=x^{(0)}+\alpha_0d^{(0)}= (0,0)^T+\frac13(2,0)^T= \left(\frac23,0\right)^T.

残差：

r^{(1)}=b-Ax^{(1)}=\left(0,\frac23\right)^T.

计算共轭修正：

\beta_0=-\frac{(r^{(1)},Ad^{(0)})}{(d^{(0)},Ad^{(0)})} =\frac19.

d^{(1)}=r^{(1)}+\beta_0d^{(0)} =\left(0,\frac23\right)^T+\frac19(2,0)^T =\left(\frac29,\frac23\right)^T.

第二步：

\alpha_1=\frac{(r^{(1)},d^{(1)})}{(d^{(1)},Ad^{(1)})} =\frac32.

x^{(2)}=x^{(1)}+\alpha_1d^{(1)} =\left(\frac23,0\right)^T+\frac32\left(\frac29,\frac23\right)^T =(1,1)^T.

因此两步得到精确解：

x^*=(1,1)^T.

应用实例：Poisson 方程与静电场#

课本最后给出一个静电场问题。

设平面上有一对正负电荷，它们周围形成电场。电势 $u(x,y)$ 满足 Poisson 方程：

\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2} =-\frac{\rho}{\varepsilon_p}.

其中：

$\rho(x,y)$ 是电荷密度
$\varepsilon_p$ 是介质电容常数
边界条件为 $u=0$

用二阶中心差分离散二阶导数：

\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^2} + \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h^2} =-\frac{\rho_{ij}}{\varepsilon_p}.

整理得五点差分格式：

u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}-4u_{i,j} =h^2\left(-\frac{\rho}{\varepsilon_p}\right)_{ij}.

把所有内部网格点未知量排成向量 $u$ ，就得到一个线性方程组：

Au=b.

这个矩阵 $A$ 具有典型结构：

规模为 $n^2\times n^2$
大部分元素为 0
块三对角结构
每个块内部也是三对角结构

因此它非常适合使用迭代法，尤其是 SOR。

课本例子中使用 SOR 计算，取 $n=30$ 时，最优松弛因子约为

\omega=1.811,

误差精度为

\varepsilon=10^{-5},

迭代约 70 步得到收敛结果。

[图占位：插入课本图 3-1， $n=30$ 时电位势数值解曲面。]

[图占位：插入课本图 3-2， $n=100$ 时电位势数值解曲面。]

这个例子说明：

偏微分方程离散后通常得到大规模稀疏线性系统。直接法容易破坏稀疏性，迭代法能够利用稀疏结构。

课堂中 lesson5 的内容分界可以这样理解：

Slide 5-259：围绕 Chapter3 的最速下降法、共轭梯度法、预条件思想展开
Slide 260 以后：进入矩阵特征值计算，是 Chapter4 开头内容，本文件不并入展开

本章知识框架#

1
线性方程组 Ax=b
2
│
3
├─ 直接法局限
4
│  ├─ O(n^3) 代价高
5
│  ├─ 稀疏矩阵会 fill-in
6
│  └─ 大规模 PDE 离散系统难以直接分解
7
│
8
├─ 一般迭代格式
9
│  ├─ x = Mx + g
10
│  ├─ x^{(k+1)} = Mx^{(k)} + g
11
│  ├─ e^{(k)} = M^k e^{(0)}
12
│  └─ 收敛充要条件：ρ(M)<1
13
│
14
├─ 基本迭代法
15
│  ├─ Jacobi：全用旧值，可并行，慢
16
│  ├─ Gauss-Seidel：用最新值，可原地更新
17
│  └─ SOR：引入松弛因子 ω，加速 GS
18
│
19
├─ 收敛判据
20
│  ├─ 谱半径 ρ(M)<1
21
│  ├─ 严格对角占优
22
│  ├─ 不可约弱对角占优
23
│  └─ 对称正定：GS 收敛，SOR 在 0<ω<2 收敛
24
│
25
├─ 停机准则
26
│  ├─ 理论上可估计 ||x^{(k)}-x^*||
27
│  ├─ 实际更推荐残差 ||b-Ax^{(k)}||
28
│  └─ 同时设置最大迭代次数
29
│
30
└─ 对称正定系统的优化观点
31
   ├─ Ax=b ⇔ min 1/2 x^T A x - b^T x
32
   ├─ 最速下降法：沿负梯度方向
33
   ├─ 共轭梯度法：沿 A-共轭方向，适合大规模稀疏 SPD 系统
34
   └─ 预条件思想：降低等价系统的条件数，加速 CG

这一章最重要的理解：

迭代法不是“精确消元”，而是构造一个不断改进的近似解序列。它牺牲有限步精确性，换取对大规模稀疏问题的可计算性。

概述#

目录#

课前过渡：条件数与迭代法动机#

条件数的意义#

Hilbert 矩阵：病态矩阵的经典例子#

非方阵方程组与线性最小二乘#

为什么要学迭代法#

课前补充：Gram-Schmidt 正交化#

二维情形：减去一个投影#

三维与一般情形#

为什么它和共轭梯度法有关#

迭代法的一般形式#

从方程组到不动点迭代#

误差递推#

Jacobi 迭代法#

分量形式#

矩阵形式#

算法步骤#

例 1：Jacobi 迭代求三元方程组#

Gauss-Seidel 迭代法#

分量形式#

矩阵形式#

例 2：Gauss-Seidel 迭代继续求例 1#

松弛法与 SOR 迭代#

从 Gauss-Seidel 到松弛修正#

SOR 的分量形式#

SOR 的矩阵形式#

例 3：SOR 迭代#

松弛因子如何选#

迭代法的收敛条件#

谱半径#

一般迭代法的充要条件#

常用判别条件#

严格对角占优#

不可约弱对角占优#

对称正定#

方程顺序会影响收敛#

停机准则：不要只看前后两步差#

最速下降法#

对称正定方程组与二次函数#

负梯度方向#

最优步长#

从一维线搜索重新推导步长#

最速下降中相邻残差正交#

为什么最速下降可能很慢#

共轭梯度法#

AAA-共轭方向#

算法公式#

从 Gram-Schmidt 到 AAA-共轭化#

为什么只需要减去前一个方向#

有限步收敛与条件数估计#

课堂补充：预条件共轭梯度的基本想法#

例 4：共轭梯度法两步得到精确解#

应用实例：Poisson 方程与静电场#

本章知识框架#

$A$ -共轭方向#

从 Gram-Schmidt 到 $A$ -共轭化#