共轭梯度法：原理、设置与验证

共轭梯度法是求解对称正定稀疏系统的默认算法，每步只需一次矩阵—向量积、两次内积和 $O(n)$ 存储。它的收敛速率不取决于迭代上限，而取决于矩阵条件数 $\kappa$ 的平方根以及谱的聚集程度——这正是预条件能把迭代数从几千降到几十的原因。本文推导 A-范数误差界，用 Lanczos 过程解释为什么谱聚集比谱跨度更重要，并给出可手算的迭代数估算。预条件 CG 的配置只有五个可调项：预条件类型、填充级别或丢弃阈值、相对与绝对容差、范数口径、最大迭代数。以 $256\times256$ 结构化网格上的二维 Poisson 为基准（$N=66049$，$h=1/257$，$\kappa\approx4/(\pi^2h^2)=2.68\times10^4$），按上界估算 $\varepsilon=10^{-8}$ 需约 1563 步；CG 的日志只给一行"迭代 180 次，相对残差 3.1e-9"，这既不能说明解可信，也不能说明预条件已经到位。要判断收敛质量，需要把残差历史换算成实测收敛因子、与 $\kappa$ 预测值对比，再用一个离散误差为零的算例把代数误差单独隔离出来。以下数据取自 $64\times64$ 网格二维 Poisson（$N=3969$，$\kappa\approx1.71\times10^3$）与一维模型问题。

共轭梯度法在最小化什么

对对称正定 $A$，求解 $A\mathbf x=\mathbf b$ 等价于最小化严格凸二次泛函

$$\phi(\mathbf x)=\tfrac12\mathbf x^{\mathsf T}A\mathbf x-\mathbf b^{\mathsf T}\mathbf x, \qquad \nabla\phi=A\mathbf x-\mathbf b=-\mathbf r .$$

最速下降沿 $-\mathbf r$ 走，但在病态矩阵上会反复折返；CG 改为构造 $A$-共轭方向族 $\mathbf p_i^{\mathsf T}A\mathbf p_j=0\ (i\neq j)$，并在每一步做精确线搜索：

$$\alpha_k=\frac{\mathbf r_k^{\mathsf T}\mathbf r_k}{\mathbf p_k^{\mathsf T}A\mathbf p_k}, \qquad \mathbf x_{k+1}=\mathbf x_k+\alpha_k\mathbf p_k .$$

共轭性保证第 $k$ 步误差在 $\mathrm{span}\{\mathbf p_0,\dots,\mathbf p_k\}$ 上最优，因此无需保存历史向量，短递推就够。这也是 CG 与 GMRES 的根本分野：CG 用三项递推换内存，代价是必须保持对称正定。

有限精度下的正交性丢失

理论上 CG 至多 $n$ 步终止。以二维 Poisson 在 $64\times64$ 网格上为例，$N=3969$，$\kappa\approx4/(\pi^2h^2)=1.71\times10^3$（$h=1/65$），代入上界得 $\varepsilon=10^{-8}$ 时约需 396 步；而浮点运算中实测约 180 步即达标，说明上界偏保守约 2.2 倍。但正交性会被舍入逐步破坏，表现为递推残差 $\|\mathbf r_k\|$ 与真实残差 $\|\mathbf b-A\mathbf x_k\|$ 逐渐分离。长时间迭代后必须周期性重算真实残差；一旦两者相差超过两个数量级，就应停止并重新考虑预条件，而不是继续加迭代数。

用 Lanczos 三对角矩阵估计条件数

CG 的副产物就是 Lanczos 分解，把每步的 $\alpha_k$、$\beta_k$ 组装成三对角矩阵

$$T_k=\begin{pmatrix}\alpha_1&\beta_1&&\\ \beta_1&\alpha_2&\ddots&\\ &\ddots&\ddots&\beta_{k-1}\\ &&\beta_{k-1}&\alpha_k\end{pmatrix},$$

它的极值特征值（Ritz 值）单调逼近 $A$ 的极值。实测在 $k=20$ 步后 Ritz 最大值已稳定到 $7.98$，最小值降到 $4.7\times10^{-3}$，比值 $1.70\times10^3$，与理论 $\kappa$ 吻合到 1% 以内。这给出一条廉价的诊断路径：不必显式调用特征值求解器，从 CG 自身的递推系数就能得到 $\kappa$ 估计，进而判断预条件是否还有改进空间。若 Ritz 最小值在 20 步后仍持续下降，说明迭代尚未"看到"真正的低频模态，此时谈收敛因子为时过早。

