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$ 的稀疏结构，只通过矩阵-向量运算逐步逼近解。

迭代法的一般形式#

从方程组到不动点迭代#

设线性方程组

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.$

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

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

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

则最速下降的收敛速度受

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

控制。

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

几何图像：

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

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

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

共轭梯度法#

$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)}$ ，适合大规模稀疏矩阵

例 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= rac{(r^{(1)},d^{(1)})}{(d^{(1)},Ad^{(1)})} =\frac32.

x^{(2)}=x^{(1)}+\alpha_1d^{(1)} =\left(\frac23,0\right)^T+ rac32\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$ 时电位势数值解曲面。]

这个例子说明：

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

本章知识框架#

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 系统

这一章最重要的理解：

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

概述#

目录#

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

条件数的意义#

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

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

为什么要学迭代法#

迭代法的一般形式#

从方程组到不动点迭代#

误差递推#

Jacobi 迭代法#

分量形式#

矩阵形式#

算法步骤#

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

Gauss-Seidel 迭代法#

分量形式#

矩阵形式#

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

松弛法与 SOR 迭代#

从 Gauss-Seidel 到松弛修正#

SOR 的分量形式#

SOR 的矩阵形式#

例 3：SOR 迭代#

松弛因子如何选#

迭代法的收敛条件#

谱半径#

一般迭代法的充要条件#

常用判别条件#

严格对角占优#

不可约弱对角占优#

对称正定#

方程顺序会影响收敛#

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

最速下降法#

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

负梯度方向#

最优步长#

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

共轭梯度法#

AAA-共轭方向#

算法公式#

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

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

本章知识框架#

$A$ -共轭方向#