定义
共轭梯度法是求解线性方程组 $Ax=b$ 的一种迭代算法,其中系数矩阵 $A$ 为 $n$ 阶实对称正定矩阵。算法从初始猜测解出发,沿着相互关于 $A$ 共轭的搜索方向逐步修正近似解,直到残差足够小[2]。它对矩阵的要求是它必须对称正定,这是该方法能保证收敛的前提[4]。
两个非零向量 $u$ 与 $v$ 关于 $A$ 共轭,指的是二者满足 $u^{T}Av=0$。当 $A$ 对称正定时,这一条件等价于在一个由 $A$ 定义的内积下二者正交,共轭关系本身是对称的[2]。在计算数学的分类中,共轭梯度法属于 克雷洛夫子空间 迭代法,是求解此类系统最早的一种方法[1]。
在不考虑舍入误差的理想情况下,共轭方向法至多经过 $n$ 次迭代即可得到精确解,其中 $n$ 为方程组未知数的个数,这一性质称为有限步终止或二次终止性[5]。
原理
从代数上看,当 $A$ 对称正定时,解 $x^{*}$ 同时也是二次泛函 $\varphi(x)=\frac{1}{2}x^{T}Ax-b^{T}x$ 的唯一极小值点,求线性方程组因此可以转化为求该泛函的极小值[2]。泛函的梯度为 $Ax-b$,沿负梯度方向下降的思想正是最速下降法,但最速下降法前后两次的搜索方向相互正交,容易产生“之”字形路径,收敛缓慢[6][5]。
共轭梯度法改用一组关于 $A$ 共轭的方向作为搜索方向。第 $k$ 步的残差记为 $r_{k}=b-Ax_{k}$,它是当前点的负梯度;新的搜索方向由当前残差与上一个搜索方向组合而成,即 $p_{k}=r_{k-1}+\beta_{k-1}p_{k-1}$[6]。由于只保留前一个方向,算法只需短递推,存储量小、每步计算量固定[2]。
每步的步长由线搜索确定,使泛函在该方向上取极小:$\alpha_{k}=\frac{r_{k-1}^{T}r_{k-1}}{p_{k}^{T}Ap_{k}}$,随后更新解与残差 $x_{k}=x_{k-1}+\alpha_{k}p_{k}$、$r_{k}=r_{k-1}-\alpha_{k}Ap_{k}$,系数 $\beta_{k}$ 常取 $\frac{r_{k}^{T}r_{k}}{r_{k-1}^{T}r_{k-1}}$[6]。实际程序通常以残差范数与右端项范数之比小于给定容差作为停止条件。
该算法的收敛速度取决于系数矩阵的条件数 $\kappa$,误差的 $A$ 范数满足 $\frac{\|x_{k}-x^{*}\|_{A}}{\|x_{0}-x^{*}\|_{A}}\leq 2\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^{k}$,条件数越接近 1 收敛越快,迭代次数大致与 $\sqrt{\kappa}$ 成正比[5][7]。
发展历程
20 世纪 50 年代初,电子计算机的出现推动了数值算法的研究。美国国家标准局(NBS)数值分析研究所的 Magnus Hestenes 与来自苏黎世联邦理工学院的访问学者 Eduard Stiefel 各自独立得到同一方法,并在 1952 年合作发表了求解线性系统的共轭梯度法论文,刊于 NBS 研究杂志[3][8]。
在最初的论文中,该方法是被当作直接法而非迭代法表述的,有限步终止性质被看得很重。20 世纪 60 年代,由于朴素实现处理更大规模问题时效果不佳,方法一度少人问津[1][8]。
1964 年,Fletcher 与 Reeves 把共轭梯度思想推广到非线性函数的极小化,给出了求解非线性最优化问题的共轭梯度算式,此后的 Polak–Ribiére–Polyak、Dai–Yuan 等方法都属于这一系列[2][9]。据文献记载,把共轭梯度法明确当作迭代法来看待,始于 1971 年 Reid 的工作,他指出该方法最重要的应用领域是稀疏正定系统[8]。
20 世纪 70 年代初,随着预处理技术的引入和稀疏矩阵应用需求的增长,共轭梯度法重新受到重视,此后成为求解大型稀疏对称系统的主流方法[1]。1999 年,《Computing in Science and Engineering》将包括共轭梯度法在内的克雷洛夫子空间迭代评为 20 世纪十大算法之一;2002 年,苏黎世联邦理工学院与 NIST 联合举办了纪念该方法发表 50 周年的研讨会[3]。
应用
共轭梯度法最主要的用途是求解偏微分方程离散后产生的大型稀疏对称正定线性系统,这类系统在有限元、有限差分分析中十分常见,未知数可达数百万个[10][11]。结构力学分析、计算流体力学、地下水污染建模、核聚变计算和化学中的 Ornstein–Zernike 方程求解都使用这类求解器[3]。
在有限元分析中,共轭梯度法的一个突出优点是全局刚度矩阵无需组装,只需存储各单元的局部矩阵,存储量正比于单元数目,比直接法低一个数量级以上[12]。此外,算法主要由矩阵向量乘、向量内积和向量加乘三类运算构成,天然适合并行和图形处理器加速[10][11]。
实际使用中最常见的形式是预处理共轭梯度法(Preconditioned Conjugate Gradient,PCG),即先引入预处理矩阵 $M$ 改善系数矩阵的条件数,再执行迭代。常用的预处理子包括 Jacobi 预处理、不完全 Cholesky 分解和对称逐次超松弛等[13][11]。在图像恢复、电磁散射计算、结构工程和信号处理等领域也能看到该方法的直接应用[14][15]。
在算法库实现方面,1993 至 1994 年出版的迭代法模板集给出了共轭梯度法及其变体的标准程序框架,田纳西大学与橡树岭国家实验室提供的迭代模板例程以及 NIST 的 IML++ 库中都包含该算法的实现。
局限
共轭梯度法的收敛性保证以系数矩阵对称正定为前提。若矩阵对称但不定,二次泛函不再是严格凸的,某些曲率量可能取非正值,标准步长公式失去定义,算法可能失效;对于非对称矩阵,该方法也不保证收敛,需要改用双共轭梯度法、BiCGSTAB 或 GMRES 等方法[4][16]。
理论上的 $n$ 步终止在实际浮点计算中通常无法实现。舍入误差会破坏残差之间的正交性,使有限步终止失去保证;由于大型问题的 $n$ 往往极大,$O(n)$ 步迭代本身也超出可接受的计算量,因此实用算法大多采用最大迭代次数配合残差范数的双重停止准则。
收敛速度对条件数高度敏感,病态问题可能导致迭代极其缓慢甚至停滞,预处理正是为缓解这一困难而引入的[17][11]。预处理子本身也存在代价与风险:简单的 Jacobi 预处理在均匀网格的二维泊松方程上只相当于把矩阵乘以常数,几乎不改变条件数;不完全 Cholesky 分解在强各向异性问题上可能因舍弃填充项而失去正定性,导致分解失败[17]。此外,选择最优预处理子被证明是 NP 难问题,实际中往往依赖经验或试算[11]。
停止准则的选取本身也没有完全客观的依据。若只用相对残差下降若干数量级作为判据,它与近似解的真实误差之间缺乏明确联系,使用者面临继续迭代浪费算力或过早停止损失可靠性的两难[18]。残差范数在迭代过程中并不单调下降,可能出现震荡甚至上升,这也增加了判断收敛的难度[19]。
参见
参考资料
- NBS-Developed Method Named `Algorithm of the Century' [引用日期2026-10-14]
- 图文详情——科普中国资源服务 [引用日期2026-10-14]
- 50th Anniversary of Seminal NBS Publication Celebrated [引用日期2026-10-14]
- CS475-Lecture-Notes(PDF) [引用日期2026-10-14]
- slides_04C_CG(PDF) [引用日期2026-10-14]
- Conjugate Gradient Method - Chapters and Articles [引用日期2026-10-14]
- \[ [引用日期2026-10-14]
- A short history of Krylov methods [引用日期2026-10-14]
- S0362546X05002725 [引用日期2026-10-14]
- Conjugate Gradients - Chapters and Articles [引用日期2026-10-14]
- 基于深度学习的共轭梯度预条件自适应选择算法 [引用日期2026-10-14]
- pdf [引用日期2026-10-14]
- PCG [引用日期2026-10-14]
- 共辘梯度法 [引用日期2026-10-14]
- 预条件共轭梯度法(PCG)-知网阅读 [引用日期2026-10-14]
- 大型稀疏系统迭代法 | Bohrium [引用日期2026-10-14]
- 共轭梯度法 | Bohrium [引用日期2026-10-14]
- siam.org 上的网页 [引用日期2026-10-14]
- Electronic Transactions on Numerical Analysis [引用日期2026-10-14]