容差、范数与停机判据

PCG 的收敛判据一般写成

$$\|\mathbf r_k\|_2 \le \max\!\left(\varepsilon_{\text{abs}},\,\varepsilon_{\text{rel}}\|\mathbf r_0\|_2\right),$$

常用 $\varepsilon_{\text{rel}}=1.0\times10^{-8}$、$\varepsilon_{\text{abs}}=1.0\times10^{-12}$。绝对容差不可省略：当右端接近零（例如增量形式的载荷步）时 $\|\mathbf r_0\|$ 本身极小，纯相对判据会让迭代在第 1 步就"收敛"。另一个易错点是范数口径——PETSc 默认用预条件后的范数 $\|\mathbf z_k\|_{M^{-1}}$，它与真实残差差一个 $M^{-1}$ 的尺度，在 IC 预条件下两者可能相差一个数量级。需要可对照的物理残差时应显式设置 -ksp_norm_type unpreconditioned，代价是每步多一次内积。

适用边界：SPD 之外一概不成立

CG 的推导依赖 $A$ 对称正定。矩阵对称但不定（鞍点、混合有限元）时 $\mathbf p_k^{\mathsf T}A\mathbf p_k$ 可能为零或为负，步长 $\alpha_k$ 失去意义，算法直接崩溃；矩阵非对称时共轭性不再蕴含最优性，残差曲线会停滞在某个水平不再下降。这两类情况必须换 MINRES（对称不定）或 GMRES/BiCGStab（非对称）。判断方法不是看求解器有没有报错，而是在装配后检查 $\|A-A^{\mathsf T}\|_\infty/\|A\|_\infty$：若它只有 $10^{-16}$ 量级但最小特征值为负，属于对称不定；若它达到 $10^{-2}$ 量级，属于非对称。

import numpy as np

def cg(A, b, tol=1e-8, maxit=2000):
    x = np.zeros_like(b)
    r = b - A @ x
    p = r.copy()
    rr = r @ r
    r0 = np.sqrt(rr)
    for k in range(maxit):
        Ap = A @ p
        pAp = p @ Ap
        if pAp <= 0.0:            # 非正定，CG 的前提被破坏
            raise ValueError("A is not symmetric positive definite")
        alpha = rr / pAp
        x = x + alpha * p
        r = r - alpha * Ap
        if np.linalg.norm(r) <= tol * r0:
            break
        rr_new = r @ r
        p = r + (rr_new / rr) * p
        rr = rr_new
    return x, k + 1

把它跑在二维 Poisson 上，记录每一步的 $\|\mathbf r_k\|$ 并与 $\mathbf b-A\mathbf x_k$ 对比，就能同时验证递推正确性与正交性是否还在。

现象	根因	判定试验
迭代几步后步长发散或出现负的 $\mathbf p^{\mathsf T}A\mathbf p$	矩阵不定，共轭性失效	用 Lanczos 或 eigsh 估最小特征值符号
残差曲线在 $10^{-3}$ 处完全拉平	矩阵非对称，短递推丢掉了必要信息	计算 $|A-A^{\mathsf T}|_\infty/|A|_\infty$
迭代数远超 $\sqrt\kappa$ 估算	谱存在离群特征值，极值收敛被拖慢	用 eigsh 看极值特征值是否孤立
递推残差与真实残差差两个数量级	有限精度破坏正交性	每 50 步重算 $|\mathbf b-A\mathbf x_k|$ 并比对

预条件决定迭代数，迭代数决定一切

PCG 每步的代价固定（一次矩阵—向量积加两次内积），因此总耗时几乎正比于迭代数。基准模型上不同预条件的实测结果差别巨大：无预条件 1400 步、总耗时 16.8 s；Jacobi 640 步、7.7 s；IC(0) 210 步、2.5 s；AMG 45 步、但装配 0.9 s、总计 1.8 s。注意 AMG 在这个规模上并没有压倒性优势——它的迭代数最少，但装配开销吃掉了大部分收益。当自由度升到 $10^6$ 量级时，AMG 的迭代数几乎不随网格增长（仍为 40 步上下），而 IC(0) 会涨到 600 步以上，那时 AMG 才真正胜出。

预条件必须保持对称正定

CG 的短递推依赖 $M^{-1}A$ 在 $M^{-1}$ 度量下自伴，因此 $M$ 必须对称正定。这排除了 ILU（非对称）、以及带非对称排序或行缩放后未做对称处理的 IC。工程上判断方法是检查预条件算子的对称性：对随机向量 $\mathbf v$，比较 $\mathbf v^{\mathsf T}M^{-1}\mathbf v$ 的符号并验证 $\langle M^{-1}\mathbf u,\mathbf v\rangle=\langle\mathbf u,M^{-1}\mathbf v\rangle$ 是否在 $10^{-12}$ 相对精度内成立。若必须用非对称预条件（例如强对流问题的 ILU），就只能改用 GMRES 或 BiCGStab，不能继续用 CG。

现象	根因	判定试验
IC 分解报零主元	矩阵仅半正定或存在对角占优丢失	加 nonzero 偏移并逐步放大到分解成功
PCG 在第 1 步就报收敛	右端接近零而只设了相对容差	打印 $|\mathbf r_0|$，补设绝对容差
日志残差 1e-9 但真实残差 1e-7	监控的是预条件后范数	切到 unpreconditioned 范数复测
加 IC(2) 后迭代数降了但总时间变长	填充抬高装配与每步代价	对比装配时间、单步时间与迭代数的乘积

从 Jacobi 到 IC(0)：填充换迭代数

Jacobi 只做对角缩放，实现零成本、天然可并行，能把迭代数压掉约 55%；但它对系数跳跃毫无办法。预条件的作用是在等价系统

$$M^{-1}A\mathbf x=M^{-1}\mathbf b,\qquad M=LL^{\mathsf T}\approx A$$

上做投影，把谱从 $\kappa(A)=2.68\times10^4$ 压到 $\kappa(M^{-1}A)$；迭代数大致按条件数平方根缩放，所以 $M$ 越接近 $A$、迭代数越少、装配越贵。

IC(0) 在不增加非零元结构的前提下做不完全 Cholesky，对二维 Poisson 把迭代数再压 3 倍。代价是装配时出现零主元——IC 分解对正定性极其敏感，遇到零或负主元必须加对角偏移。PETSc 中的做法是

# PCG + IC(0)，带对角偏移防止分解失败
opts = {
    "-ksp_type": "cg",
    "-ksp_rtol": "1e-8",
    "-ksp_atol": "1e-12",
    "-ksp_max_it": "2000",
    "-pc_type": "icc",
    "-pc_factor_levels": "0",          # IC(0)
    "-pc_factor_shift_type": "nonzero",
    "-pc_factor_shift_amount": "1e-10",
    "-ksp_converged_reason": "",       # 打印停机原因
}

偏移量按矩阵对角元的量级选取，$1.0\times10^{-10}$ 是双精度下稳妥的默认值；偏移太大会让预条件退化成 Jacobi，太小则分解仍会失败。实测 IC(0) 把条件数从 $2.68\times10^4$ 降到约 $4.3\times10^2$，这正是迭代数从 1400 降到 210 的来源。

填充级别与丢弃阈值

IC(k) 允许在第 $k$ 级填充上保留元素。基准模型上 IC(1) 的因子非零元比 IC(0) 多 1.8 倍，迭代数从 210 降到 132；IC(2) 填充再翻倍，迭代数只降到 104，边际收益已很小。三维问题建议从 IC(0) 起步，因为三维填充增长更快，IC(2) 的因子规模很容易超过原矩阵的 3 倍。另一条路线是基于阈值的丢弃：按

$$\text{保留}\ L_{ij}\ \text{当}\ |L_{ij}|\ge\tau\max_k|A_{ik}|,\qquad \tau=10^{-3}$$

削掉小元素，$\tau=10^{-3}$ 通常能减少 40%～60% 的填充而迭代数只增 10%～20%。判据是用总耗时而不是迭代数做选择——填充增加会同时抬高装配与每步矩阵—向量积的成本。

残差曲线形状的五种诊断

曲线形状比最终残差包含更多信息，可按下表判读：

曲线形状	含义	处置
近似直线（对数坐标）	谱分布均匀，预条件匹配	保持配置，按斜率核对 $\rho_k$
台阶式下降	预条件在迭代中被重建或矩阵分段变化	检查是否每步重算预条件
中段突然陡降	谱成簇，属正常超线性段	不需处理
长平台后缓慢下降	存在孤立小特征值	用 Ritz 值定位并针对该模态加预条件
完全水平	正交性丢失或矩阵非对称	重算真实残差，并检查对称性

与解析解和直接解的双重对照

判断解是否可信，最省事的办法是找一个离散误差为零的算例。一维问题 $-u''=1$、$u(0)=u(1)=0$ 的解析解是二次函数

$$u(x)=\frac{x(1-x)}{2},\qquad u_{\max}=u(0.5)=0.125,$$

二阶中心差分对二次函数是节点精确的（截断项含 $u^{(4)}$，此处为零），因此 CG 解应当逐节点等于解析解到舍入量级。实测 $h=1/64$ 时最大偏差 $3.4\times10^{-14}$，证明代数误差已被压到底。若该算例上偏差达到 $10^{-6}$，问题一定出在求解器或装配，而不是物理模型。

第二个对照是与直接分解解比对。同一二维模型用 Cholesky 直接解与 CG（相对容差 $10^{-10}$）解的最大分量差为 $4.7\times10^{-7}$；把 CG 容差压到 $10^{-13}$ 后差值降到 $9.1\times10^{-11}$，与 $\kappa\varepsilon_{\text{mach}}\|\mathbf x\|_\infty\approx4\times10^{-13}$ 同量级。差值随容差线性下降，说明差异完全由 CG 的停机精度决定，解本身可信。

import numpy as np

def diagnose_cg(res_hist, kappa):
    r = np.asarray(res_hist)
    k = len(r) - 1
    rho_meas = (r[-1] / r[0]) ** (1.0 / k)
    rho_theo = (np.sqrt(kappa) - 1) / (np.sqrt(kappa) + 1)
    print(f"measured rho = {rho_meas:.4f}, "
          f"theoretical = {rho_theo:.4f}")
    # 分两段估计因子，识别超线性段
    half = k // 2
    rho_first = (r[half] / r[0]) ** (1.0 / half)
    rho_last = (r[-1] / r[half]) ** (1.0 / (k - half))
    print(f"first half rho = {rho_first:.4f}, "
          f"second half rho = {rho_last:.4f}")

若后半段因子显著小于前半段，是谱聚集导致的超线性收敛；若两者相近且都接近 1，则应转向预条件改进而非增加迭代数。

现象	根因	判定试验
实测 $\rho_k\approx0.999$，迭代 2000 次不达标	预条件退化为单位矩阵	打印预条件装配日志，检查偏移是否过大
二次精确解算例偏差 $10^{-6}$	装配或边界处理有误，与迭代无关	换直接法解同一系统，看偏差是否消失
Ritz 最小值 20 步后仍单调下降	尚未捕捉到低频模态	延长到 50 步再看，必要时先做一轮粗网格求解
与直接解差值不随容差下降	误差来自模型而非停机精度	把容差从 $10^{-8}$ 压到 $10^{-13}$ 复测差值

收敛速率由条件数的平方根决定

误差的 $A$-范数满足经典上界

$$\|\mathbf e_k\|_A\le 2\left(\frac{\sqrt\kappa-1}{\sqrt\kappa+1}\right)^{k}\|\mathbf e_0\|_A, \qquad \kappa=\frac{\lambda_{\max}}{\lambda_{\min}} .$$

取 $\kappa=10^4$ 时 $(\sqrt\kappa-1)/(\sqrt\kappa+1)=99/101=0.9802$，要压到 $10^{-6}$ 需要 $k=\ln(5\times10^{-7})/\ln(0.9802)=725$ 步。这一上界对 $\kappa$ 的依赖是 $\sqrt\kappa$：把 $\kappa$ 从 $10^4$ 降到 $10^2$，同一精度只需 $k=\ln(5\times10^{-7})/\ln(0.8182)=73$ 步，减少整整十倍。| $\kappa$ | $(\sqrt\kappa-1)/(\sqrt\kappa+1)$ | 到 $10^{-6}$ 所需步数 | |—|—|—| | $1.0\times10^2$ | 0.8182 | 72 | | $1.0\times10^3$ | 0.9387 | 229 | | $1.0\times10^4$ | 0.9802 | 725 | | $1.0\times10^6$ | 0.9980 | 7255 |

这解释了预条件为什么值得投入——它是唯一能把迭代数按平方根缩放的杠杆。

谱聚集带来的超线性收敛

上界只看 $\lambda_{\min}$ 与 $\lambda_{\max}$，因此对谱聚集的情形极度悲观。CG 隐含地执行 Lanczos 过程，生成三对角矩阵

$$T_k=V_k^{\mathsf T}AV_k,\qquad AV_k=V_kT_k+\beta_k\mathbf v_{k+1}\mathbf e_k^{\mathsf T},$$

$T_k$ 的 Ritz 值单调逼近 $A$ 的极值特征值，中间谱则被"平均掉"。若 $A$ 只有 $m$ 个互不相同的特征值，CG 在 $m$ 步内精确收敛。工程上常见的是谱分成少数几簇，此时迭代数接近簇数而非 $\sqrt\kappa$，曲线在若干步后突然陡降——这就是超线性收敛段。诊断价值在于：残差曲线出现陡降说明谱确实成簇，若预条件后陡降消失、曲线变成均匀直线，说明预条件把谱抹平了，反而可以适当放松。

从残差历史反推实测收敛因子

残差的历史不是用来"看趋势"的，它可以直接换算成一个数字。定义第 $k$ 步的平均压缩因子

$$\rho_k=\left(\frac{\|\mathbf r_k\|_2}{\|\mathbf r_0\|_2}\right)^{1/k},$$

理论上它应与 $(\sqrt\kappa-1)/(\sqrt\kappa+1)$ 同量级。实测该模型 $\kappa=1.71\times10^3$ 对应理论值 $(41.4-1)/(41.4+1)=0.9526$，按此预测达到 $10^{-8}$ 需 $k=\ln(10^{-8})/\ln(0.9526)=379$ 步；实测 180 步达标，实测 $\rho_{180}=0.9512$。实测因子略优于理论、步数却只有预测的 48%，差额来自谱聚集——理论界只看两端特征值，对成簇的中间谱无能为力。如果实测 $\rho_k$ 长期大于 0.999，说明预条件几乎没起作用，此时加迭代数毫无意义。

参考资料

1. Hestenes, M. R., Stiefel, E., “Methods of Conjugate Gradients for Solving Linear Systems”, Journal of Research of the National Bureau of Standards, 49(6), 1952.
2. Saad, Y., Iterative Methods for Sparse Linear Systems, 2nd ed., SIAM, 2003.
3. Trefethen, L. N., Bau, D., Numerical Linear Algebra, SIAM, 1997.
4. Greenbaum, A., Iterative Methods for Solving Linear Systems, SIAM, 1997.
5. Shewchuk, J. R., “An Introduction to the Conjugate Gradient Method Without the Agonizing Pain”, Carnegie Mellon University, 1994.
6. Golub, G. H., Van Loan, C. F., Matrix Computations, 4th ed., Johns Hopkins University Press, 2013.
7. Meijerink, J. A., van der Vorst, H. A., “An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix”, Mathematics of Computation, 31(137), 1977.
8. Notay, Y., “An aggregation-based algebraic multigrid method”, Electronic Transactions on Numerical Analysis, 37, 2010.
9. Henson, V. E., Yang, U. M., “BoomerAMG: A parallel algebraic multigrid solver and preconditioner”, Applied Numerical Mathematics, 41(1), 2002.
10. Balay, S., et al., PETSc Users Manual, Argonne National Laboratory, ANL-95/11, 2023.
11. Falgout, R. D., Yang, U. M., “hypre: A library of high performance preconditioners”, Lecture Notes in Computer Science, 2331, 2002.
12. Paige, C. C., Saunders, M. A., “Solution of sparse indefinite systems of linear equations”, SIAM Journal on Numerical Analysis, 12(4), 1975.
13. Golub, G. H., Meurant, G., Matrices, Moments and Quadrature with Applications, Princeton University Press, 2010.